(* ::Package:: *)

(*
 * Compute equivalent operators from a c-number expression and visualize a wave function
 * Copyright Dec. 2007 (C) Hiroaki Kusunose <kusu@phys.sci.ehime-u.ac.jp>
 * Reference: H. Kusunose: J. Phys. Soc. Jpn. 77 (2008) 064710, http://arxiv.org/abs/0803.4358
 *
 * This version is for Mathematica 7 or later.
 *
 * Usage: see, multipole_example.nb
 * (1) put this package into directories included in the search path.
 * (2) load this package by using "<< Multipole`".
 *
 *)

(* user functions *)
BeginPackage["Multipole`"];

CreateJMatrix::usage =
"CreateJMatrix[J] creates a set of {\!\(\*SubscriptBox[\(J\), \(x\)]\),\!\(\*SubscriptBox[\(J\), \(y\)]\),\!\(\*SubscriptBox[\(J\), \
\(z\)]\)} matrices (basis in descending order of M) for given J.";

ReplaceOperator::usage =
"ReplaceOperator[ex,j,r] replaces a polynomial expression, ex, in terms of r={x,y,z} by a corresponding equivalent operator with \
j={\!\(\*SubscriptBox[\(J\), \(x\)]\),\!\(\*SubscriptBox[\(J\), \(y\)]\),\!\(\*SubscriptBox[\(J\), \(z\)]\)}.";

Multipole::usage =
"Multipole[\"name\",r] gives a polynomial expression in terms of r={x,y,z} for a multipole specified by name. For further details, evaluate Multipole[].";

VectorSphericalHarmonicY::usage =
"VectorSphericalHarmonicY[p,q,l,\[Theta],\[Phi]] gives a vector spherical harmonics, \!\(\*SubsuperscriptBox[OverscriptBox[\(Y\), \(\[RightVector]\)], \(pq\), \(l\)]\)(\[Theta],\[Phi])."\
;

StevensFactor::usage =
"StevensFactor[n,p] gives a generalized Stevens factor of rank p for a Hund-rule ground multiplet in \!\(\*SuperscriptBox[\(f\), \(n\)]\) configuration.";

ElectricChargeDensity::usage =
"ElectricChargeDensity[n,u,\[Theta],\[Phi]] gives an angle dependence of the electric charge density, \!\(\*SubscriptBox[\(\[Rho]\), \
\(electric\)]\)(\[Theta],\[Phi]) for a given wave function, u in the M basis.";

MagneticChargeDensity::usage =
"MagneticChargeDensity[n,u,\[Theta],\[Phi]] gives an angle dependence of the magnetic charge density, \!\(\*SubscriptBox[\(\[Rho]\), \
\(magnetic\)]\)(\[Theta],\[Phi]) for a given wave function, u in the M basis.";

PlotColorBarH::usage =
"PlotColorBarH[gmap] gives a horizontal color bar of gmap.";

PlotColorBarV::usage =
"PlotColorBarV[gmap] gives a vertical color bar of gmap.";

PlotWaveFunction::usage =
"PlotWaveFunction[n,u,vp={2.5,2.5,2.5},gmap=MagneticMap] gives a wavefunction plot for a set of u in the M basis with \
\!\(\*SuperscriptBox[\(f\), \(n\)]\) configuration. A radial extent and a surface color using gmap are determined by \!\(\*SubscriptBox[\(\[Rho]\), \
\(electric\)]\)(\[Theta],\[Phi]) and \!\(\*SubscriptBox[\(\[Rho]\), \(magnetic\)]\)(\[Theta],\[Phi]). vp specifies a view point.";


(* detailed definitions *)
Begin["`Private`"];

