For[i=0,i<50,i++,
uu[i]= u /. FindRoot[u-Sin[u]/Sqrt[2] == .04 Pi i,{u, .04 Pi i}] ;
pp[i]=N[2 ArcTan[ Cos[uu[i]/2], Sqrt[(Sqrt[2]+1)/(Sqrt[2]-1)] Sin[uu[i]/2]]] ]
xy=Table[{.5 Cos[pp[i]]/(1+Cos[pp[i]]/Sqrt[2.]),.5 Sin[pp[i]]/(1+Cos[pp[i]]/Sqrt[2.])},{i,0,49}]
ListPlot[%,AspectRatio->Automatic,PlotStyle->{RGBColor[1,0,0]}]
vxy=Table[{-Sin[pp[i]],Cos[pp[i]]+1/Sqrt[2.]},{i,0,49}]
ListPlot[%,AspectRatio->Automatic,PlotStyle->{RGBColor[1,0,0]}]
---
e=N[1/Sqrt[2]]
phi[u_]=2 ArcTan[ Cos[u/2], Sqrt[(1+e)/(1-e)] Sin[u/2]]
r[u_]=(1-e^2)/(1+e Cos[phi[u]])
v[u_]=Sqrt[1+e^2+2 e Cos[phi[u]]]
vr[u_]=Sqrt[v[u]^2-(1-e^2)^2/r[u]^2]Sign[Sin[u]]
wt[u_]=25(u-e Sin[u])/Pi

ParametricPlot[{wt[u],r[u]},{u,0,4 Pi},PlotRange->{0,2},PlotPoints->60]
ParametricPlot[{wt[u],v[u]},{u,0,4 Pi},PlotRange->{0,2},PlotPoints->60]
ParametricPlot[{wt[u],vr[u]},{u,0,4 Pi},PlotPoints->60]

ParametricPlot[{wt[u],r[u]},{u,-uu[5],uu[5]},PlotRange->{.2,.85},PlotPoints->60]
ParametricPlot[{wt[u],v[u]},{u,-uu[5],uu[5]},PlotRange->{.8,1.8},PlotPoints->60]
ParametricPlot[{wt[u],vr[u]},{u,-uu[5],uu[5]},PlotRange->{-.8,.8},PlotPoints->60]
---
Plot[-1/x +.25/x^2,{x,.1,2.5},PlotPoints->100,PlotStyle->{RGBColor[1,0,0]}]
Plot[-.5,{x,0,2.5},PlotStyle->{RGBColor[1,0,1]}]
Show[%,%%,PlotRange->{-1.25,.25}]
Show[%%%,%%,PlotRange->{-1.25,.25}]
----
Plot[{-.5,-.5/4,-.5/9,-.5/16,-.5/25},{x,-.15,.15}]
Plot[{-.5/4,-.5/9,-.5/16,-.5/25},{x,.85,1.15}]
Plot[{-.5/9,-.5/16,-.5/25},{x,1.85,2.15}]
Plot[{-.5/16,-.5/25},{x,2.85,3.15}]
Plot[{-.5/25},{x,3.85,4.15}]
Show[%%%%%,%%%%,%%%,%%,%,PlotRange->{-.6,.1},Ticks->{None,Automatic}]
---
pr[n_,l_,r_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]]) (2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n] 
v[l_,r_]=l(l+1)/(2 r^2)-1/r
Plot[{pr[1,0,x]-.5,pr[2,0,x]-.125,pr[3,0,x]-.5/9,pr[4,0,x]-.5/16,-.5,-.125,-.5/9,-.5/16,
v[0,x]},{x,0,20},PlotRange->{-1,.5}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Integrate[pr[1,0,r] pr[1,0,r] r^2, {r,0,Infinity}]
Integrate[pr[2,0,r] pr[2,0,r] r^2, {r,0,Infinity}]
Integrate[pr[3,0,r] pr[3,0,r] r^2, {r,0,Infinity}]
Integrate[pr[2,1,r] pr[2,1,r] r^2, {r,0,Infinity}]
Integrate[pr[3,1,r] pr[3,1,r] r^2, {r,0,Infinity}]
Integrate[pr[3,2,r] pr[3,2,r] r^2, {r,0,Infinity}]
Integrate[pr[3,1,r] pr[2,1,r] r^2, {r,0,Infinity}]

Plot[{-.5,-.125,-.5/9,-.5/16,
v[0,x]},{x,0,20},PlotRange->{-1,.1}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{-.125,-.5/9,-.5/16,
v[1,x]},{x,0,20},PlotRange->{-1,.1}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{-.5/16,
v[3,x]},{x,0,20},PlotRange->{-1,.1}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[pr[1,0,x],{x,0.01,5},PlotRange->{-.1,2}]
Plot[pr[4,0,x],{x,0.01,50},PlotRange->{-.04,.25}]
Plot[pr[3,0,x],{x,0.01,30},PlotRange->{-.06,.4}]
Plot[pr[2,0,x],{x,0.01,15},PlotRange->{-.1,.7}]

Plot[pr[4,1,x],{x,0.01,50}]
Plot[pr[3,1,x],{x,0.01,30}]
Plot[pr[2,1,x],{x,0.01,15}]

Plot[pr[4,3,x],{x,0.01,50}]

Plot[x pr[30,10,x],{x,0.01,2000}]
Plot[pr[30,10,x],{x,0.01,2000},PlotRange->{-.00015,.0003}]
Plot[{v[10,x],-.5/30^2},{x,0.01,2000},PlotRange->{-.005,.001},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,1]}}]
----
p[n_,l_,m_,r_,theta_,phi_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]])(2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n]  SphericalHarmonicY[l,m,theta,phi] 

