(* load the package. "Multipole.m", which must be located in the same directory of this example file. *)
SetDirectory[NotebookDirectory[]];
<<"Multipole.m";
Names["Multipole`*"]

?StevensFactor

?Multipole

?VectorSphericalHarmonicY

?CreateJMatrix

?ReplaceOperator

?ElectricChargeDensity

?MagneticChargeDensity

?PlotColorBarH

?PlotColorBarV

?PlotWaveFunction

Multipole[]

(*
 * example of equivalent operator
 *)
r={x,y,z};
J=CreateJMatrix[4];
O4=ReplaceOperator[Multipole["O40",r]+5Multipole["O44",r],J,r];
O6c=ReplaceOperator[Multipole["O60",r]-21Multipole["O64",r],J,r];
O6t=ReplaceOperator[Multipole["O62",r]-Multipole["O66",r],J,r];
Hcef[W_,x_,y_]:=W(x/60 O4+(1-Abs[x])/1260O6c+ y/30O6t);
Eigenvalues[Hcef[W,x,y]]//MatrixForm

(*
 * example of vector spherical harmonics
 *)
{VectorSphericalHarmonicY[0,0,0,th,ph],
Table[VectorSphericalHarmonicY[1,q,0,th,ph],{q,-1,1}],
VectorSphericalHarmonicY[0,0,1,th,ph]}//FullSimplify

(*
 * color bars
 *)
{PlotColorBarH[],PlotColorBarV[ColorData["GrayTones"]],PlotColorBarH[ColorData["RedGreenSplit"]],PlotColorBarV[ColorData["DeepSeaColors"]]}

(*
 * example of visualization
 *)
uG7p={Sqrt[1/6],0,0,0,-Sqrt[5/6],0};
uG7m={0,-Sqrt[5/6],0,0,0,Sqrt[1/6]};
uG8ap={Sqrt[5/6],0,0,0,Sqrt[1/6],0};
uG8am={0,Sqrt[1/6],0,0,0,Sqrt[5/6]};
uG8bp={0,0,1,0,0,0};
uG8bm={0,0,0,1,0,0};

u5g3p=(uG8ap-I uG8bp)/Sqrt[2];
u5g3m=-(uG8am+I uG8bm)/Sqrt[2];
u2up=(uG8ap-I uG8bp)/Sqrt[2];
u2um=(uG8am-I uG8bm)/Sqrt[2];
u5u=((7+4Sqrt[2]I)^(1/4)/Sqrt[3]uG8ap+(Sqrt[2]-I)/Sqrt[3]uG8am-(-1)^(1/4)uG8bp+uG8bm)/2;
u5up={-((I (Sqrt[5]) (((7-4 I (Sqrt[2])))^(1/4)))/(6 (Sqrt[2]))),((I+(Sqrt[2]))/(6 (Sqrt[2]))),-(1/2) (((-1))^(1/4)),-(1/2),-((I (((7-4 I (Sqrt[2])))^(1/4)))/(6 (Sqrt[2]))),(1/6) (Sqrt[(5/2)]) (I+(Sqrt[2]))};

{
{PlotWaveFunction[1,{uG7p}],PlotWaveFunction[1,{uG7m}]},
{PlotWaveFunction[1,{uG8ap}],PlotWaveFunction[1,{uG8am}],PlotWaveFunction[1,{uG8bp}],PlotWaveFunction[1,{uG8bm}]},
{PlotWaveFunction[1,{uG8ap,uG8am,uG8bp,uG8bm},{0,0,1}]},
{PlotWaveFunction[1,{u5g3p,u5g3m},{0,0,1}]},
{PlotWaveFunction[1,{u2up,u2um},{0,0,1}]},
{PlotWaveFunction[1,{u5u},{0,0,1}]}
}

(*
 * example of Pr Gamma3
 *)
uG3a={Sqrt[7/24],0,0,0,-Sqrt[5/12],0,0,0,Sqrt[7/24]};
uG3b={0,0,1/Sqrt[2],0,0,0,1/Sqrt[2],0,0};

u2up=(uG3a+I uG3b)/Sqrt[2];
u2um=(uG3a-I uG3b)/Sqrt[2];

{
{PlotWaveFunction[2,{uG3a}],PlotWaveFunction[2,{uG3b}]},
{PlotWaveFunction[2,{u2up,u2um},{0,0,1}]},
{PlotWaveFunction[2,{u2up},{0,0,1}]}
}