CreateJMatrix[J_] := Module[{Jp, Jm, Jz, Jvec, i}, Jvec = Table[i, {i, J, -J, -1}]; Jp[j1_, j2_] := KroneckerDelta[j1, j2 + 1]*Sqrt[(J - j2)*(J + j2 + 1)]; 
     Jm[j1_, j2_] := KroneckerDelta[j1, j2 - 1]*Sqrt[(J + j2)*(J - j2 + 1)]; Jz[j1_, j2_] := KroneckerDelta[j1, j2]*j2; 
     {(1/2)*(Outer[Jp, Jvec, Jvec] + Outer[Jm, Jvec, Jvec]), (Outer[Jp, Jvec, Jvec] - Outer[Jm, Jvec, Jvec])/(2*I), Outer[Jz, Jvec, Jvec]}
];

ReplaceOperator[ex_, {jx_, jy_, jz_}, {x_, y_, z_}] := Module[{Lbar, GetExp, ReplaceTerm, term, j1, j2, j3, f}, 
    Lbar[L_] := Plus @@ Apply[Dot, Permutations[L], 1]; GetExp[(c_)?NumericQ*(val_)] := {c, Exponent[val, x], Exponent[val, y], Exponent[val, z]}; 
     GetExp[val_] := {1, Exponent[val, x], Exponent[val, y], Exponent[val, z]}; 
     ReplaceTerm[et_] := Module[{c, l, m, n, ls}, {c, l, m, n} = GetExp[et]; ls = Join[Table[j1, {l}], Table[j2, {m}], Table[j3, {n}]]; 
        (c*l!*m!*n!*Lbar[ls])/(l + m + n)!]; f[ex1_] := ReplaceTerm[ex1]; f[(ex1_) + (ex2_)] := ReplaceTerm[ex1] + f[ex2]; 
     f[Expand[ex]] /. {j1 -> jx, j2 -> jy, j3 -> jz}
];

VectorSphericalHarmonicY[p_, q_, l_, th_, ph_] := Module[{ev, m, vsh}, ev = {(-Sqrt[2]^(-1))*{1, I, 0}, {0, 0, 1}, (1/Sqrt[2])*{1, -I, 0}}; 
     vsh = (-1)^(l + q + 1)*Sqrt[2*p + 1]*Sum[ThreeJSymbol[{p, -q}, {l, q - m}, {1, m}]*SphericalHarmonicY[l, q - m, th, ph]*ev[[2 - m]], 
        {m, Max[q - l, -1], Min[q + l, 1]}]; vsh
];

StevensFactor[n_, rank_] := gpn[[n,rank]];

ElectricChargeDensity[n_, u_, th_, phi_] := Module[{p, q, rho, J, x, y, z, r}, J = CreateJMatrix[nj[[n]]]; r = {x, y, z}; 
     rho = n + Sum[(2*p + 1)*StevensFactor[n, p]*Sum[Re[Conjugate[u] . ReplaceOperator[Multipole[multipoles[[p,q]], r], J, r] . u]*
           Multipole[multipoles[[p,q]], r], {q, Length[multipoles[[p]]]}], {p, 2, 6, 2}]; rho /. {x -> Sin[th]*Cos[phi], y -> Sin[th]*Sin[phi], z -> Cos[th]}
];

MagneticChargeDensity[n_, u_, th_, phi_] := Module[{p, q, rho, J, x, y, z, r}, J = CreateJMatrix[nj[[n]]]; r = {x, y, z}; 
     rho = Sum[(2*p + 1)*StevensFactor[n, p]*Sum[Re[Conjugate[u] . ReplaceOperator[Multipole[multipoles[[p,q]], r], J, r] . u]*Multipole[multipoles[[p,q]], r], 
         {q, Length[multipoles[[p]]]}], {p, 1, 5, 2}]; rho /. {x -> Sin[th]*Cos[phi], y -> Sin[th]*Sin[phi], z -> Cos[th]}
];

NElectricChargeDensity[n_, u_, th_, phi_] := Module[{p, q, rho, J, x, y, z, r}, J = N[CreateJMatrix[nj[[n]]]]; r = {x, y, z}; 
     rho = n + Sum[(2*p + 1)*StevensFactor[n, p]*Sum[Chop[Re[Conjugate[N[u]] . ReplaceOperator[Multipole[multipoles[[p,q]], r], J, r] . N[u]]]*
           Multipole[multipoles[[p,q]], r], {q, Length[multipoles[[p]]]}], {p, 2, 6, 2}]; rho /. {x -> Sin[th]*Cos[phi], y -> Sin[th]*Sin[phi], z -> Cos[th]}
];