ContourPlot[Abs[p[5,2,2,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-60,60},{z,-60,60},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00001}]
Plot[Abs[p[5,2,2,r,Pi/2,0]]^2,{r,0,60},PlotPoints->60,PlotRange->{0,.00008}]
ContourPlot[Abs[p[5,2,2,Sqrt[x^2+y^2],Pi/2,0]]^2,{x,-60,60},{y,-60,60},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00001}]

ContourPlot[Abs[p[5,2,0,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-60,60},{z,-60,60},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00004}]
Plot[Abs[p[5,2,0,r,0,0]]^2,{r,0,60},PlotPoints->60,PlotRange->{0,.0002}]

ContourPlot[Abs[p[5,2,1,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-60,60},{z,-60,60},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00002}]

ContourPlot[Abs[p[4,1,1,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-35,35},{z,-35,35},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00005}]
ContourPlot[Abs[p[4,1,1,Sqrt[x^2+y^2],Pi/2,0]]^2,{x,-35,35},{y,-35,35},PlotPoints->100,
ContourLines->False,PlotRange->{0,.00005}]
Plot[Abs[p[4,1,1,r,Pi/2,0]]^2,{r,0,35},PlotPoints->60,PlotRange->{0,.0004}]

ContourPlot[Abs[p[4,1,0,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-35,35},{z,-35,35},PlotPoints->100,
ContourLines->False,PlotRange->{0,.0001}]
Plot[Abs[p[4,1,0,r,0,0]]^2,{r,0,35},PlotPoints->60,PlotRange->{0,.0008}]
----
pr[n_,l_,r_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]]) (2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n] 

Integrate[pr[1,0,r2]^2 r2 Integrate[pr[2,1,r1]^2 r1^2, {r1,0,r2}],{r2,0,Infinity}]
=13/729

Integrate[pr[1,0,r2]^2 r2^2 Integrate[pr[2,1,r1]^2 r1, {r1,r2,Infinity}],{r2,0,Infinity}]
=164/729

Total= 59/243 = 0.242798

Other way around:

Integrate[pr[2,1,r2]^2 r2 Integrate[pr[1,0,r1]^2 r1^2, {r1,0,r2}],{r2,0,Infinity}]
Integrate[pr[2,1,r2]^2 r2^2 Integrate[pr[1,0,r1]^2 r1, {r1,r2,Infinity}],{r2,0,Infinity}]


--

Integrate[pr[1,0,r2]^2 r2 Integrate[pr[2,0,r1]^2 r1^2, {r1,0,r2}],{r2,0,Infinity}]
=7/243