NMagneticChargeDensity[n_, u_, th_, phi_] := Module[{p, q, rho, J, x, y, z, r}, J = N[CreateJMatrix[nj[[n]]]]; r = {x, y, z}; 
     rho = Sum[(2*p + 1)*StevensFactor[n, p]*Sum[Chop[Re[Conjugate[N[u]] . ReplaceOperator[Multipole[multipoles[[p,q]], r], J, r] . N[u]]]*Multipole[multipoles[[p,q]], r], 
         {q, Length[multipoles[[p]]]}], {p, 1, 5, 2}]; rho /. {x -> Sin[th]*Cos[phi], y -> Sin[th]*Sin[phi], z -> Cos[th]}
];

MagneticMap[x_] := If[x < 0.5, Hue[0.65, 1 - 2*x, 1], Hue[0, 2*x - 1, 1]];

PlotColorBarH[gmap_:MagneticMap] := Module[{x, y}, DensityPlot[x, {x, 0, 1}, {y, 0, 1}, ColorFunction -> (gmap[#1] & ), PlotPoints -> 50, AspectRatio -> 0.12, 
     Mesh -> False, Frame -> False]
];

PlotColorBarV[gmap_:MagneticMap] := Module[{x, y}, DensityPlot[y, {x, 0, 1}, {y, 0, 1}, ColorFunction -> (gmap[#1] & ), PlotPoints -> 50, AspectRatio -> 1/0.12, 
     Mesh -> False, Frame -> False]
];

SphericalContourPlot3D[fe_, fm_, th_, phi_, vp_, gmap_] := Module[{pn = 35, data, fmmax, ScaledMap, x, y, z, u, v, f},
   data = Table[N[fm], {th, 0, Pi, Pi/pn}, {phi, 0, 2*Pi, 2*(Pi/pn)}];
   fmmax = Max[Abs[data]];
   ScaledMap[f_] := If[fmmax < 10^(-6), 0.5, Rescale[f,{-fmmax,fmmax}]];
   SphericalPlot3D[fe, th, phi,
      Boxed -> False, Axes -> False,
      Mesh -> pn, PlotPoints -> pn,
      ImageSize -> 400,
      ViewPoint -> vp,
      Lighting -> {{"Ambient", White}},
      ColorFunction -> Function[{x, y, z, u, v, f}, gmap[ScaledMap[fm /. {th -> Pi*u, phi -> 2*Pi*v}]]]]
];

PlotWaveFunction[n_, u_, vp_:{2.5, 2.5, 2.5}, gmap_:MagneticMap] := Module[{i, th, phi, nu, rhoe, rhom},
   nu = Length[u];
   rhoe = Sum[NElectricChargeDensity[n, N[u[[i]]], th, phi], {i, nu}];
   rhom = Sum[NMagneticChargeDensity[n, N[u[[i]]], th, phi], {i, nu}];
   SphericalContourPlot3D[rhoe, rhom, th, phi, vp, gmap]
];

Multipole[___] := Module[{}, Print["The following multipoles are available."]; Print["(i) cartesian multipoles:"]; Print[multipoles]; 
     Print["(ii) cubic multipoles (rank,Gamma,multiplicity,gamma):"]; Print[cubicMultipoles]; Print["(iii) stevens operators:"]; Print[stevensMultipoles];
];

multipoles = {{"c11", "s11", "j10"}, {"c21", "c22", "s21", "s22", "j20"}, {"c31", "c32", "c33", "s31", "s32", "s33", "j30"}, 
    {"c41", "c42", "c43", "c44", "s41", "s42", "s43", "s44", "j40"}, {"c51", "c52", "c53", "c54", "c55", "s51", "s52", "s53", "s54", "s55", "j50"}, 
    {"c61", "c62", "c63", "c64", "c65", "c66", "s61", "s62", "s63", "s64", "s65", "s66", "j60"}};

Multipole["c11", {x_, y_, z_}] := x; 
Multipole["s11", {x_, y_, z_}] := y; 
Multipole["j10", {x_, y_, z_}] := z; 
Multipole["c21", {x_, y_, z_}] := Sqrt[3]*x*z; 
Multipole["c22", {x_, y_, z_}] := (1/2)*Sqrt[3]*(x^2 - y^2); 
Multipole["s21", {x_, y_, z_}] := Sqrt[3]*y*z; 
Multipole["s22", {x_, y_, z_}] := Sqrt[3]*x*y; 
Multipole["j20", {x_, y_, z_}] := (1/2)*(-x^2 - y^2 + 2*z^2); 
Multipole["c31", {x_, y_, z_}] := (-(1/2))*Sqrt[3/2]*x*(x^2 + y^2 - 4*z^2); 
Multipole["c32", {x_, y_, z_}] := (1/2)*Sqrt[15]*(x^2 - y^2)*z; 
Multipole["c33", {x_, y_, z_}] := (1/2)*Sqrt[5/2]*(x^3 - 3*x*y^2); 
Multipole["s31", {x_, y_, z_}] := (-(1/2))*Sqrt[3/2]*y*(x^2 + y^2 - 4*z^2); 
Multipole["s32", {x_, y_, z_}] := Sqrt[15]*x*y*z; 
Multipole["s33", {x_, y_, z_}] := (1/2)*Sqrt[5/2]*(3*x^2*y - y^3); 
Multipole["j30", {x_, y_, z_}] := -((3*x^2*z)/2) - (3*y^2*z)/2 + z^3; 
Multipole["c41", {x_, y_, z_}] := (-(1/2))*Sqrt[5/2]*x*z*(3*x^2 + 3*y^2 - 4*z^2); 
Multipole["c42", {x_, y_, z_}] := (-(1/4))*Sqrt[5]*(x^2 - y^2)*(x^2 + y^2 - 6*z^2); 
Multipole["c43", {x_, y_, z_}] := (1/2)*Sqrt[35/2]*(x^3 - 3*x*y^2)*z; 
Multipole["c44", {x_, y_, z_}] := (1/8)*Sqrt[35]*(x^4 - 6*x^2*y^2 + y^4); 
Multipole["s41", {x_, y_, z_}] := (-(1/2))*Sqrt[5/2]*y*z*(3*x^2 + 3*y^2 - 4*z^2); 
Multipole["s42", {x_, y_, z_}] := (-(1/2))*Sqrt[5]*x*y*(x^2 + y^2 - 6*z^2); 
Multipole["s43", {x_, y_, z_}] := (1/2)*Sqrt[35/2]*(3*x^2*y - y^3)*z; 
Multipole["s44", {x_, y_, z_}] := (1/2)*Sqrt[35]*x*y*(x^2 - y^2); 
Multipole["j40", {x_, y_, z_}] := (1/8)*(35*z^4 - 30*z^2*(x^2 + y^2 + z^2) + 3*(x^2 + y^2 + z^2)^2); 
Multipole["c51", {x_, y_, z_}] := (1/8)*Sqrt[15]*x*(x^4 + y^4 - 12*y^2*z^2 + 8*z^4 + 2*x^2*(y^2 - 6*z^2)); 
Multipole["c52", {x_, y_, z_}] := (-(1/4))*Sqrt[105]*(x^2 - y^2)*z*(x^2 + y^2 - 2*z^2); 
Multipole["c53", {x_, y_, z_}] := (-(1/8))*Sqrt[35/2]*(x^3 - 3*x*y^2)*(x^2 + y^2 - 8*z^2); 
Multipole["c54", {x_, y_, z_}] := (3/8)*Sqrt[35]*(x^4 - 6*x^2*y^2 + y^4)*z; 
Multipole["c55", {x_, y_, z_}] := (3/8)*Sqrt[7/2]*(x^5 - 10*x^3*y^2 + 5*x*y^4); 
Multipole["s51", {x_, y_, z_}] := (1/8)*Sqrt[15]*y*(x^4 + y^4 - 12*y^2*z^2 + 8*z^4 + 2*x^2*(y^2 - 6*z^2)); 
Multipole["s52", {x_, y_, z_}] := (-(1/2))*Sqrt[105]*x*y*z*(x^2 + y^2 - 2*z^2); 
Multipole["s53", {x_, y_, z_}] := (-(1/8))*Sqrt[35/2]*(3*x^2*y - y^3)*(x^2 + y^2 - 8*z^2); 
Multipole["s54", {x_, y_, z_}] := (3/2)*Sqrt[35]*x*y*(x^2 - y^2)*z; 
Multipole["s55", {x_, y_, z_}] := (3/8)*Sqrt[7/2]*(5*x^4*y - 10*x^2*y^3 + y^5); 
Multipole["j50", {x_, y_, z_}] := (1/8)*(63*z^5 - 70*z^3*(x^2 + y^2 + z^2) + 15*z*(x^2 + y^2 + z^2)^2); 
Multipole["c61", {x_, y_, z_}] := (1/8)*Sqrt[21]*x*z*(5*x^4 + 5*y^4 - 20*y^2*z^2 + 8*z^4 + 10*x^2*(y^2 - 2*z^2)); 
Multipole["c62", {x_, y_, z_}] := (1/16)*Sqrt[105/2]*(x^2 - y^2)*(x^4 + y^4 - 16*y^2*z^2 + 16*z^4 + 2*x^2*(y^2 - 8*z^2)); 
Multipole["c63", {x_, y_, z_}] := (-(1/8))*Sqrt[105/2]*(x^3 - 3*x*y^2)*z*(3*x^2 + 3*y^2 - 8*z^2); 
Multipole["c64", {x_, y_, z_}] := (-(3/16))*Sqrt[7]*(x^4 - 6*x^2*y^2 + y^4)*(x^2 + y^2 - 10*z^2); 
Multipole["c65", {x_, y_, z_}] := (3/8)*Sqrt[77/2]*(x^5 - 10*x^3*y^2 + 5*x*y^4)*z; 
Multipole["c66", {x_, y_, z_}] := (1/16)*Sqrt[231/2]*(x^6 - 15*x^4*y^2 + 15*x^2*y^4 - y^6); 
Multipole["s61", {x_, y_, z_}] := (1/8)*Sqrt[21]*y*z*(5*x^4 + 5*y^4 - 20*y^2*z^2 + 8*z^4 + 10*x^2*(y^2 - 2*z^2)); 
Multipole["s62", {x_, y_, z_}] := (1/8)*Sqrt[105/2]*x*y*(x^4 + y^4 - 16*y^2*z^2 + 16*z^4 + 2*x^2*(y^2 - 8*z^2)); 
Multipole["s63", {x_, y_, z_}] := (-(1/8))*Sqrt[105/2]*(3*x^2*y - y^3)*z*(3*x^2 + 3*y^2 - 8*z^2); 
Multipole["s64", {x_, y_, z_}] := (-(3/4))*Sqrt[7]*x*y*(x^2 - y^2)*(x^2 + y^2 - 10*z^2); 
Multipole["s65", {x_, y_, z_}] := (3/8)*Sqrt[77/2]*(5*x^4*y - 10*x^2*y^3 + y^5)*z; 
Multipole["s66", {x_, y_, z_}] := (1/8)*Sqrt[231/2]*x*y*(3*x^4 - 10*x^2*y^2 + 3*y^4); 
Multipole["j60", {x_, y_, z_}] := (1/16)*(231*z^6 - 315*z^4*(x^2 + y^2 + z^2) + 105*z^2*(x^2 + y^2 + z^2)^2 - 5*(x^2 + y^2 + z^2)^3);

cubicMultipoles = {{"x1411", "x1412", "x1413"}, {"x2311", "x2312", "x2511", "x2512", "x2513"}, {"x3211", "x3411", "x3412", "x3413", "x3511", "x3512", "x3513"}, 
    {"x4111", "x4311", "x4312", "x4411", "x4412", "x4413", "x4511", "x4512", "x4513"}, {"x5311", "x5312", "x5411", "x5412", "x5413", "x5421", "x5422", "x5423", 
     "x5511", "x5512", "x5513"}, {"x6111", "x6211", "x6311", "x6312", "x6411", "x6412", "x6413", "x6511", "x6512", "x6513", "x6521", "x6522", "x6523"}}; 
Multipole["x1411", {x_, y_, z_}] := x; 
Multipole["x1412", {x_, y_, z_}] := y; 
Multipole["x1413", {x_, y_, z_}] := z; 
Multipole["x2311", {x_, y_, z_}] := (1/2)*(-x^2 - y^2 + 2*z^2); 
Multipole["x2312", {x_, y_, z_}] := (1/2)*Sqrt[3]*(x^2 - y^2); 
Multipole["x2511", {x_, y_, z_}] := Sqrt[3]*y*z; 
Multipole["x2512", {x_, y_, z_}] := Sqrt[3]*x*z; 
Multipole["x2513", {x_, y_, z_}] := Sqrt[3]*x*y; 
Multipole["x3211", {x_, y_, z_}] := Sqrt[15]*x*y*z; 
Multipole["x3411", {x_, y_, z_}] := x^3 - (3/2)*x*(y^2 + z^2); 
Multipole["x3412", {x_, y_, z_}] := -((3*x^2*y)/2) + y^3 - (3*y*z^2)/2; 
Multipole["x3413", {x_, y_, z_}] := -((3*x^2*z)/2) - (3*y^2*z)/2 + z^3; 
Multipole["x3511", {x_, y_, z_}] := (1/2)*Sqrt[15]*x*(y^2 - z^2); 
Multipole["x3512", {x_, y_, z_}] := (1/2)*Sqrt[15]*y*(-x^2 + z^2); 
Multipole["x3513", {x_, y_, z_}] := (1/2)*Sqrt[15]*(x^2 - y^2)*z; 
Multipole["x4111", {x_, y_, z_}] := (1/2)*Sqrt[7/3]*(x^4 + y^4 - 3*y^2*z^2 + z^4 - 3*x^2*(y^2 + z^2)); 
Multipole["x4311", {x_, y_, z_}] := (-(1/4))*Sqrt[5/3]*(x^4 + y^4 + 6*y^2*z^2 - 2*z^4 + 6*x^2*(-2*y^2 + z^2)); 
Multipole["x4312", {x_, y_, z_}] := (1/4)*Sqrt[5]*(x^2 - y^2)*(x^2 + y^2 - 6*z^2); 
Multipole["x4411", {x_, y_, z_}] := (1/2)*Sqrt[35]*y*z*(y^2 - z^2); 
Multipole["x4412", {x_, y_, z_}] := (1/2)*Sqrt[35]*x*z*(-x^2 + z^2); 
Multipole["x4413", {x_, y_, z_}] := (1/2)*Sqrt[35]*x*y*(x^2 - y^2); 
Multipole["x4511", {x_, y_, z_}] := (-(1/2))*Sqrt[5]*y*z*(-6*x^2 + y^2 + z^2); 
Multipole["x4512", {x_, y_, z_}] := (-(1/2))*Sqrt[5]*x*z*(x^2 - 6*y^2 + z^2); 
Multipole["x4513", {x_, y_, z_}] := (-(1/2))*Sqrt[5]*x*y*(x^2 + y^2 - 6*z^2); 
Multipole["x5311", {x_, y_, z_}] := (3/2)*Sqrt[35]*x*y*(x^2 - y^2)*z; 
Multipole["x5312", {x_, y_, z_}] := (1/2)*Sqrt[105]*x*y*z*(x^2 + y^2 - 2*z^2); 
Multipole["x5411", {x_, y_, z_}] := (1/8)*x*(8*x^4 - 40*x^2*(y^2 + z^2) + 15*(y^2 + z^2)^2); 
Multipole["x5412", {x_, y_, z_}] := (1/8)*y*(15*x^4 + 8*y^4 - 40*y^2*z^2 + 15*z^4 + x^2*(-40*y^2 + 30*z^2)); 
Multipole["x5413", {x_, y_, z_}] := (1/8)*z*(15*x^4 + 15*y^4 - 40*y^2*z^2 + 8*z^4 + 10*x^2*(3*y^2 - 4*z^2)); 
Multipole["x5421", {x_, y_, z_}] := (3/8)*Sqrt[35]*x*(y^4 - 6*y^2*z^2 + z^4); 
Multipole["x5422", {x_, y_, z_}] := (3/8)*Sqrt[35]*y*(x^4 - 6*x^2*z^2 + z^4); 
Multipole["x5423", {x_, y_, z_}] := (3/8)*Sqrt[35]*(x^4 - 6*x^2*y^2 + y^4)*z; 
Multipole["x5511", {x_, y_, z_}] := (1/4)*Sqrt[105]*x*(-y^4 + z^4 + 2*x^2*(y^2 - z^2)); 
Multipole["x5512", {x_, y_, z_}] := (1/4)*Sqrt[105]*y*(x^2 - z^2)*(x^2 - 2*y^2 + z^2); 
Multipole["x5513", {x_, y_, z_}] := (-(1/4))*Sqrt[105]*(x^2 - y^2)*z*(x^2 + y^2 - 2*z^2); 
Multipole["x6111", {x_, y_, z_}] := (2*x^6 + 2*y^6 - 15*y^4*z^2 - 15*y^2*z^4 + 2*z^6 - 15*x^4*(y^2 + z^2) - 15*x^2*(y^4 - 12*y^2*z^2 + z^4))/(4*Sqrt[2]); 
Multipole["x6211", {x_, y_, z_}] := (1/4)*Sqrt[1155/2]*(x^2 - y^2)*(x^2 - z^2)*(y^2 - z^2); 
Multipole["x6311", {x_, y_, z_}] := (-(1/4))*Sqrt[7/2]*(x^6 + y^6 - 15*x^4*z^2 - 15*y^4*z^2 + 15*x^2*z^4 + 15*y^2*z^4 - 2*z^6); 
Multipole["x6312", {x_, y_, z_}] := (1/4)*Sqrt[21/2]*(x^2 - y^2)*(x^4 + y^4 - 5*y^2*z^2 + 5*z^4 - x^2*(9*y^2 + 5*z^2)); 
Multipole["x6411", {x_, y_, z_}] := (-(3/4))*Sqrt[7]*y*z*(y^2 - z^2)*(-10*x^2 + y^2 + z^2); 
Multipole["x6412", {x_, y_, z_}] := (3/4)*Sqrt[7]*x*z*(x^2 - z^2)*(x^2 - 10*y^2 + z^2); 
Multipole["x6413", {x_, y_, z_}] := (-(3/4))*Sqrt[7]*x*y*(x^2 - y^2)*(x^2 + y^2 - 10*z^2); 
Multipole["x6511", {x_, y_, z_}] := (1/8)*Sqrt[231/2]*y*z*(3*y^4 - 10*y^2*z^2 + 3*z^4); 
Multipole["x6512", {x_, y_, z_}] := (1/8)*Sqrt[231/2]*x*z*(3*x^4 - 10*x^2*z^2 + 3*z^4); 
Multipole["x6513", {x_, y_, z_}] := (1/8)*Sqrt[231/2]*x*y*(3*x^4 - 10*x^2*y^2 + 3*y^4); 
Multipole["x6521", {x_, y_, z_}] := (1/8)*Sqrt[105/2]*y*z*(16*x^4 - 16*x^2*(y^2 + z^2) + (y^2 + z^2)^2); 
Multipole["x6522", {x_, y_, z_}] := (1/8)*Sqrt[105/2]*x*z*(x^4 + 16*y^4 - 16*y^2*z^2 + z^4 + 2*x^2*(-8*y^2 + z^2)); 
Multipole["x6523", {x_, y_, z_}] := (1/8)*Sqrt[105/2]*x*y*(x^4 + y^4 - 16*y^2*z^2 + 16*z^4 + 2*x^2*(y^2 - 8*z^2)); 
stevensMultipoles = {{"O20", "O22"}, {"O40", "O42", "O43", "O43s", "O44", "O44s"}, {"O60", "O62", "O63", "O64", "O66"}}; 
Multipole["O20", {x_, y_, z_}] := -x^2 - y^2 + 2*z^2; 
Multipole["O22", {x_, y_, z_}] := x^2 - y^2; 
Multipole["O40", {x_, y_, z_}] := 3*x^4 + 3*y^4 - 24*y^2*z^2 + 8*z^4 + 6*x^2*(y^2 - 4*z^2); 
Multipole["O42", {x_, y_, z_}] := (-(x^2 - y^2))*(x^2 + y^2 - 6*z^2); 
Multipole["O43", {x_, y_, z_}] := x*(x^2 - 3*y^2)*z; 
Multipole["O43s", {x_, y_, z_}] := (-y)*(-3*x^2 + y^2)*z; 
Multipole["O44", {x_, y_, z_}] := x^4 - 6*x^2*y^2 + y^4; 
Multipole["O44s", {x_, y_, z_}] := 4*x*y*(x^2 - y^2); 
Multipole["O60", {x_, y_, z_}] := -5*x^6 - 5*y^6 + 90*y^4*z^2 - 120*y^2*z^4 + 16*z^6 - 15*x^4*(y^2 - 6*z^2) - 15*x^2*(y^4 - 12*y^2*z^2 + 8*z^4); 
Multipole["O62", {x_, y_, z_}] := (x^2 - y^2)*(x^4 + y^4 - 16*y^2*z^2 + 16*z^4 + 2*x^2*(y^2 - 8*z^2)); 
Multipole["O63", {x_, y_, z_}] := (-x)*(x^2 - 3*y^2)*z*(3*x^2 + 3*y^2 - 8*z^2); 
Multipole["O64", {x_, y_, z_}] := (-(x^4 - 6*x^2*y^2 + y^4))*(x^2 + y^2 - 10*z^2); 
Multipole["O66", {x_, y_, z_}] := x^6 - 15*x^4*y^2 + 15*x^2*y^4 - y^6;

gpn = {{6/7, -(2/35), -(2/35), 2/315, 4/693, 0}, {4/5, -(52/2475), 26/2475, -(4/5445), -(4/7623), 272/4459455}, 
    {8/11, -(7/1089), 70/4719, -(136/467181), 680/1401543, -(1615/42513471)}, {3/5, 14/1815, 7/605, 952/2335905, 68/467181, 2584/42513471}, 
    {2/7, 13/315, -(26/945), 26/10395, -(52/7623), 0}, {0, 0, 0, 0, 0, 0}, {2, 0, 0, 0, 0, 0}, {3/2, -(1/99), -(7/660), 2/16335, 13/78408, -(1/891891)}, 
    {4/3, -(2/315), -(4/2457), -(8/135135), -(188/2675673), 4/3864861}, {5/4, -(1/450), 1/500, -(1/30030), 1/135135, -(5/3864861)}, 
    {6/5, 4/1575, 244/75075, 2/45045, 124/3270267, 8/3864861}, {7/6, 1/99, 1/990, 8/49005, -(7/58806), -(5/891891)}, 
    {8/7, 2/63, -(16/455), -(2/1155), 40/27027, 4/27027}};

nj = {5/2, 4, 9/2, 4, 5/2, 0, 7/2, 6, 15/2, 8, 15/2, 6, 7/2};

End[];

Protect[
   CreateJMatrix, ReplaceOperator, VectorSphericalHarmonicY, StevensFactor, ElectricChargeDensity, MagneticChargeDensity, PlotColorBarH, PlotColorBarV, 
   PlotWaveFunction, Multipole
];

EndPackage[];