Integrate[pr[1,0,r2]^2 r2^2 Integrate[pr[2,0,r1]^2 r1, {r1,r2,Infinity}],{r2,0,Infinity}]
=44/243

Total=17/81=0.209877

Other way around:

Integrate[pr[2,0,r2]^2 r2 Integrate[pr[1,0,r1]^2 r1^2, {r1,0,r2}],{r2,0,Infinity}]
Integrate[pr[2,0,r2]^2 r2^2 Integrate[pr[1,0,r1]^2 r1, {r1,r2,Infinity}],{r2,0,Infinity}]
--

Integrate[pr[1,0,r2]^2 r2 Integrate[pr[1,0,r1]^2 r1^2, {r1,0,r2}],{r2,0,Infinity}]
=5/16

Integrate[pr[1,0,r2]^2 r2^2 Integrate[pr[1,0,r1]^2 r1, {r1,r2,Infinity}],{r2,0,Infinity}]
=5/16

Total=5/8=0.625 checked

--Exchange:
Integrate[pr[1,0,r2]pr[2,0,r2]  r2 Integrate[pr[2,0,r1]pr[1,0,r1] r1^2, {r1,0,r2}],{r2,0,Infinity}]
=8/729

Integrate[pr[1,0,r2]pr[2,0,r2]  r2^2 Integrate[pr[2,0,r1]pr[1,0,r1]  r1, {r1,r2,Infinity}],{r2,0,Infinity}]
= must be same as above but get error:
Integrate[pr[2,0,r1]pr[1,0,r1]  r1, {r1,r2,Infinity}]=
4 Sqrt[2]/27 - Sqrt[2] r2^2 Exp [-3 r2/2]/3 - 4 Sqrt[2] Gamma[2, 0, 3 r2/2]/27

All except last term: Integrate[pr[1,0,r2]pr[2,0,r2]  r2^2 (4 Sqrt[2]/27 - Sqrt[2] r2^2 Exp [-3 r2/2]/3),{r2,0,Infinity}]
yields: -8/729

Gamma[2, 0, 3 r2/2] = 1-(1+3 r2/2)Exp[-3 r2/2]
Integrate[pr[1,0,r2]pr[2,0,r2]  r2^2 4 Sqrt[2] (1-(1+3 r2/2)Exp[-3 r2/2])/27,{r2,0,Infinity}]
-16/729 (so with "-", is 16/729

SubTotal=-8/729+16/729=8/729

Total= 16/729= 0.0219479

Integrate[pr[1,0,r2]pr[2,1,r2]   Integrate[pr[2,1,r1]pr[1,0,r1] r1^3, {r1,0,r2}],{r2,0,Infinity}]
=56/2187

Integrate[pr[1,0,r2]pr[2,1,r2]  r2^3 Integrate[pr[2,1,r1]pr[1,0,r1]  , {r1,r2,Infinity}],{r2,0,Infinity}]
=56/2187

Total= 112/2187 = 0.0512117

-(1+1/4)/2+(0.209877-0.0219479)/2
-(1+1/4)/2+(0.242798+0.0170706)/2
-----

Integrate[pr[2,1,r]pr[2,0,r]r^3,{r,0,Infinity}]

----
p[n_,l_,m_,r_,theta_,phi_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]])(2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n]  SphericalHarmonicY[l,m,theta,phi] 
sp[x_,z_]=(p[2,0,0,Sqrt[x^2+z^2],ArcTan[z,x],0]+p[2,1,0,Sqrt[x^2+z^2],ArcTan[z,x],0])/Sqrt[2]
sp2[x_,z_]=(2 - Sqrt[x^2  + z^2 ] + z)/(8 Sqrt[Pi] Exp[Sqrt[x^2  + z^2 ]/2])
sp1[x_,z_]=(2 - Sqrt[x^2  + z^2 ] - z)/(8 Sqrt[Pi] Exp[Sqrt[x^2  + z^2 ]/2])

ContourPlot[Abs[sp2[x,z]]^2,{x,-10,10},{z,-10,10},PlotPoints->100,
ContourLines->False,PlotRange->{0,.01}]
ContourPlot[Abs[sp1[x,z]]^2,{x,-10,10},{z,-10,10},PlotPoints->100,
ContourLines->False,PlotRange->{0,.01}]
Plot[Abs[sp2[0,r]]^2,{r,-10,10},PlotPoints->60,PlotRange->{0,.02}]
Plot[Abs[sp1[0,r]]^2,{r,-10,10},PlotPoints->60,PlotRange->{0,.02}]
---
Plot[{-30,-15,0,15,30},{x,-.15,.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-45/2,-15/2,15/2,45/2},{x,.85,1.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-15,0,15},{x,1.85,2.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-15/2,15/2},{x,2.85,3.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{0},{x,3.85,4.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-45/2,-15/2,15/2,45/2},{x,-.85,-1.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-15,0,15},{x,-1.85,-2.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{-15/2,15/2},{x,-2.85,-3.15},PlotStyle->{RGBColor[1,0,0]}]
Plot[{0},{x,-3.85,-4.15},PlotStyle->{RGBColor[1,0,0]}]
Show[%%%%%%%%%,%%%%%%%%,%%%%%%%,%%%%%%,%%%%%,%%%%,%%%,%%,%]
---

ephi[x_,y_,z_]=1/Sqrt[(x-1)^2+y^2+z^2]+1/Sqrt[(x+1)^2+y^2+z^2]+
1/Sqrt[(y-1)^2+x^2+z^2]+1/Sqrt[(y+1)^2+x^2+z^2]+
1/Sqrt[(z-1)^2+y^2+x^2]+1/Sqrt[(z+1)^2+y^2+x^2]

Series[ephi[x,y,z],{x,0,6},{y,0,6},{z,0,6}]
o4xyz[x_,y_,z_]=7(z^4+y^4+x^4)/2 -(21/2)(z^2 y^2 +y^2 x^2+z^2x^2)
o6xyz[x_,y_,z_]=3(z^6+y^6+x^6)/4-45(z^4 y^2+ z^2y^4+z^4x^2+y^4x^2+z^2x^4+y^2x^4)/8 + 135(z^2y^2x^2)/2

o42[theta_,phi_]=o4xyz[Sin[theta] Cos[phi],Sin[theta] Sin[phi],Cos[theta]]
o62[theta_,phi_]=o6xyz[Sin[theta] Cos[phi],Sin[theta] Sin[phi],Cos[theta]]

Series[ephi[r Sin[theta] Cos[phi],r Sin[theta] Sin[phi], r Cos[theta]],{r,0,6}] 
Simplify[%]
o4[theta_,phi_]=(7 (18 + 30 Cos[4 phi] - 20 Cos[4 phi - 2 theta] +
5 Cos[4 (phi - theta)] + 40 Cos[2 theta] + 70 Cos[4 theta] +
5 Cos[4 (phi + theta)] - 20 Cos[4 phi + 2 theta]))/256

Integrate[Sin[theta]Integrate[o42[theta,phi] SphericalHarmonicY[4,4,theta,phi] ,{phi,0,2 Pi}],{theta,0,Pi}]
4,4=35 Sqrt[Pi/70]/3
4,3=0
4,2=0
4,1=0
4,0=7 Sqrt[Pi]/3
4,-1=0
4,-2=0
4,-3=0
4,-4=35 Sqrt[Pi/70]/3

Simplify[(35 Sqrt[Pi/70]/3)(SphericalHarmonicY[4,4,theta,phi]+SphericalHarmonicY[4,-4,theta,phi])+
(7 Sqrt[Pi]/3) SphericalHarmonicY[4,0,theta,phi]-o4[theta,phi]]

o6[theta_,phi_]=(3 (100 - 420 Cos[4 phi] - 231 Cos[4 phi - 6 theta] -
105 Cos[4 phi - 2 theta] + 546 Cos[4 (phi - theta)] + 
210 Cos[2 theta] + 252 Cos[4 theta] + 462 Cos[6 theta] +
546 Cos[4 (phi + theta)] - 105 Cos[4 phi + 2 theta] -
231 Cos[4 phi + 6 theta]))/4096
Integrate[Sin[theta]Integrate[o6[theta,phi] SphericalHarmonicY[6,6,theta,phi] ,{phi,0,2 Pi}],{theta,0,Pi}]
6,6=0
6,5=0
6,4=-21 Sqrt[31 Pi/14]/26
6,3=0
6,2=0
6,1=0
6,0=3 Sqrt[13 Pi]/26

Simplify[(-21 Sqrt[13 Pi/14]/26)(SphericalHarmonicY[6,4,theta,phi]+SphericalHarmonicY[6,-4,theta,phi])+
(3 Sqrt[13 Pi]/26) SphericalHarmonicY[6,0,theta,phi]-o6[theta,phi]]


Series[ephi[Sin[theta] Cos[phi]/r,Sin[theta] Sin[phi]/r, Cos[theta]/r],{r,0,8}] 
Simplify[%]

phismall[x_,y_,z_]=6+7(z^4+y^4+x^4)/2 -(21/2)(z^2 y^2 +y^2 x^2+z^2x^2)+
3(z^6+y^6+x^6)/4-45(z^4 y^2+ z^2y^4+z^4x^2+y^4x^2+z^2x^4+y^2x^4)/8 + 135(z^2y^2x^2)/2

phibig[x_,y_,z_]=(6+(7(z^4+y^4+x^4)/2 -(21/2)(z^2 y^2 +y^2 x^2+z^2x^2))/(x^2+y^2+z^2)^(4)+
(3(z^6+y^6+x^6)/4-45(z^4 y^2+ z^2y^4+z^4x^2+y^4x^2+z^2x^4+y^2x^4)/8 + 135(z^2y^2x^2)/2)/(x^2+y^2+z^2)^(6))/Sqrt[x^2+y^2+z^2]

phiapprox[x_,y_,z_]:=If[x^2+y^2+z^2<1,phismall[x,y,z],phibig[x,y,z]]

Plot[{phiapprox[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]],ephi[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]]},{x,0,2},PlotRange->{0,6}]
Plot[{phiapprox[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]],ephi[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]]},{x,0,2},PlotRange->{0,6},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,0]}}]
Plot[{phiapprox[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]],ephi[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]]},{x,.8,1.2},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,0]}}]
Plot[ephi[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]],{x,0,2},PlotRange->{0,6}]

Plot[{phiapprox[x/Sqrt[2],x/Sqrt[2],0],ephi[x/Sqrt[2],x/Sqrt[2],0]},{x,0,2},PlotRange->{0,6},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,0]}}]

Plot[{phiapprox[x,0,0],ephi[x,0,0]},{x,0,2},PlotRange->{0,12},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,0]}}]

ContourPlot[ephi[x,y,0],{x,-2,2},{y,-2,2},PlotPoints->100,
Contours->25,PlotRange->{0,9},ContourShading->False]
ContourPlot[ephi[x/Sqrt[2],x/Sqrt[2],z],{x,-2,2},{z,-2,2},PlotPoints->100,
Contours->25,PlotRange->{0,9},ContourShading->False]
ContourPlot[ephi[x/Sqrt[2],x/Sqrt[2],z],{x,-2,2},{z,-2,2},PlotPoints->100,
Contours->{9,8.5,8,7.5,7,6.5,6.05,6,5.95,5.5,5,4.5,4,3.5,3,2.5,2},PlotRange->{0,9},ContourShading->False]
ContourPlot[ephi[x,y,0],{x,-2,2},{y,-2,2},PlotPoints->100,
Contours->{9,8.5,8,7.5,7,6.5,6.05,6,5.95,5.5,5,4.5,4,3.5,3,2.5,2},PlotRange->{0,9},ContourShading->False]

ContourPlot[phiapprox[x,y,0],{x,-2,2},{y,-2,2},PlotPoints->100,
Contours->{9,8.5,8,7.5,7,6.5,6.05,6,5.95,5.5,5,4.5,4,3.5,3,2.5,2},PlotRange->{0,9},ContourShading->False]
ContourPlot[phiapprox[x/Sqrt[2],x/Sqrt[2],z],{x,-2,2},{z,-2,2},PlotPoints->100,
Contours->{9,8.5,8,7.5,7,6.5,6.05,6,5.95,5.5,5,4.5,4,3.5,3,2.5,2},PlotRange->{0,9},ContourShading->False]
----

Integrate[Sin[theta]Integrate[SphericalHarmonicY[2,2,theta,phi]SphericalHarmonicY[2,2,theta,phi] \
(SphericalHarmonicY[4,-4,theta,phi]+SphericalHarmonicY[4,4,theta,phi]),{phi,0,2 Pi}],{theta,0,Pi}]

int2[l1_,m1_,l2_,m2_,l3_,m3_]:=(-1)^(-m1)Sqrt[(2l1+1)(2l2+1)(2l3+1)/(4 Pi)]\
ThreeJSymbol[{l1,0},{l2,0},{l3,0}]ThreeJSymbol[{l1,-m1},{l2,m2},{l3,m3}]
int2[2,-2,4,4,2,2]+int2[2,-2,4,-4,2,2]
Sqrt[5/(14 Pi)]

int2[2,2,4,0,2,2]=1/(14 Sqrt[Pi])
int2[2,1,4,0,2,1]=-2/(7 Sqrt[Pi])
int2[2,0,4,0,2,0] = 3/(7 Sqrt[Pi])
int2[2,-1,4,0,2,-1] = -2/(7 Sqrt[Pi])

int2[2,-2,4,0,2,-2] = 1/(14 Sqrt[Pi])

m40 = {{1/(14 Sqrt[Pi]),0,0,0,0},\
{0,-2/(7 Sqrt[Pi]),0,0,0},\
{0,0,3/(7 Sqrt[Pi]),0,0},\
{0,0,0,-2/(7 Sqrt[Pi]),0},\
{0,0,0,0,1/(14 Sqrt[Pi])}}

m44={{0,0,0,0,Sqrt[5/(14 Pi)]},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{Sqrt[5/(14 Pi)],0,0,0,0}}

m=(7Sqrt[Pi]/3)m40+Sqrt[35 Pi/18]m44

Eigensystem[m]

----
p[n_,l_,m_,r_,theta_,phi_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]])(2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n]  SphericalHarmonicY[l,m,theta,phi] 

dxy[x_,y_,z_]=(p[3,2,2,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]]-
p[3,2,-2,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]])/(Sqrt[2] I)
dx2y2[x_,y_,z_]=(p[3,2,2,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]]+
p[3,2,-2,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]])/(Sqrt[2])
dxz[x_,y_,z_]=(p[3,2,1,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]]+
p[3,2,-1,Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]])/(Sqrt[2])

N[dxy[4,4,0]]
Plot[Abs[dxy[x/Sqrt[2],x/Sqrt[2],0]]^2,{x,0,15}]

ContourPlot[Abs[dxy[x/Sqrt[2],x/Sqrt[2],z]]^2,{x,-16,16},{z,-16,16},PlotPoints->100,
ContourLines->False,PlotRange->{0,.0006}]
ContourPlot[Abs[dxy[x,y,0]]^2,{x,-16,16},{y,-16,16},PlotPoints->100,
ContourLines->False,PlotRange->{0,.0006}]
ContourPlot[Abs[dx2y2[x,y,0]]^2,{x,-16,16},{y,-16,16},PlotPoints->100,
ContourLines->False,PlotRange->{0,.0006}]

ContourPlot[Abs[dxz[x,0,z]]^2,{x,-16,16},{z,-16,16},PlotPoints->100,
ContourLines->False,PlotRange->{0,.0006}]
-----
p[n_,l_,m_,r_,theta_,phi_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]])(2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n]  SphericalHarmonicY[l,m,theta,phi] 
px[r_,theta_,phi_]=(-p[2,1,1,r,theta,phi]+p[2,1,-1,r,theta,phi])/Sqrt[2]
py[r_,theta_,phi_]=(p[2,1,1,r,theta,phi]+p[2,1,-1,r,theta,phi]) I/Sqrt[2]
pz[r_,theta_,phi_]=p[2,1,0,r,theta,phi]
s[r_,theta_,phi_]=p[2,0,0,r,theta,phi]
teta[r_,theta_,phi_]=(-s[r,theta,phi]+px[r,theta,phi]+py[r,theta,phi]+pz[r,theta,phi])/2
teta2[x_,y_,z_]=Simplify[teta[Sqrt[x^2+y^2+z^2],ArcTan[z,Sqrt[x^2+y^2]],ArcTan[x,y]]]

Plot[teta[r,.9553,Pi/4],{r,0,5}]
Plot[teta[r,Pi-.9553,Pi/4+Pi],{r,0,5}]

ContourPlot[Abs[teta2[x/Sqrt[2],x/Sqrt[2],z]]^2,{x,-8,8},{z,-8,8},PlotPoints->100,
ContourLines->False,PlotRange->{0,.006}]

Plot[teta2[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]],{x,-10,10}]
Plot[Abs[teta2[x/Sqrt[3],x/Sqrt[3],x/Sqrt[3]]]^2,{x,-10,10}]
----

LaguerreL[5,9,x]

--
Integrate[Exp[-2 a r]r^2,{r,0,Infinity}]=1/(4a^3)
Integrate[Exp[-2 a r](-a^2+(-2+2a)/r)r^2,{r,0,Infinity}]=(-2+a)/(4a^2)

Integrate[Exp[-2 a r2]r2 Integrate[Exp[-2 a r]r^2,{r,0,r2}],{r2,0,Infinity}]+
Integrate[Exp[-2 a r2]r2^2 Integrate[Exp[-2 a r]r,{r,r2,Infinity}],{r2,0,Infinity}]=5/(128a^5)

a(a-2)+5a/(z 8) /. a->(1-5/(16 z))
----
ephi[x_,y_,z_]=1/Sqrt[(x-1)^2+(y-1)^2+(z-1)^2]+1/Sqrt[(x+1)^2+(y+1)^2+(z-1)^2]+
1/Sqrt[(y-1)^2+(x+1)^2+(z+1)^2]+1/Sqrt[(y+1)^2+(x-1)^2+(z+1)^2]

Series[ephi[r Sin[theta] Cos[phi],r Sin[theta] Sin[phi], r Cos[theta]],{r,0,6}] 
Simplify[%]

SeriesData[r, 0, {4/3^(1/2), 0, 0, 
    (10*Cos[theta]*Sin[2*phi]*Sin[theta]^2)/(9*3^(1/2)), 
    (7*(-18 - 30*Cos[4*phi] + 20*Cos[4*phi - 2*theta] - 
         5*Cos[4*(phi - theta)] - 40*Cos[2*theta] - 70*Cos[4*theta] - 
         5*Cos[4*(phi + theta)] + 20*Cos[4*phi + 2*theta]))/(5184*3^(1/2)), 
    0, (100 - 420*Cos[4*phi] - 231*Cos[4*phi - 6*theta] - 
       105*Cos[4*phi - 2*theta] + 546*Cos[4*(phi - theta)] + 
       210*Cos[2*theta] + 252*Cos[4*theta] + 462*Cos[6*theta] + 
       546*Cos[4*(phi + theta)] - 105*Cos[4*phi + 2*theta] - 
       231*Cos[4*phi + 6*theta])/(31104*3^(1/2))}, 0, 7, 1]

Note: O3 term has no effect

o4[theta_,phi_]=(7 (18 + 30 Cos[4 phi] - 20 Cos[4 phi - 2 theta] +
5 Cos[4 (phi - theta)] + 40 Cos[2 theta] + 70 Cos[4 theta] +
5 Cos[4 (phi + theta)] - 20 Cos[4 phi + 2 theta]))/(5184*3^(1/2))

Integrate[Sin[theta]Integrate[o4[theta,phi] SphericalHarmonicY[4,4,theta,phi] ,{phi,0,2 Pi}],{theta,0,Pi}]
4,4=(4*((35*Pi)/6)^(1/2))/243
4,3=0
4,2=0
4,1=0
4,0=(28*(Pi/3)^(1/2))/243
4,-1=0
4,-2=0
4,-3=0
4,-4=(4*((35*Pi)/6)^(1/2))/243

