Plot[{Pi^2,4 Pi^2,9 Pi^2,16 Pi^2,25 Pi^2,36 Pi^2},{x,0,1},
PlotRange->{0,400}]

---
u0=120
kg[1]=2.5
kg[2]=5
kg[3]=8
kg[4]=10

feven[k_]=k Cot[k/2]+Sqrt[u0-k^2]
fodd[k_]=k Tan[k/2]-Sqrt[u0-k^2]

k[4]=x /. N[FindRoot[feven[x],{x,kg[4]},AccuracyGoal->13],12]
k[3]=x /. N[FindRoot[fodd[x],{x,kg[3]},AccuracyGoal->13],12]
k[2]=x /. N[FindRoot[feven[x],{x,kg[2]},AccuracyGoal->13],12]
k[1]=x /. N[FindRoot[fodd[x],{x,kg[1]},AccuracyGoal->13],12]

For[i=1,i<5,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Cos[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Cos[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

peven[i_,x_]=If[Abs[x]<.5, a[i]Cos[k[i]x], a[i]ba[i]Exp[-kappa[i]Abs[x]]]

For[i=2,i<5,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Sin[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Sin[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

podd[i_,x_]=If[Abs[x]<.5, a[i]Sin[k[i]x], Sign[x]a[i]ba[i]Exp[-kappa[i]Abs[x]]]

p[i_,x_]:=If[Mod[i,2]==0, podd[i,x], peven[i,x]]
u[x_]:=If[Abs[x]<.5,0,u0]
p2[i_,x_]:=If[Mod[i,2]==0, 5 podd[i,x],5 peven[i,x]]

Plot[{(p2[1,x])+k[1]^2,(p2[2,x])+k[2]^2,(p2[3,x])+k[3]^2,(p2[4,x])+k[4]^2,
u[x],k[1]^2,k[2]^2,k[3]^2,k[4]^2},{x,-1,1}, PlotPoints->100, \
PlotRange->{0,130},PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},
 {RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,1]},     
 {RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]}}]
----
Plot[p[4,x],{x,-1,1}]
Plot[{p[4,x],k[4]a[4]Cos[k[4]/2](x-.5)+a[4]Sin[k[4]/2]},{x,.45,.6},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,1]}},Axes->False]
----
u0=100
kg[1]=2.5
kg[2]=5
kg[3]=8
kg[4]=9.5

feven[k_]=k Cot[k/2]+Sqrt[u0-k^2]
fodd[k_]=k Tan[k/2]-Sqrt[u0-k^2]

k[4]=x /. N[FindRoot[feven[x],{x,kg[4]},AccuracyGoal->13],12]
k[3]=x /. N[FindRoot[fodd[x],{x,kg[3]},AccuracyGoal->13],12]
k[2]=x /. N[FindRoot[feven[x],{x,kg[2]},AccuracyGoal->13],12]
k[1]=x /. N[FindRoot[fodd[x],{x,kg[1]},AccuracyGoal->13],12]

For[i=1,i<5,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Cos[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Cos[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

peven[i_,x_]=If[Abs[x]<.5, a[i]Cos[k[i]x], a[i]ba[i]Exp[-kappa[i]Abs[x]]]

For[i=2,i<5,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Sin[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Sin[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

podd[i_,x_]=If[Abs[x]<.5, a[i]Sin[k[i]x], Sign[x]a[i]ba[i]Exp[-kappa[i]Abs[x]]]

p[i_,x_]:=If[Mod[i,2]==0, podd[i,x], peven[i,x]]
u[x_]:=If[Abs[x]<.5,0,u0]
p2[i_,x_]:=If[Mod[i,2]==0, 5 podd[i,x],5 peven[i,x]]

Plot[{u[x],(p2[1,x])+k[1]^2,(p2[2,x])+k[2]^2,(p2[3,x])+k[3]^2,(p2[4,x])+k[4]^2,
k[1]^2,k[2]^2,k[3]^2,k[4]^2},{x,-1,1}, PlotPoints->100, \
PlotRange->{0,130},PlotStyle->{{RGBColor[1,0,0]},{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]}}]

--
u0=80
kg[1]=2.5
kg[2]=5
kg[3]=8

feven[k_]=k Cot[k/2]+Sqrt[u0-k^2]
fodd[k_]=k Tan[k/2]-Sqrt[u0-k^2]

k[3]=x /. N[FindRoot[fodd[x],{x,kg[3]},AccuracyGoal->13],12]
k[2]=x /. N[FindRoot[feven[x],{x,kg[2]},AccuracyGoal->13],12]
k[1]=x /. N[FindRoot[fodd[x],{x,kg[1]},AccuracyGoal->13],12]

For[i=1,i<4,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Cos[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Cos[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

peven[i_,x_]=If[Abs[x]<.5, a[i]Cos[k[i]x], a[i]ba[i]Exp[-kappa[i]Abs[x]]]

For[i=2,i<4,i=i+2,kappa[i]=Sqrt[u0-k[i]^2];
ba[i]=Sin[k[i]/2]/(Exp[-kappa[i]/2]);
a[i]=.5/(Integrate[Sin[k[i]x]^2,{x,0,.5}]+ba[i]^2 Integrate[Exp[-2 kappa[i]x],{x,.5,Infinity}])]

podd[i_,x_]=If[Abs[x]<.5, a[i]Sin[k[i]x], Sign[x]a[i]ba[i]Exp[-kappa[i]Abs[x]]]

p[i_,x_]:=If[Mod[i,2]==0, podd[i,x], peven[i,x]]
u[x_]:=If[Abs[x]<.5,0,u0]
p2[i_,x_]:=If[Mod[i,2]==0, 5 podd[i,x],5 peven[i,x]]

Plot[{(p2[1,x])+k[1]^2,(p2[2,x])+k[2]^2,(p2[3,x])+k[3]^2,
u[x],k[1]^2,k[2]^2,k[3]^2},{x,-1,1}, PlotPoints->100, \
PlotRange->{0,130},PlotStyle->{{RGBColor[0,0,0]},
 {RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[1,0,0]},     
 {RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]}}]

aint[i]=aint[i-1]+NIntegrate[AiryAi[x]^2,{x,a[i],a[i-1]}];
n[i]=1/Sqrt[aint[i]]]
----
uieven[k_]=(Cos[k/2]/k)^2
uiodd[k_]=(Sin[k/2]/k)^2
ParametricPlot[{uieven[k],k/Pi},{k,1 Pi,0 Pi}]
ParametricPlot[{uieven[k],k/Pi},{k,3 Pi,2 Pi}]
ParametricPlot[{uieven[k],k/Pi},{k,5 Pi,4 Pi}]
ParametricPlot[{uieven[k],k/Pi},{k,7 Pi,6 Pi}]
ParametricPlot[{uiodd[k],k/Pi},{k,2 Pi,1 Pi}]
ParametricPlot[{uiodd[k],k/Pi},{k,4 Pi,3 Pi}]
ParametricPlot[{uiodd[k],k/Pi},{k,6 Pi,5 Pi}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,PlotRange->{{0,.03},{0,7.5}}]
Show[%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,PlotRange->{{0,.03},{0,7.5}}]
----
x/: Conjugate[x]=x
k/: Conjugate[k]=k
q/: Conjugate[q]=q

psiL[x_]=Exp[I k x]+r Exp[-I k x]
psi0[x_]=a Cos[q x]+b Sin[q x]
psiR[x_]=t Exp[I k x]

jL=(hbar/(2 m I))(Conjugate[psiL[x]]D[psiL[x],x]-Conjugate[D[psiL[x],x]] psiL[x])
Simplify[%]
hbar k (1 - r Conjugate[r])/m

j0=(hbar/(2 m I))(Conjugate[psi0[x]]D[psi0[x],x]-Conjugate[D[psi0[x],x]] psi0[x])
Simplify[%]
%  /. Conjugate[c_ + d_]->Conjugate[c] + Conjugate[d]
(-I/2*hbar*q*(b*Conjugate[a] - a*Conjugate[b]))/m

jR=(hbar/(2 m I))(Conjugate[psiR[x]]D[psiR[x],x]-Conjugate[D[psiR[x],x]] psiR[x])
hbar k t Conjugate[t]/m
----
r[k_,l_]=If[l==0,0,I(1-k^2)Sin[l]/(2k Cos[l]-I(k^2+1)Sin[l])]
ContourPlot[Abs[r[Sqrt[k],l]]^2,{l,0,3 Pi},{k,0,1},PlotPoints->100,ContourLines->False]
ContourPlot[Abs[r[e/(e+u),Sqrt[e+u]]]^2,{u,0,100},{e,0,100},PlotPoints->100,ContourLines->False]
ParametricPlot[{Sqrt[e+9],e/(e+9)},{e,0 ,(3 Pi)^2-9},PlotStyle->{RGBColor[1,0,0]},
PlotRange->{{0,3 Pi},{0,1}}]
Show[%%,%]
----
e=2 Pi^2
u=7 Pi^2
r=I(u)Sin[Sqrt[e+u]]/(2 Sqrt[e(e+u)] Cos[Sqrt[e+u]]-I(2e+u)Sin[Sqrt[e+u]])
a=r+1
b=I(1-r)Sqrt[e/(e+u)]
t=Exp[-I Sqrt[e]](a Cos[Sqrt[e+u]] + b Sin[Sqrt[e+u]])
p[x_]:=If[x<0,Exp[I Sqrt[e] x]+r Exp[-I Sqrt[e] x],If[x>1,t Exp[I Sqrt[e] x], a Cos[Sqrt[e+u] x]+b Sin[Sqrt[e+u] x]]]
Plot[Re[p[x]],{x,-3,4}]
Plot[Abs[p[x]]^2,{x,-3,4},PlotRange->{0,1.2}]
Plot[Im[p[x]],{x,-3,4}]
--
e=.5 Pi^2
u=7 Pi^2
r=I(u)Sin[Sqrt[e+u]]/(2 Sqrt[e(e+u)] Cos[Sqrt[e+u]]-I(2e+u)Sin[Sqrt[e+u]])
a=r+1
b=I(1-r)Sqrt[e/(e+u)]
t=Exp[-I Sqrt[e]](a Cos[Sqrt[e+u]] + b Sin[Sqrt[e+u]])
p[x_]:=If[x<0,Exp[I Sqrt[e] x]+r Exp[-I Sqrt[e] x],If[x>1,t Exp[I Sqrt[e] x], a Cos[Sqrt[e+u] x]+b Sin[Sqrt[e+u] x]]]
Plot[Re[p[x]],{x,-3,4}]
Plot[Re[p[x]],{x,-.02,.02}]
Plot[Abs[p[x]]^2,{x,-3,4},PlotRange->{0,3.5}]
Plot[Im[p[x]],{x,-3,4}]

--
e=2 Pi^2+1
u=7 Pi^2
r=I(u)Sin[Sqrt[e+u]]/(2 Sqrt[e(e+u)] Cos[Sqrt[e+u]]-I(2e+u)Sin[Sqrt[e+u]])
a=r+1
b=I(1-r)Sqrt[e/(e+u)]
t=Exp[-I Sqrt[e]](a Cos[Sqrt[e+u]] + b Sin[Sqrt[e+u]])
p[x_]:=If[x<0,Exp[I Sqrt[e] x]+r Exp[-I Sqrt[e] x],If[x>1,t Exp[I Sqrt[e] x], a Cos[Sqrt[e+u] x]+b Sin[Sqrt[e+u] x]]]
Plot[Re[p[x]],{x,-3,4}]
----
p[x_,y_]=Sin[.5 Pi x]^2 Sin[4 Pi y/3]^2
ContourPlot[p[x,y],{x,0,4},{y,0,3},PlotPoints->100,ContourLines->False,AspectRatio->.75]
---
Plot[{BesselJ[0,x],BesselJ[1,x],BesselJ[2,x],BesselJ[3,x]},{x,0,10},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{BesselY[0,x],BesselY[1,x],BesselY[2,x],BesselY[3,x]},{x,0,10},
PlotRange->{-2,1},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{BesselI[0,x],BesselI[1,x],BesselI[2,x],BesselI[3,x]},{x,0,5},
PlotRange->{0,10},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{BesselK[0,x],BesselK[1,x],BesselK[2,x],BesselK[3,x]},{x,0,3},
PlotRange->{0,10},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

<<Graphics`Master`

LogPlot[{BesselK[0,x],BesselK[1,x],BesselK[2,x],BesselK[3,x]},{x,0,5},
PlotRange->{.01,100},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

----
<<NumericalMath`BesselZeros`
Plot[Evaluate[BesselJZeros[0,4]^2],{x,-.15,.15}]
Plot[Evaluate[BesselJZeros[1,3]^2],{x,.85,1.15}]
Plot[Evaluate[BesselJZeros[1,3]^2],{x,-1.15,-.85}]
Plot[Evaluate[BesselJZeros[2,3]^2],{x,1.85,2.15}]
Plot[Evaluate[BesselJZeros[2,3]^2],{x,-2.15,-1.85}]
Plot[Evaluate[BesselJZeros[3,2]^2],{x,2.85,3.15}]
Plot[Evaluate[BesselJZeros[3,2]^2],{x,-3.15,-2.85}]
Plot[Evaluate[BesselJZeros[4,2]^2],{x,3.85,4.15}]
Plot[Evaluate[BesselJZeros[4,2]^2],{x,-4.15,-3.85}]
Plot[Evaluate[BesselJZeros[5,1]^2],{x,4.85,5.15}]
Plot[Evaluate[BesselJZeros[5,1]^2],{x,-5.15,-4.85}]
Plot[Evaluate[BesselJZeros[6,1]^2],{x,5.85,6.15}]
Plot[Evaluate[BesselJZeros[6,1]^2],{x,-6.15,-5.85}]
Plot[Evaluate[BesselJZeros[7,1]^2],{x,6.85,7.15}]
Plot[Evaluate[BesselJZeros[7,1]^2],{x,-7.15,-6.85}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,%%%%%%%%%%%%,%%%%%%%%%%%%%,%%%%%%%%%%%%%%,%%%%%%%%%%%%%%%,PlotRange->{0,150}]
Show[%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,%%%%%%%%%%%%,%%%%%%%%%%%%%,%%%%%%%%%%%%%%,%%%%%%%%%%%%%%%,
%%%%%%%%%%%%%%%%,PlotRange->{0,100}]
Show[%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,%%%%%%%%%%%%,%%%%%%%%%%%%%,%%%%%%%%%%%%%%,%%%%%%%%%%%%%%%,
%%%%%%%%%%%%%%%%,%%%%%%%%%%%%%%%%%,PlotRange->{0,80}]
----
<<NumericalMath`BesselZeros`
roots={BesselJZeros[0,11],BesselJZeros[1,11],BesselJZeros[2,11],BesselJZeros[3,11],BesselJZeros[4,11],
BesselJZeros[5,11],BesselJZeros[6,11],BesselJZeros[7,11],BesselJZeros[8,11],BesselJZeros[9,11],
BesselJZeros[10,11]}
m /: IntegerQ[m]=True
nr /: IntegerQ[nr]=True

r[m_,nr_,r_]=BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2./(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])]
r2[m_,nr_,r_]=BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2./(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])]
m=0
nr=10
r010[r_]=BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2/(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])]
r2010[r_]=If[r<1,(BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2/(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])])^2,0]
ContourPlot[r2010[Sqrt[x^2+y^2]],{x,-1,1},{y,-1,1},PlotPoints->100,ContourLines->False]
m=10
nr=5
r105[r_]=BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2/(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])]
r2105[r_]=If[r<1,(BesselJ[m,roots[[m+1,nr+1]] r] Sqrt[2./(-BesselJ[m-1,roots[[m+1,nr+1]] ]BesselJ[m+1,roots[[m+1,nr+1]] ])])^2,0]
ContourPlot[r2105[Sqrt[x^2+y^2]],{x,-1,1},{y,-1,1},PlotPoints->100,ContourLines->False]

---
v[m_,r_]=If[r<1,m^2/r^2,201]
k00=roots[[1,1]]
k01=roots[[1,2]]
k02=roots[[1,3]]
k03=roots[[1,4]]
r00[r_]=BesselJ[0,k00 r] Sqrt[2./(-BesselJ[-1,k00]BesselJ[1,k00 ])]
r01[r_]=BesselJ[0,k01 r] Sqrt[2./(-BesselJ[-1,k01]BesselJ[1,k01 ])]
r02[r_]=BesselJ[0,k02 r] Sqrt[2./(-BesselJ[-1,k02]BesselJ[1,k02 ])]
r03[r_]=BesselJ[0,k03 r] Sqrt[2./(-BesselJ[-1,k03]BesselJ[1,k03 ])]
a=5
Plot[{a r00[r]+k00^2,a r01[r]+k01^2,a r02[r]+k02^2,a r03[r]+k03^2},{r,0,1},PlotRange->{0,200}]
Plot[{k00^2,k01^2,k02^2,k03^2,v[0,r]},{r,0,1.1},PlotRange->{0,200}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

k10=roots[[2,1]]
k11=roots[[2,2]]
k12=roots[[2,3]]
k13=roots[[2,4]]
r10[r_]=BesselJ[1,k10 r] Sqrt[2./(-BesselJ[0,k10]BesselJ[2,k10 ])]
r11[r_]=BesselJ[1,k11 r] Sqrt[2./(-BesselJ[0,k11]BesselJ[2,k11 ])]
r12[r_]=BesselJ[1,k12 r] Sqrt[2./(-BesselJ[0,k12]BesselJ[2,k12 ])]
r13[r_]=BesselJ[1,k13 r] Sqrt[2./(-BesselJ[0,k13]BesselJ[2,k13 ])]
a=5
Plot[{a r10[r]+k10^2,a r11[r]+k11^2,a r12[r]+k12^2,a r13[r]+k13^2},{r,0,1},PlotRange->{0,200}]
Plot[{k10^2,k11^2,k12^2,k13^2,v[1,r]},{r,0,1.1},PlotRange->{0,200}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

k20=roots[[3,1]]
k21=roots[[3,2]]
k22=roots[[3,3]]
r20[r_]=BesselJ[2,k20 r] Sqrt[2./(-BesselJ[1,k20]BesselJ[3,k20 ])]
r21[r_]=BesselJ[2,k21 r] Sqrt[2./(-BesselJ[1,k21]BesselJ[3,k21 ])]
r22[r_]=BesselJ[2,k22 r] Sqrt[2./(-BesselJ[1,k22]BesselJ[3,k22 ])]
a=5
Plot[{a r20[r]+k20^2,a r21[r]+k21^2,a r22[r]+k22^2},{r,0,1},PlotRange->{0,200}]
Plot[{k20^2,k21^2,k22^2,v[2,r]},{r,0,1.1},PlotRange->{0,200}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

Show[GraphicsArray[{%32,%52,%69}]]

---
<<NumericalMath`BesselZeros`

g0=BesselJZeros[0,2]-.1
g1=BesselJZeros[1,2]-.1
g2=BesselJZeros[2,2]-.1
g3=BesselJZeros[3,2]-.1
g4=BesselJZeros[4,2]-.1

u0=50

besselJPrime[m_,k_]=D[BesselJ[m,k],k]
besselKPrime[m_,k_]=D[BesselK[m,k],k]

f[m_,k_]=k besselJPrime[m,k]/BesselJ[m,k]-Sqrt[u0-k^2] besselKPrime[m, Sqrt[u0-k^2]]/BesselK[m,Sqrt[u0-k^2]]
Plot[f[0,x],{x,0,Sqrt[50]}]
Plot[f[0,x],{x,7,Sqrt[50]}]
Plot[f[5,x],{x,0,Sqrt[50]}]

k[0,0]=x /. N[FindRoot[f[0,x],{x,g0[[1]]},AccuracyGoal->13],12]
k[0,1]=x /. N[FindRoot[f[0,x],{x,g0[[2]]},AccuracyGoal->13],12]
k[0,2]=x /. N[FindRoot[f[0,x],{x,Sqrt[50]-.0001},AccuracyGoal->13],12]
k[1,0]=x /. N[FindRoot[f[1,x],{x,g1[[1]]},AccuracyGoal->13],12]
k[1,1]=x /. N[FindRoot[f[1,x],{x,g1[[2]]},AccuracyGoal->13],12]
k[2,0]=x /. N[FindRoot[f[2,x],{x,g2[[1]]},AccuracyGoal->13],12]
k[2,1]=x /. N[FindRoot[f[2,x],{x,Sqrt[49.999]},AccuracyGoal->13],12]
k[3,0]=x /. N[FindRoot[f[3,x],{x,g3[[1]]},AccuracyGoal->13],12]
k[4,0]=x /. N[FindRoot[f[4,x],{x,7},AccuracyGoal->13],12]

un[m_,nr_,r_]=If[r<1,BesselJ[m,k[m,nr] r], BesselJ[m,k[m,nr]] BesselK[m,Sqrt[u0-k[m,nr]^2] r]/BesselK[m,Sqrt[u0-k[m,nr]^2]]]
Plot[un[0,0,x],{x,0,1.4}]
Plot[un[2,1,x],{x,.9,1.1}]


norm[0,0]=1/Sqrt[Integrate[BesselJ[0,k[0,0] r]^2 r,{r,0,1}]+ 
(BesselJ[0,k[0,0]]/BesselK[0,Sqrt[u0-k[0,0]^2]])^2 Integrate[BesselK[0,Sqrt[u0-k[0,0]^2] r]^2 r,{r,1,Infinity}]]
norm[0,1]=1/Sqrt[Integrate[BesselJ[0,k[0,1] r]^2 r,{r,0,1}]+ 
(BesselJ[0,k[0,1]]/BesselK[0,Sqrt[u0-k[0,1]^2]])^2 Integrate[BesselK[0,Sqrt[u0-k[0,1]^2] r]^2 r,{r,1,Infinity}]]
norm[0,2]=1/Sqrt[Integrate[BesselJ[0,k[0,2] r]^2 r,{r,0,1}]+ 
(BesselJ[0,k[0,2]]/BesselK[0,Sqrt[u0-k[0,2]^2]])^2 Integrate[BesselK[0,Sqrt[u0-k[0,2]^2] r]^2 r,{r,1,Infinity}]]

norm[1,0]=1/Sqrt[Integrate[BesselJ[1,k[1,0] r]^2 r,{r,0,1}]+ 
(BesselJ[1,k[1,0]]/BesselK[1,Sqrt[u0-k[1,0]^2]])^2 Integrate[BesselK[1,Sqrt[u0-k[1,0]^2] r]^2 r,{r,1,Infinity}]]
norm[1,1]=1/Sqrt[Integrate[BesselJ[1,k[1,1] r]^2 r,{r,0,1}]+ 
(BesselJ[1,k[1,1]]/BesselK[1,Sqrt[u0-k[1,1]^2]])^2 Integrate[BesselK[1,Sqrt[u0-k[1,1]^2] r]^2 r,{r,1,Infinity}]]

norm[2,0]=1/Sqrt[Integrate[BesselJ[2,k[2,0] r]^2 r,{r,0,1}]+ 
(BesselJ[2,k[2,0]]/BesselK[2,Sqrt[u0-k[2,0]^2]])^2 Integrate[BesselK[2,Sqrt[u0-k[2,0]^2] r]^2 r,{r,1,Infinity}]]
norm[2,1]=1/Sqrt[Integrate[BesselJ[2,k[2,1] r]^2 r,{r,0,1}]+ 
(BesselJ[2,k[2,1]]/BesselK[2,Sqrt[u0-k[2,1]^2]])^2 Integrate[BesselK[2,Sqrt[u0-k[2,1]^2] r]^2 r,{r,1,Infinity}]]

norm[3,0]=1/Sqrt[Integrate[BesselJ[3,k[3,0] r]^2 r,{r,0,1}]+ 
(BesselJ[3,k[3,0]]/BesselK[3,Sqrt[u0-k[3,0]^2]])^2 Integrate[BesselK[3,Sqrt[u0-k[3,0]^2] r]^2 r,{r,1,Infinity}]]

norm[4,0]=1/Sqrt[Integrate[BesselJ[4,k[4,0] r]^2 r,{r,0,1}]+ 
(BesselJ[4,k[4,0]]/BesselK[4,Sqrt[u0-k[4,0]^2]])^2 Integrate[BesselK[4,Sqrt[u0-k[4,0]^2] r]^2 r,{r,1,Infinity}]]

psi[m_,nr_,r_]=norm[m,nr] If[r<1,BesselJ[m,k[m,nr] r], BesselJ[m,k[m,nr]] BesselK[m,Sqrt[u0-k[m,nr]^2] r]/BesselK[m,Sqrt[u0-k[m,nr]^2]]]

NIntegrate[psi[3,0,r]^2 r,{r,0,1.5}]
v[m_,r_]=If[r<1,m^2/r^2,m^2/r^2+50]

a=5
Plot[{a psi[0,0,r]+k[0,0]^2,a psi[0,1,r]+k[0,1]^2,a psi[0,2,r]+k[0,2]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[0,0]^2,k[0,1]^2,k[0,2]^2,v[0,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50}}]

Plot[{a psi[1,0,r]+k[1,0]^2,a psi[1,1,r]+k[1,1]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[1,0]^2,k[1,1]^2,v[1,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50}}]

Plot[{a psi[2,0,r]+k[2,0]^2,a psi[2,1,r]+k[2,1]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[2,0]^2,k[2,1]^2,v[2,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50}}]

Show[GraphicsArray[{%102,%105,%99}]]

---
Plot[Evaluate[BesselJZeros[0,4]^2],{x,-.15,.15}]
Plot[Evaluate[BesselJZeros[1,3]^2],{x,.85,1.15}]
Plot[Evaluate[BesselJZeros[1,3]^2],{x,-1.15,-.85}]
Plot[Evaluate[BesselJZeros[2,3]^2],{x,1.85,2.15}]
Plot[Evaluate[BesselJZeros[2,3]^2],{x,-2.15,-1.85}]
Plot[Evaluate[BesselJZeros[3,2]^2],{x,2.85,3.15}]
Plot[Evaluate[BesselJZeros[3,2]^2],{x,-3.15,-2.85}]
Plot[Evaluate[BesselJZeros[4,2]^2],{x,3.85,4.15}]
Plot[Evaluate[BesselJZeros[4,2]^2],{x,-4.15,-3.85}]
Plot[Evaluate[BesselJZeros[5,1]^2],{x,4.85,5.15}]
Plot[Evaluate[BesselJZeros[5,1]^2],{x,-5.15,-4.85}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,PlotRange->{0,80}]
Plot[{k[0,0]^2,k[0,1]^2,k[0,2]^2},{x,-.15,.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[1,0]^2,k[1,1]^2},{x,.85,1.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[1,0]^2,k[1,1]^2},{x,-1.15,-.85},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[2,0]^2,k[2,1]^2},{x,1.85,2.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[2,0]^2,k[2,1]^2},{x,-2.15,-1.85},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[3,0]^2},{x,2.85,3.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[3,0]^2},{x,-3.15,-2.85},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[4,0]^2},{x,3.85,4.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[4,0]^2},{x,-4.15,-3.85},PlotStyle->{{RGBColor[1,0,0]}}]
Show[%%%%%%%%%%,%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
PlotRange->{0,80}]
---
besselJPrime[m_,k_]=D[BesselJ[m,k],k]
besselKPrime[m_,k_]=D[BesselK[m,k],k]

m=0
Limit[x besselKPrime[m, x]/BesselK[m,x],{x->0}]
m=1
FindRoot[k besselJPrime[m,k] + m BesselJ[m,k],{k,2}]
m=2
FindRoot[k besselJPrime[m,k] + m BesselJ[m,k],{k,5}]

f[m_,k_,u0_]=k besselJPrime[m,k]/BesselJ[m,k]-Sqrt[u0-k^2] besselKPrime[m, Sqrt[u0-k^2]]/BesselK[m,Sqrt[u0-k^2]]
f1[m_,k_,u0_]=k besselJPrime[m,k]/BesselJ[m,k]
f2[m_,k_,u0_]=Sqrt[u0-k^2] besselKPrime[m, Sqrt[u0-k^2]]/BesselK[m,Sqrt[u0-k^2]]

Plot[{f1[0,k,.2],f2[0,k,.2]},{k,0,Sqrt[.2]}]
Plot[{f1[1,k,50],f2[1,k,50]},{k,0,Sqrt[50]},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}},PlotRange->{-10,5}]
Plot[{f1[2,k,50],f2[2,k,50]},{k,0,Sqrt[50]},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}},PlotRange->{-10,10}]
Plot[{f1[1,k,5],f2[1,k,5]},{k,0,Sqrt[5]}]
Plot[{f1[0,k,.1],f2[0,k,.1]},{k,.3161,Sqrt[.1]}]
---
ls[w_]=(k[0,0]^2-w 2)^2+(k[0,1]^2-w 6)^2+(k[0,2]^2-w 10)^2+2((k[1,0]^2-w 4)^2+(k[1,1]^2-w 8)^2+
(k[2,0]^2-w 6)^2+(k[2,1]^2-w 10)^2+(k[3,0]^2-w 8)^2+(k[4,0]^2-w 10)^2)
Plot[ls[w],{w,1,10}]
D[ls[w],w]
FindRoot[%,{w, 4.2}]
w=4.31182
Plot[{w 2, w 6, w 10},{x,-.15,.15},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 4, w 8, w 12},{x,.85,1.15},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 4, w 8, w 12},{x,-1.15,-.85},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 6, w 10},{x,1.85,2.15},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 6, w 10},{x,-2.15,-1.85},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 8, w 12},{x,2.85,3.15},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[{w 8, w 12},{x,-3.15,-2.85},PlotStyle->{{RGBColor[0,1,0]},{RGBColor[0,1,0]}}]
Plot[w 10,{x,3.85,4.15},PlotStyle->{{RGBColor[0,1,0]}}]
Plot[w 10,{x,-4.15,-3.85},PlotStyle->{{RGBColor[0,1,0]}}]
Plot[w 12,{x,4.85,5.15},PlotStyle->{{RGBColor[0,1,0]}}]
Plot[w 12,{x,-5.15,-4.85},PlotStyle->{{RGBColor[0,1,0]}}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,PlotRange->{0,80}]
Plot[{k[0,0]^2,k[0,1]^2,k[0,2]^2},{x,-.15,.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[1,0]^2,k[1,1]^2},{x,.85,1.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[1,0]^2,k[1,1]^2},{x,-1.15,-.85},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[2,0]^2,k[2,1]^2},{x,1.85,2.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[2,0]^2,k[2,1]^2},{x,-2.15,-1.85},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[3,0]^2},{x,2.85,3.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[3,0]^2},{x,-3.15,-2.85},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[4,0]^2},{x,3.85,4.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[4,0]^2},{x,-4.15,-3.85},PlotStyle->{{RGBColor[1,0,0]}}]
Show[%%%%%%%%%%,%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
PlotRange->{0,80}]
-----
<<NumericalMath`BesselZeros`

g0=BesselJZeros[0,5]
nd0[k_]=If[k<g0[[1]],0,If[k<g0[[3]],1,If[k<g0[[5]],2,3]]]
g1=BesselJZeros[1,5]
nd1[k_]=If[k<g1[[1]],0,If[k<g1[[3]],1,If[k<g1[[5]],2,3]]]
g2=BesselJZeros[2,5]
nd2[k_]=If[k<g2[[1]],0,If[k<g2[[3]],1,If[k<g2[[5]],2,3]]]

d[m_,k_]=ArcTan[-BesselY[m,k],-BesselJ[m,k]]
d0[k_]=d[0,k]-2 Pi nd0[k]
d1[k_]=d[1,k]-2 Pi nd1[k]
d2[k_]=d[2,k]-2 Pi nd2[k]

Plot[{d[0,k],Pi/(2(Log[k]-.11593))},{k,0,.2},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d0[k],-(k+.25 Pi)},{k,0,4},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{d[1,k],-Pi (k/2)^2},{k,0,.8},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d1[k],-(k-.25 Pi)},{k,0,4},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{d[2,k],-Pi/2 (k/2)^4},{k,0,1},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d2[k],-(k-.75 Pi)},{k,0,5},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

f[phi_,k_,mmax_]=Sqrt[2/(Pi k I)]((Exp[2 I d[0,k]]-1)/2 +
Sum[(Exp[2 I d[m,k]]-1) Cos[m phi],{m,1,mmax}])

Plot[{Abs[f[phi,.1,1]]^2,Abs[f[phi,.1,2]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,.3,1]]^2,Abs[f[phi,.3,2]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,1.,3]]^2,Abs[f[phi,1.,6]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,3.,5]]^2,Abs[f[phi,3.,10]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,10.,15]]^2,Abs[f[phi,10.,30]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,.1,2]]^2,Abs[f[phi,.3,2]]^2,Abs[f[phi,1.,6]]^2,
Abs[f[phi,3.,10]]^2,Abs[f[phi,10.,15]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]}},PlotRange->{0,2.5}]


Plot[{Abs[f[phi,30.,50]]^2,Abs[f[phi,30.,100]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]}}]

Plot[.5 Sin[phi/2]+Cot[phi/2]^2 Sin[30 phi]^2/(2 Pi 30),{phi,0,Pi},
PlotStyle->{RGBColor[0,1,0]}]
-----
For[i=0,i<16,i++,dd[i]=d[i,10.]]

p0[r_,phi_]=Exp[I 10. r Cos[phi]]+(Exp[2 I dd[0]]-1)/2 (
BesselJ[0,10. r] + I BesselY[0,10. r])+
Sum[I^m (Exp[2 I dd[m]]-1) (BesselJ[m,10. r] + I BesselY[m,10. r])Cos[m phi],{m,1,15}]

p2[x_,y_]=If[x^2+y^2<1,0,Abs[p0[Sqrt[x^2+y^2],ArcTan[x,y]]]^2]

ContourPlot[p2[x,y],{x,-3,3},{y,-3,3},PlotPoints->100,ContourLines->False]
---
For[i=0,i<16,i++,dd[i]=d[i,1.]]

p0[r_,phi_]=Exp[I 1. r Cos[phi]]+(Exp[2 I dd[0]]-1)/2 (
BesselJ[0,1. r] + I BesselY[0,1. r])+
Sum[I^m (Exp[2 I dd[m]]-1) (BesselJ[m,1. r] + I BesselY[m,1. r])Cos[m phi],{m,1,4}]

p2[x_,y_]=If[x^2+y^2<1,0,Abs[p0[Sqrt[x^2+y^2],ArcTan[x,y]]]^2]

ContourPlot[p2[x,y],{x,-10,10},{y,-10,10},PlotPoints->100,ContourLines->False]
---
For[i=0,i<16,i++,dd[i]=d[i,.1]]

p0[r_,phi_]=Exp[I .1 r Cos[phi]]+(Exp[2 I dd[0]]-1)/2 (
BesselJ[0,.1 r] + I BesselY[0,.1 r])+
Sum[I^m (Exp[2 I dd[m]]-1) (BesselJ[m,.1 r] + I BesselY[m,.1 r])Cos[m phi],{m,1,2}]

p2[x_,y_]=If[x^2+y^2<1,0,Abs[p0[Sqrt[x^2+y^2],ArcTan[x,y]]]^2]

ContourPlot[p2[x,y],{x,-100,100},{y,-100,100},PlotPoints->100,ContourLines->False]

----
phi0[b_]=Sign[b](Pi-2 ArcSin[Abs[b]])
Plot[phi0[b],{b,-1,1}]

dsig[phi_]=.5 Sin[phi/2]
Plot[dsig[phi],{phi,0,Pi}]
---

For[i=0,i<16,i++,lnp[i]=besselJPrime[i,Sqrt[50]]/BesselJ[i,Sqrt[50]]]
Table[N[besselJPrime[i,Sqrt[50]]/BesselJ[i,Sqrt[50]]],{i,0,100}]
---
besselJPrime[m_,k_]=D[BesselJ[m,k],k]
besselYPrime[m_,k_]=D[BesselY[m,k],k]
u0=50

d2[m_,k_]=ArcTan[
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselY[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselYPrime[m,k]),
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselJ[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselJPrime[m,k])]

k0 =k /. N[FindRoot[-(Sqrt[k^2+u0](besselJPrime[0,Sqrt[k^2+u0]])BesselJ[0,k]
-BesselJ[0,Sqrt[k^2+u0]] k besselJPrime[0,k]),{k,6},AccuracyGoal->13],12]
k1 =k /. N[FindRoot[-(Sqrt[k^2+u0](besselJPrime[1,Sqrt[k^2+u0]])BesselJ[1,k]
-BesselJ[1,Sqrt[k^2+u0]] k besselJPrime[1,k]),{k,7},AccuracyGoal->13],12]
k2 =k /. N[FindRoot[-(Sqrt[k^2+u0](besselJPrime[2,Sqrt[k^2+u0]])BesselJ[2,k]
-BesselJ[2,Sqrt[k^2+u0]] k besselJPrime[2,k]),{k,5.5},AccuracyGoal->13],12]
k31 =k /. N[FindRoot[-(Sqrt[k^2+u0](besselJPrime[3,Sqrt[k^2+u0]])BesselJ[3,k]
-BesselJ[3,Sqrt[k^2+u0]] k besselJPrime[3,k]),{k,3.5},AccuracyGoal->13],12]
k32 =k /. N[FindRoot[-(Sqrt[k^2+u0](besselJPrime[3,Sqrt[k^2+u0]])BesselJ[3,k]
-BesselJ[3,Sqrt[k^2+u0]] k besselJPrime[3,k]),{k,6},AccuracyGoal->13],12]

fd0[k_]=If[k<k0,d2[0,k]+2 Pi,d2[0,k]]
fd1[k_]=If[k<k1,d2[1,k]+2 Pi,d2[1,k]]
fd2[k_]=If[k<k2,d2[2,k]+2 Pi,d2[2,k]]
fd3[k_]=If[(k-k31)*(k-k32)<0,d2[3,k]+2 Pi,d2[3,k]]

 Plot[{d2[3,k],d2[4,k],d2[5,k]},{k,.1,30},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{fd0[k],fd1[k],fd2[k],fd3[k],d2[4,k],d2[5,k]},{k,.1,30},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]}},PlotRange->{0,10}]

Plot[{d2[6,k],d2[7,k],d2[8,k],d2[9,k],d2[10,k]},{k,1,10},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]}},PlotRange->{0,10}]


----
For[i=0,i<21,i++,dd[i]=d2[i,10.]]
For[i=0,i<21,i++,n[i]=.5((Exp[2 I dd[i]]+1)BesselJ[i,10.] + 
I (Exp[2 I dd[i]]-1)BesselY[i,10.] )/BesselJ[i,Sqrt[150.]] ]

p1[r_,phi_]= n[0] BesselJ[0,Sqrt[150.] r] + 
Sum[2 I^m n[m] (BesselJ[m,Sqrt[150.] r] )Cos[m phi],{m,1,20}]

p0[r_,phi_]=Exp[I 10. r Cos[phi]]+(Exp[2 I dd[0]]-1)/2 (
BesselJ[0,10. r] + I BesselY[0,10. r])+
Sum[I^m (Exp[2 I dd[m]]-1) (BesselJ[m,10. r] + I BesselY[m,10. r])Cos[m phi],{m,1,15}]

p2[x_,y_]=If[x^2+y^2<1,Abs[p1[Sqrt[x^2+y^2],ArcTan[x,y]]]^2,Abs[p0[Sqrt[x^2+y^2],ArcTan[x,y]]]^2]

ContourPlot[p2[x,y],{x,-3,3},{y,-3,3},PlotPoints->100,ContourLines->False]

Plot[p2[x,0],{x,-3,10}]
Plot[p2[x,.5],{x,-3,3}]
----
For[i=0,i<21,i++,dd[i]=d2[i,1.]]
For[i=0,i<21,i++,n[i]=.5((Exp[2 I dd[i]]+1)BesselJ[i,1.] + 
I (Exp[2 I dd[i]]-1)BesselY[i,1.] )/BesselJ[i,Sqrt[51.]] ]

p1[r_,phi_]= n[0] BesselJ[0,Sqrt[51.] r] + 
Sum[2 I^m n[m] (BesselJ[m,Sqrt[51.] r] )Cos[m phi],{m,1,8}]

p0[r_,phi_]=Exp[I  r Cos[phi]]+(Exp[2 I dd[0]]-1)/2 (
BesselJ[0,1. r] + I BesselY[0,1. r])+
Sum[I^m (Exp[2 I dd[m]]-1) (BesselJ[m,1. r] + I BesselY[m,1. r])Cos[m phi],{m,1,5}]

p2[x_,y_]=If[x^2+y^2<1,Abs[p1[Sqrt[x^2+y^2],ArcTan[x,y]]]^2,Abs[p0[Sqrt[x^2+y^2],ArcTan[x,y]]]^2]

ContourPlot[p2[x,y],{x,-10,10},{y,-10,10},PlotPoints->100,ContourLines->False]
ContourPlot[p2[x,y],{x,-3,3},{y,-3,3},PlotPoints->100,ContourLines->False]

Plot[p2[x,0],{x,-10,10}]
Plot[p2[x,.5],{x,-10,10}]
----
kr2 =k /. N[FindRoot[d2[5,k]-Pi/2,{k,2.4},AccuracyGoal->12],11]
2.3987659659
kr0 =k /. N[FindRoot[d2[0,k]-Pi/2,{k,.001},AccuracyGoal->12],11]
0.086138291913
---
dd[5]=d2[5,2.4]
n[5]=.5((Exp[2 I dd[5]]+1)BesselJ[5,2.4] + 
I (Exp[2 I dd[5]]-1)BesselY[5,2.4] )/BesselJ[5,Sqrt[2.4^2+50]]
p5[r_]=If[r<1,Abs[n[5]]^2 BesselJ[5,Sqrt[2.4^2+50]r],Re[Conjugate[n[5]] .5((Exp[2 I dd[5]]+1)BesselJ[5,2.4 r] + 
I (Exp[2 I dd[5]]-1)BesselY[5,2.4 r] )]]
v5[r_]=If[r<1,-50+25/r^2,25/r^2]
Plot[{v5[r],.2 p5[r]+2.4^2,2.4^2},{r,0,5},PlotStyle->{{RGBColor[1,0,0]},
 {RGBColor[0,0,0]}, {RGBColor[1,0,1]}},PlotRange->{-30,30}]
-
dd[5]=d2[5,2.5]
n[5]=.5((Exp[2 I dd[5]]+1)BesselJ[5,2.5] + 
I (Exp[2 I dd[5]]-1)BesselY[5,2.5] )/BesselJ[5,Sqrt[2.5^2+50]]
p5[r_]=If[r<1,Abs[n[5]]^2 BesselJ[5,Sqrt[2.5^2+50]r],Re[Conjugate[n[5]] .5((Exp[2 I dd[5]]+1)BesselJ[5,2.5 r] + 
I (Exp[2 I dd[5]]-1)BesselY[5,2.5 r] )]]
v5[r_]=If[r<1,-50+25/r^2,25/r^2]
Plot[{v5[r],2 p5[r]+2.5^2,2.5^2},{r,0,5},PlotStyle->{{RGBColor[1,0,0]},
 {RGBColor[0,0,0]}, {RGBColor[1,0,1]}},PlotRange->{-30,30}]
-----

f[phi_,k_,mmax_]=Sqrt[2/(Pi k I)]((Exp[2 I d2[0,k]]-1)/2 +
Sum[(Exp[2 I d2[m,k]]-1) Cos[m phi],{m,1,mmax}])

Plot[{Abs[f[phi,.086,5]]^2,Abs[f[phi,.086,3]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,1,5]]^2,Abs[f[phi,1,3]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,2.4,5]]^2,Abs[f[phi,2.4,6]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,3,5]]^2,Abs[f[phi,3,6]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,10,15]]^2,Abs[f[phi,10,20]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,30,45]]^2,Abs[f[phi,30,60]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[phi,2.3,5]]^2,Abs[f[phi,2.4,5]]^2,Abs[f[phi,2.5,5]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]



Plot[{Abs[f[phi,.1,4]]^2,Abs[f[phi,.3,4]]^2,Abs[f[phi,1.,4]]^2,
Abs[f[phi,3.,6]]^2,Abs[f[phi,10.,15]]^2},{phi,0,Pi},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]}},PlotRange->{0,7}]

---
sp[i_,n_,b_]= (Pi-2ArcSin[b])-i(Pi-2ArcSin[b/n]) 
In[168]:= FindMinimum[-sp[2,Sqrt[1.5],b],{b,.8}]

Out[168]= {2.07784, {b -> 0.912871}}
Plot[{sp[1,Sqrt[1.5],b],sp[2,Sqrt[1.5],b]},{b,0,1}]
Plot[{sp[1,Sqrt[1+50/9],b],sp[2,Sqrt[1+50/9],b],sp[3,Sqrt[1+50/9],b]},{b,0,1}]
FindMinimum[-sp[3,Sqrt[1+50/9],b],{b,.5}]
Plot[{sp[0,1.33,b],sp[1,1.33,b],sp[2,1.33,b],sp[3,1.33,b]},{b,0,1},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]
FindMinimum[-sp[2,1.33,b],{b,.95}]
FindMinimum[-sp[3,1.33,b],{b,.95}]
FindMinimum[-sp[4,1.33,b],{b,.98}]
---

R[n_,theta_]=Abs[(Cos[theta]-Sqrt[n^2-Sin[theta]^2])/(Cos[theta]+Sqrt[n^2-Sin[theta]^2])]^2
sp[i_,b_]= (Pi-2ArcSin[b])-i(Pi-2ArcSin[b/1.33]) 
dsp[i_,b_]= Abs[1/D[sp[i,b],b]]

ParametricPlot[{{Pi-Abs[Mod[sp[0,b],2 Pi]-Pi],R[1.33,ArcSin[b]] dsp[0,b]},
{Pi-Abs[Mod[sp[1,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])(1-R[1/1.33,ArcSin[b/1.33]])dsp[1,b]},
{Pi-Abs[Mod[sp[2,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])R[1/1.33,ArcSin[b/1.33]](1-R[1/1.33,ArcSin[b/1.33]])dsp[2,b]},
{Pi-Abs[Mod[sp[3,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])R[1/1.33,ArcSin[b/1.33]]^2(1-R[1/1.33,ArcSin[b/1.33]])dsp[3,b]}},
{b,0,1},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}},
PlotRange->{0,1}]

ParametricPlot[{{Pi-Abs[Mod[sp[0,b],2 Pi]-Pi],R[1.33,ArcSin[b]] dsp[0,b]},
{Pi-Abs[Mod[sp[1,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])(1-R[1/1.33,ArcSin[b/1.33]])dsp[1,b]},
{Pi-Abs[Mod[sp[2,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])R[1/1.33,ArcSin[b/1.33]](1-R[1/1.33,ArcSin[b/1.33]])dsp[2,b]},
{Pi-Abs[Mod[sp[3,b],2 Pi]-Pi],(1-R[1.33,ArcSin[b]])R[1/1.33,ArcSin[b/1.33]]^2(1-R[1/1.33,ArcSin[b/1.33]])dsp[3,b]}},
{b,0,1},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}},
PlotRange->{{2,3},{0,.45}}]
---
b=0.86238
theta=ArcSin[b]
alpha=ArcSin[b/1.33]
phi=2.39954-Pi
Graphics[{{RGBColor[1,0,0],Circle[{0,0},1]},Line[{{2.5,b},{Cos[theta],Sin[theta]},
{Cos[theta+Pi-2 alpha],Sin[theta+Pi-2 alpha]},
{Cos[theta+2Pi-4 alpha],Sin[theta+2Pi-4 alpha]},
{Cos[theta+2Pi-4 alpha]+1.5 Cos[phi],Sin[theta+2Pi-4 alpha]+1.5 Sin[phi]}}]} ]
Show[%,AspectRatio->Automatic]
---
b=0.95073
theta=ArcSin[b]
alpha=ArcSin[b/1.33]
phi=4.01602-Pi
Graphics[{{RGBColor[1,0,0],Circle[{0,0},1]},Line[{{2.5,b},{Cos[theta],Sin[theta]},
{Cos[theta+Pi-2 alpha],Sin[theta+Pi-2 alpha]},
{Cos[theta+2Pi-4 alpha],Sin[theta+2Pi-4 alpha]},
{Cos[theta+3Pi-6 alpha],Sin[theta+3Pi-6 alpha]},
{Cos[theta+3Pi-6 alpha]+1.5 Cos[phi],Sin[theta+3Pi-6 alpha]+1.5 Sin[phi]}}]} ]
Show[%,AspectRatio->Automatic]

---
u0=700
d2[m_,k_]=ArcTan[
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselY[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselYPrime[m,k]),
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselJ[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselJPrime[m,k])]
For[i=0,i<51,i++,dd[i]=d2[i,30.]]
f30[phi_]=Sqrt[2/(Pi 30. I)]((Exp[2 I dd[0]]-1)/2 +
Sum[(Exp[2 I dd[m]]-1) Cos[m phi],{m,1,45}])
Plot[Abs[f30[phi]]^2,{phi,0,Pi}]
Plot[Abs[f30[phi]]^2,{phi,2,3},PlotPoints->200,PlotRange->{0,.45}]

---
u0=7700
d2[m_,k_]=ArcTan[
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselY[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselYPrime[m,k]),
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselJ[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselJPrime[m,k])]
For[i=0,i<151,i++,dd[i]=d2[i,100.]]
f100[phi_]=Sqrt[2/(Pi 100. I)]((Exp[2 I dd[0]]-1)/2 +
Sum[(Exp[2 I dd[m]]-1) Cos[m phi],{m,1,150}])
Plot[Abs[f100[phi]]^2,{phi,0,Pi},PlotPoints->200,PlotRange->{0,1}]
Plot[Abs[f100[phi]]^2,{phi,2,3},PlotPoints->200]
n=Sqrt[1+u0/100^2]
sp[i_,b_]= (Pi-2ArcSin[b])-i(Pi-2ArcSin[b/n]) 
dsp[i_,b_]=1/Abs[D[sp[i,b],b]] 
ParametricPlot[{{Mod[sp[i,b],2 Pi],dsp[1,b]},{Mod[sp[2,b],2 Pi],dsp[2,b]},{Mod[sp[3,b],2 Pi],dsp[3,b]},
{Mod[sp[4,b],2 Pi],dsp[4,b]},{Mod[sp[5,b],2 Pi],dsp[5,b]}},{b,0,1}]

ParametricPlot[{{Pi-Abs[Mod[sp[i,b],2 Pi]-Pi],dsp[1,b]},{Pi-Abs[Mod[sp[2,b],2 Pi]-Pi],dsp[2,b]},{Pi-Abs[Mod[sp[3,b],2 Pi]-Pi],dsp[3,b]},
{Pi-Abs[Mod[sp[4,b],2 Pi]-Pi],dsp[4,b]},{Pi-Abs[Mod[sp[5,b],2 Pi]-Pi],dsp[5,b]},
{Pi-Abs[Mod[sp[6,b],2 Pi]-Pi],dsp[4,b]},{Pi-Abs[Mod[sp[7,b],2 Pi]-Pi],dsp[5,b]}},{b,0,1}]

---
u0=69000
d2[m_,k_]=ArcTan[
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselY[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselYPrime[m,k]),
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselJ[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselJPrime[m,k])]
For[i=0,i<451,i++,dd[i]=d2[i,300.]]
f300[phi_]=Sqrt[2/(Pi 300. I)]((Exp[2 I dd[0]]-1)/2 +
Sum[(Exp[2 I dd[m]]-1) Cos[m phi],{m,1,330}])
Plot[Abs[f300[phi]]^2,{phi,0,Pi},PlotPoints->200,PlotRange->{0,1}]
Plot[Abs[f300[phi]]^2,{phi,2,3},PlotPoints->200]
---
besselJPrime[m_,k_]=D[BesselJ[m,k],k]
besselYPrime[m_,k_]=D[BesselY[m,k],k]
u0=50

d2[m_,k_]=ArcTan[
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselY[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselYPrime[m,k]),
-(Sqrt[k^2+u0](besselJPrime[m,Sqrt[k^2+u0]])BesselJ[m,k]
-BesselJ[m,Sqrt[k^2+u0]] k besselJPrime[m,k])]

sigma[k_]=(4/k)(Sin[d2[0,k]]^2 + 2 Sum[Sin[d2[m,k]]^2,{m,1,11}] )

Plot[sigma[k],{k,.01,10},PlotPoints->100,PlotRange->{0,14},Frame->True]

sigma[k_]=(4/k)(Sin[d2[0,k]]^2 + 2 Sum[Sin[d2[m,k]]^2,{m,1,45}] )

Plot[sigma[k],{k,10,30},PlotRange->{0,14},Frame->True]

 Show[GraphicsArray[{%23,%24}]]
---
points={{RGBColor[1,0,0],Point[{1,1,1}]},{RGBColor[1,117/256,0],Point[{2,1,1}]},
{RGBColor[1,117/256,0],Point[{1,2,1}]},{RGBColor[1,117/256,0],Point[{1,1,2}]},
{RGBColor[0,176/256,0],Point[{2,2,1}]},{RGBColor[0,176/256,0],Point[{2,1,2}]},
{RGBColor[0,176/256,0],Point[{1,2,2}]},
{RGBColor[0,117/256,117/256],Point[{3,1,1}]},
{RGBColor[0,117/256,117/256],Point[{1,3,1}]},{RGBColor[0,117/256,117/256],Point[{1,1,3}]},
{RGBColor[0,0,1],Point[{2,2,2}]},{RGBColor[.5,0,1],Point[{1,2,3}]},
{RGBColor[.5,0,1],Point[{2,3,1}]},{RGBColor[.5,0,1],Point[{3,1,2}]},
{RGBColor[.5,0,1],Point[{2,1,3}]},{RGBColor[.5,0,1],Point[{3,2,1}]},
{RGBColor[.5,0,1],Point[{1,3,2}]},
{Point[{3,3,3}],Point[{2,3,3}],Point[{3,2,3}],Point[{3,3,2}],
Point[{1,3,3}],Point[{3,1,3}],Point[{3,3,1}],
Point[{2,2,3}],Point[{2,3,2}],Point[{3,2,2}] }}

}
Show[Graphics3D[{PointSize[.05],points}],FaceGrids->All,PlotRange->{{0,3.1},{0,3.1},{0,3.1}},
ViewPoint->{.3,-4,1},Axes->True,AxesStyle->Thickness[.01],
AxesEdge->{{-1,-1},{-1,-1},{-1,-1}},AxesLabel->{"x","y","z"}]

---
<<NumericalMath`BesselZeros`

g0=BesselJZeros[1/2,5]
g1=BesselJZeros[3/2,5]
g2=BesselJZeros[5/2,5]
g3=BesselJZeros[7/2,5]
g4=BesselJZeros[9/2,5]
g5=BesselJZeros[11/2,5]
g6=BesselJZeros[13/2,5]

Plot[{g0[[1]]^2,g0[[2]]^2,g0[[3]]^2},{x,-.15,.15},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{g1[[1]]^2,g1[[2]]^2},{x,.85,1.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{g2[[1]]^2,g2[[2]]^2},{x,1.85,2.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{g3[[1]]^2},{x,2.85,3.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{g4[[1]]^2},{x,3.85,4.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{g5[[1]]^2},{x,4.85,5.15},PlotStyle->{{RGBColor[1,0,0]}}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,PlotRange->{0,100}]
---
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
sBessely[l_,x_]=Sqrt[Pi/(2 x)] BesselY[l+(1/2),x]
sBesseli[l_,x_]=Sqrt[Pi/(2 x)] BesselI[l+(1/2),x]
sBesselk[l_,x_]=Sqrt[Pi/(2 x)] BesselK[l+(1/2),x]

Plot[{sBesselj[0,x],sBesselj[1,x],sBesselj[2,x],sBesselj[3,x]},{x,0,10},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{sBessely[0,x],sBessely[1,x],sBessely[2,x],sBessely[3,x]},{x,0,10},
PlotRange->{-2,.5},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{sBesseli[0,x],sBesseli[1,x],sBesseli[2,x],sBesseli[3,x]},{x,0,5},
PlotRange->{0,10},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]

Plot[{sBesselk[0,x],sBesselk[1,x],sBesselk[2,x],sBesselk[3,x]},{x,0,3},
PlotRange->{0,10},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]
<<Graphics`Master`

LogPlot[{sBesselk[0,x],sBesselk[1,x],sBesselk[2,x],sBesselk[3,x]},{x,0,5},
PlotRange->{.01,100},PlotPoints->100,
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,.5,0]},{RGBColor[0,1,0]},{RGBColor[0,0,1]}}]
---
<<NumericalMath`BesselZeros`
roots={BesselJZeros[1/2,11],BesselJZeros[3/2,11],BesselJZeros[5/2,11],BesselJZeros[7/2,11],BesselJZeros[9/2,11],
BesselJZeros[11/2,11],BesselJZeros[13/2,11],BesselJZeros[15/2,11],BesselJZeros[17/2,11],BesselJZeros[19/2,11],
BesselJZeros[21/2,11]}
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
v[l_,r_]=If[r<1,l(l+1)/r^2,231]
k00=roots[[1,1]]
k01=roots[[1,2]]
k02=roots[[1,3]]
k03=roots[[1,4]]
r00[r_]=sBesselj[0,k00 r] Sqrt[2./(sBesselj[1,k00])^2]
r01[r_]=sBesselj[0,k01 r] Sqrt[2./(sBesselj[1,k01])^2]
r02[r_]=sBesselj[0,k02 r] Sqrt[2./(sBesselj[1,k02])^2]
r03[r_]=sBesselj[0,k03 r] Sqrt[2./(sBesselj[1,k03])^2]

a=3
Plot[{a r00[r]+k00^2,a r01[r]+k01^2,a r02[r]+k02^2,a r03[r]+k03^2},{r,0,1},PlotRange->{0,230}]
Plot[{k00^2,k01^2,k02^2,k03^2,v[0,r]},{r,0,1.1},PlotRange->{0,230}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

k10=roots[[2,1]]
k11=roots[[2,2]]
k12=roots[[2,3]]
k13=roots[[2,4]]
r10[r_]=sBesselj[1,k10 r] Sqrt[2./(sBesselj[2,k10]^2)]
r11[r_]=sBesselj[1,k11 r] Sqrt[2./(sBesselj[2,k11]^2)]
r12[r_]=sBesselj[1,k12 r] Sqrt[2./(sBesselj[2,k12]^2)]
r13[r_]=sBesselj[1,k13 r] Sqrt[2./(sBesselj[2,k13]^2)]

a=3
Plot[{a r10[r]+k10^2,a r11[r]+k11^2,a r12[r]+k12^2,a r13[r]+k13^2},{r,0,1},PlotRange->{0,230}]
Plot[{k10^2,k11^2,k12^2,k13^2,v[1,r]},{r,0,1.1},PlotRange->{0,230}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

k20=roots[[3,1]]
k21=roots[[3,2]]
k22=roots[[3,3]]
r20[r_]=sBesselj[2,k20 r] Sqrt[2./(sBesselj[3,k20]^2)]
r21[r_]=sBesselj[2,k21 r] Sqrt[2./(sBesselj[3,k21]^2)]
r22[r_]=sBesselj[2,k22 r] Sqrt[2./(sBesselj[3,k22]^2)]

a=3
Plot[{a r20[r]+k20^2,a r21[r]+k21^2,a r22[r]+k22^2},{r,0,1},PlotRange->{0,230}]
Plot[{k20^2,k21^2,k22^2,v[2,r]},{r,0,1.1},PlotRange->{0,230}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio]

Show[GraphicsArray[{%49,%52,%55}]]
----
p12[r_]=If[r<1,r12[r],0]
m=0
Plot[Abs[r12[Sqrt[0^2+z^2]] SphericalHarmonicY[1,m,ArcTan[z,0],0] ]^2,{z,0,1}]
ContourPlot[Abs[p12[Sqrt[x^2+z^2]] SphericalHarmonicY[1,0,ArcTan[z,x],0] ]^2]
{x,-1,1},{z,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,5}]
m=1
Plot[Abs[r12[x] SphericalHarmonicY[1,1,Pi/2,0] ]^2,{x,0,1}]
ContourPlot[Abs[p12[Sqrt[x^2+z^2]] SphericalHarmonicY[1,1,ArcTan[z,x],0] ]^2,
{x,-1,1},{z,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,2.5}]
Plot[Abs[r12[r] SphericalHarmonicY[1,1,Pi/2,0] ]^2,{r,0,1}]
ContourPlot[Abs[p12[Sqrt[x^2+y^2]] SphericalHarmonicY[1,1,Pi/2,0] ]^2,
{x,-1,1},{y,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,2.5}]

p22[r_]=If[r<1,r22[r],0]
m=0
Plot[Abs[r22[z] SphericalHarmonicY[2,0,0,0] ]^2,{z,0,1},PlotRange->{0,12}]
ContourPlot[Abs[p22[Sqrt[x^2+z^2]] SphericalHarmonicY[2,0,ArcTan[z,x],0] ]^2,
{x,-1,1},{z,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,5}]
m=1
ContourPlot[Abs[p22[Sqrt[x^2+z^2]] SphericalHarmonicY[2,1,ArcTan[z,x],0] ]^2,
{x,-1,1},{z,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,5}]

Plot[Abs[r22[r] SphericalHarmonicY[2,2,Pi/2,0] ]^2,{r,0,1},PlotRange->{0,5}]
ContourPlot[Abs[p22[Sqrt[x^2+z^2]] SphericalHarmonicY[2,2,ArcTan[z,x],0] ]^2,
{x,-1,1},{z,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,2.5}]
ContourPlot[Abs[p22[Sqrt[x^2+y^2]] SphericalHarmonicY[2,2,Pi/2,0] ]^2,
{x,-1,1},{y,-1,1},PlotPoints->100,
ContourLines->False,PlotRange->{0,2.5}]

---
k010=roots[[1,11]]
r010[r_]=sBesselj[0,k010 r] Sqrt[2./(sBesselj[1,k010]^2)]
Plot[r r010[r],{r,0,1}]

k105=roots[[11,6]]
r105[r_]=sBesselj[10,k105 r] Sqrt[2./(sBesselj[11,k105]^2)]
Plot[r r105[r],{r,0.1,1}]
----

u0=50
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
sBessely[l_,x_]=Sqrt[Pi/(2 x)] BesselY[l+(1/2),x]
sBesseli[l_,x_]=Sqrt[Pi/(2 x)] BesselI[l+(1/2),x]
sBesselk[l_,x_]=Sqrt[Pi/(2 x)] BesselK[l+(1/2),x]

sbesseljPrime[m_,k_]=D[sBesselj[m,k],k]
sbesselkPrime[m_,k_]=D[sBesselk[m,k],k]

f[l_,k_]=k sbesseljPrime[l,k]/sBesselj[l,k]-Sqrt[u0-k^2] sbesselkPrime[l, Sqrt[u0-k^2]]/sBesselk[l,Sqrt[u0-k^2]]
Plot[f[1,x],{x,0,Sqrt[50]}]
Plot[f[0,x],{x,7,Sqrt[50]}]
Plot[f[5,x],{x,0,Sqrt[50]}]

Plot[k sbesseljPrime[1,k]/sBesselj[1,k],{k,0,8},PlotStyle->{RGBColor[1,0,0]},
PlotRange->{-10,5}]
Plot[Sqrt[u0-k^2] sbesselkPrime[1, Sqrt[u0-k^2]]/sBesselk[1,Sqrt[u0-k^2]],
{k,0,Sqrt[50]},PlotStyle->{RGBColor[0,0,1]},
PlotRange->{-10,5}]
Show[%,%%]


k[0,0]=x /. N[FindRoot[f[0,x],{x,2.7},AccuracyGoal->13],12]
k[0,1]=x /. N[FindRoot[f[0,x],{x,5.4},AccuracyGoal->13],12]
k[1,0]=x /. N[FindRoot[f[1,x],{x,3.9},AccuracyGoal->13],12]
k[1,1]=x /. N[FindRoot[f[1,x],{x,6.5},AccuracyGoal->13],12]
k[2,0]=x /. N[FindRoot[f[2,x],{x,5},AccuracyGoal->13],12]
k[3,0]=x /. N[FindRoot[f[3,x],{x,6},AccuracyGoal->13],12]
k[4,0]=x /. N[FindRoot[f[4,x],{x,7},AccuracyGoal->13],12]

un[l_,nr_,r_]=If[r<1,sBesselj[l,k[l,nr] r], sBesselj[l,k[l,nr]] sBesselk[l,Sqrt[u0-k[l,nr]^2] r]/sBesselk[l,Sqrt[u0-k[l,nr]^2]]]
Plot[un[0,0,x],{x,0,1.4}]
Plot[un[0,0,x],{x,.95,1.05}]

Plot[un[2,1,x],{x,.9,1.1}]


norm[0,0]=1/Sqrt[Integrate[sBesselj[0,k[0,0] r]^2 r^2,{r,0,1}]+ 
(sBesselj[0,k[0,0]]/sBesselk[0,Sqrt[u0-k[0,0]^2]])^2 Integrate[sBesselk[0,Sqrt[u0-k[0,0]^2] r]^2 r^2,{r,1,Infinity}]]
norm[0,1]=1/Sqrt[Integrate[sBesselj[0,k[0,1] r]^2 r^2,{r,0,1}]+ 
(sBesselj[0,k[0,1]]/sBesselk[0,Sqrt[u0-k[0,1]^2]])^2 Integrate[sBesselk[0,Sqrt[u0-k[0,1]^2] r]^2 r^2,{r,1,Infinity}]]
%3.61222, 6.9297


norm[1,0]=1/Sqrt[NIntegrate[sBesselj[1,k[1,0] r]^2 r^2,{r,0,1}]+ 
(sBesselj[1,k[1,0]]/sBesselk[1,Sqrt[u0-k[1,0]^2]])^2 Integrate[sBesselk[1,Sqrt[u0-k[1,0]^2] r]^2 r^2,{r,1,Infinity}]]
norm[1,1]=1/Sqrt[NIntegrate[sBesselj[1,k[1,1] r]^2 r^2,{r,0,1}]+ 
(sBesselj[1,k[1,1]]/sBesselk[1,Sqrt[u0-k[1,1]^2]])^2 Integrate[sBesselk[1,Sqrt[u0-k[1,1]^2] r]^2 r^2,{r,1,Infinity}]]
%8.03135, 6.9297

norm[2,0]=1/Sqrt[NIntegrate[sBesselj[2,k[2,0] r]^2 r^2,{r,0,1}]+ 
(sBesselj[2,k[2,0]]/sBesselk[2,Sqrt[u0-k[2,0]^2]])^2 NIntegrate[sBesselk[2,Sqrt[u0-k[2,0]^2] r]^2 r^2,{r,1,5}]]

norm[3,0]=1/Sqrt[NIntegrate[sBesselj[3,k[3,0] r]^2 r^2,{r,0,1}]+ 
(sBesselj[3,k[3,0]]/sBesselk[3,Sqrt[u0-k[3,0]^2]])^2 NIntegrate[sBesselk[3,Sqrt[u0-k[3,0]^2] r]^2 r^2,{r,1,5}]]

norm[4,0]=1/Sqrt[NIntegrate[sBesselj[4,k[4,0] r]^2 r^2,{r,0,1}]+ 
(sBesselj[4,k[4,0]]/sBesselk[4,Sqrt[u0-k[4,0]^2]])^2 NIntegrate[sBesselk[4,Sqrt[u0-k[4,0]^2] r]^2 r^2,{r,1,15}]]

psi[l_,nr_,r_]=norm[l,nr] If[r<1,sBesselj[l,k[l,nr] r], sBesselj[l,k[l,nr]] sBesselk[l,Sqrt[u0-k[l,nr]^2] r]/sBesselk[l,Sqrt[u0-k[l,nr]^2]]]

NIntegrate[psi[3,0,r]^2 r^2,{r,0,1.5}]
v[l_,r_]=If[r<1,l(l+1)/r^2,l(l+1)/r^2+50]

a=3
Plot[{a psi[0,0,r]+k[0,0]^2,a psi[0,1,r]+k[0,1]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[0,0]^2,k[0,1]^2,v[0,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50,60}}]

Plot[{a psi[1,0,r]+k[1,0]^2,a psi[1,1,r]+k[1,1]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[1,0]^2,k[1,1]^2,v[1,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50,60}}]

Plot[{a psi[2,0,r]+k[2,0]^2},{r,0,1.5},PlotRange->{0,60}]
Plot[{k[2,0]^2,v[2,r]},{r,0,1.5},PlotRange->{0,60}, \
PlotStyle->{{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]
Show[%,%%,AspectRatio->GoldenRatio,Ticks->{{0,.5,1,1.5},{10,20,30,40,50,60}}]

Show[GraphicsArray[{%67,%70,%73}]]

---

<<NumericalMath`BesselZeros`

g0=BesselJZeros[1/2,5]
g1=BesselJZeros[3/2,5]
g2=BesselJZeros[5/2,5]
g3=BesselJZeros[7/2,5]
g4=BesselJZeros[9/2,5]
g5=BesselJZeros[11/2,5]
g6=BesselJZeros[13/2,5]

Plot[{g0[[1]]^2,g0[[2]]^2,g0[[3]]^2},{x,-.15,.15}]
Plot[{g1[[1]]^2,g1[[2]]^2},{x,.85,1.15}]
Plot[{g2[[1]]^2,g2[[2]]^2},{x,1.85,2.15}]
Plot[{g3[[1]]^2},{x,2.85,3.15}]
Plot[{g4[[1]]^2},{x,3.85,4.15}]
Plot[{g5[[1]]^2},{x,4.85,5.15}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%]
Plot[{k[0,0]^2,k[0,1]^2},{x,-.15,.15},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[1,0]^2,k[1,1]^2},{x,.85,1.15},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[{k[2,0]^2},{x,1.85,2.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[3,0]^2},{x,2.85,3.15},PlotStyle->{{RGBColor[1,0,0]}}]
Plot[{k[4,0]^2},{x,3.85,4.15},PlotStyle->{{RGBColor[1,0,0]}}]
Show[%,%%,%%%,%%%%,%%%%%,PlotRange->{0,100}]

Show[%,%99]

-----
<<Graphics`Shapes`
In[5]:= ?RotateShape
RotateShape[graphics3D, phi, theta, psi] rotates the three-dimensional
   graphics object by the specified Euler angles.

<<Graphics`ParametricPlot3D`
surf=SphericalPlot3D[{1,GrayLevel[If[Abs[theta-Pi/4]<.05,0,1]]},
{theta, 0, Pi},{phi,0, 2 Pi},Lighting->False,PlotPoints->30]

Show[RotateShape[surf,0,Pi/2,Pi/2],Lighting->False,Axes->False]

ParametricPlot3D[{Sin[Pi/4-.15]Cos[p],Sin[Pi/4-.15]Sin[p],Cos[Pi/4-.15],RGBColor[1,0,0]},{p,0,2 Pi}]
ParametricPlot3D[{Sin[Pi/4-.16]Cos[p],Sin[Pi/4-.16]Sin[p],Cos[Pi/4-.16],RGBColor[1,0,0]},{p,0,2 Pi}]
ParametricPlot3D[{Sin[Pi/4-.155]Cos[p],Sin[Pi/4-.155]Sin[p],Cos[Pi/4-.155],RGBColor[1,0,0]},{p,0,2 Pi}]

Show[{%7,RotateShape[%%,0,Pi/2,Pi/2],RotateShape[%%%,0,Pi/2,Pi/2],RotateShape[%,0,Pi/2,Pi/2]},Lighting->False,Axes->False]

---
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
sBessely[l_,x_]=Sqrt[Pi/(2 x)] BesselY[l+(1/2),x]
<<NumericalMath`BesselZeros`

g0=BesselJZeros[1/2,5]
nd0[k_]=If[k<g0[[1]],0,If[k<g0[[3]],1,If[k<g0[[5]],2,3]]]
g1=BesselJZeros[3/2,5]
nd1[k_]=If[k<g1[[1]],0,If[k<g1[[3]],1,If[k<g1[[5]],2,3]]]
g2=BesselJZeros[5/2,5]
nd2[k_]=If[k<g2[[1]],0,If[k<g2[[3]],1,If[k<g2[[5]],2,3]]]

d[l_,k_]=ArcTan[-sBessely[l,k],-sBesselj[l,k]]
d0[k_]=d[0,k]-2 Pi nd0[k]
d1[k_]=d[1,k]-2 Pi nd1[k]
d2[k_]=d[2,k]-2 Pi nd2[k]

Plot[{d[0,k],-2(k/2)},{k,0,.3},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d0[k],-(k)},{k,0,4},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{d[1,k],-(8/3) (k/2)^3},{k,0,.8},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d1[k],-(k- Pi/2)},{k,0,4},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{d[2,k],-(32/45) (k/2)^5},{k,0,1},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]
Plot[{d[2,k],-(k- Pi)},{k,0,4},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Show[%12,%8,%11,PlotRange->{-.3,.02},Frame->True]
Show[%41,%42,%43,PlotRange->{-4,Pi}]

f[theta_,k_,lmax_]=1/(2 k I) (Sum[(Exp[2 I d[l,k]]-1) (2l+1)
LegendreP[l, Cos[theta]],{l,0,lmax}])

Plot[{Abs[f[theta,.1,1]]^2,Abs[f[theta,.1,2]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,.3,1]]^2,Abs[f[theta,.3,2]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,1.,3]]^2,Abs[f[theta,1.,6]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,3.,5]]^2,Abs[f[theta,3.,10]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,10.,15]]^2,Abs[f[theta,10.,30]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,.1,2]]^2,Abs[f[theta,.3,2]]^2,Abs[f[theta,1.,3]]^2,
Abs[f[theta,3.,6]]^2,Abs[f[theta,10.,15]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]}},PlotRange->{0,1.5}]


Plot[{Abs[f[theta,30.,45]]^2,Abs[f[theta,30.,60]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]}}]

Plot[.25(1+Cot[theta/2]^2 BesselJ[1,30 Sin[theta]]^2),{theta,0,Pi},
PlotStyle->{RGBColor[0,1,0]}]


Plot[.5 Sin[theta/2]+Cot[theta/2]^2 Sin[30 theta]^2/(2 Pi 30),{theta,0,Pi},
PlotStyle->{RGBColor[0,1,0]}]
-----
For[i=0,i<16,i++,dd[i]=d[i,10.]]

p0[r_,theta_]=Exp[I 10. r Cos[theta]]+
Sum[I^l (2l+1) (Exp[2 I dd[l]]-1) (sBesselj[l,10. r] + I sBessely[l,10. r])
LegendreP[l, Cos[theta]],{l,0,15}]/2

p2[x_,z_]=If[x^2+z^2<1,0,Abs[p0[Sqrt[x^2+z^2],ArcTan[z,x]]]^2]

ContourPlot[p2[x,z],{z,-3,3},{x,-3,3},PlotPoints->100,ContourLines->False]
---
For[i=0,i<16,i++,dd[i]=d[i,1.]]
p0[r_,theta_]=Exp[I 1. r Cos[theta]]+
Sum[I^l (2l+1) (Exp[2 I dd[l]]-1) (sBesselj[l,1. r] + I sBessely[l,1. r])
LegendreP[l, Cos[theta]],{l,0,4}]/2

p2[x_,z_]=If[x^2+z^2<1,0,Abs[p0[Sqrt[x^2+z^2],ArcTan[z,x]]]^2]

ContourPlot[p2[x,z],{z,-10,10},{x,-10,10},PlotPoints->100,ContourLines->False]
Plot[p2[0,z],{z,-10,10}]
---
For[i=0,i<16,i++,dd[i]=d[i,.1]]
p0[r_,theta_]=Exp[I .1 r Cos[theta]]+
Sum[I^l (2l+1) (Exp[2 I dd[l]]-1) (sBesselj[l,.1 r] + I sBessely[l,.1 r])
LegendreP[l, Cos[theta]],{l,0,2}]/2

p2[x_,z_]=If[x^2+z^2<1,0,Abs[p0[Sqrt[x^2+z^2],ArcTan[z,x]]]^2]

ContourPlot[p2[x,z],{z,-100,100},{x,-100,100},PlotPoints->100,ContourLines->False]
---
tn[i_,n_,b_]= (Pi-2ArcSin[b])-i(Pi-2ArcSin[b/n]) 

Plot[{tn[2,1.5,b],tn[2,1.6,b],tn[2,1.7,b],tn[2,1.8,b],tn[2,1.9,b],tn[2,2,b]},
{b,0,1},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]},}]

Plot[{tn[1,1.5,b],tn[1,1.6,b],tn[1,1.7,b],tn[1,1.8,b],tn[1,1.9,b],tn[1,2,b]},
{b,0,1},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]},}]

FindRoot[tn[2,1.75,b]+Pi,{b,.85}]
b=0.847215
theta=ArcSin[b]
alpha=ArcSin[b/1.75]
phi=0
Graphics[{{RGBColor[1,0,0],Circle[{0,0},1]},Line[{{2.5,b},{Cos[theta],Sin[theta]},
{Cos[theta+Pi-2 alpha],Sin[theta+Pi-2 alpha]},
{Cos[theta+2Pi-4 alpha],Sin[theta+2Pi-4 alpha]},
{Cos[theta+2Pi-4 alpha]+1.5 Cos[phi],Sin[theta+2Pi-4 alpha]+1.5 Sin[phi]}}]} ]
Show[%,AspectRatio->Automatic]
-----
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
sBessely[l_,x_]=Sqrt[Pi/(2 x)] BesselY[l+(1/2),x]
sbesseljPrime[l_,k_]=D[sBesselj[l,k],k]
sbesselyPrime[l_,k_]=D[sBessely[l,k],k]
u0=50

d2[l_,k_]=ArcTan[
-(Sqrt[k^2+u0](sbesseljPrime[l,Sqrt[k^2+u0]])sBessely[l,k]
-sBesselj[l,Sqrt[k^2+u0]] k sbesselyPrime[l,k]),
-(Sqrt[k^2+u0](sbesseljPrime[l,Sqrt[k^2+u0]])sBesselj[l,k]
-sBesselj[l,Sqrt[k^2+u0]] k sbesseljPrime[l,k])]

k0 =k /. N[FindRoot[-(Sqrt[k^2+u0](sbesseljPrime[0,Sqrt[k^2+u0]])sBesselj[0,k]
-sBesselj[0,Sqrt[k^2+u0]] k sbesseljPrime[0,k]),{k,6.3},AccuracyGoal->13],12]
k1 =k /. N[FindRoot[-(Sqrt[k^2+u0](sbesseljPrime[1,Sqrt[k^2+u0]])sBesselj[1,k]
-sBesselj[1,Sqrt[k^2+u0]] k sbesseljPrime[1,k]),{k,6.8},AccuracyGoal->13],12]
k2 =k /. N[FindRoot[-(Sqrt[k^2+u0](sbesseljPrime[2,Sqrt[k^2+u0]])sBesselj[2,k]
-sBesselj[2,Sqrt[k^2+u0]] k sbesseljPrime[2,k]),{k,5.8},AccuracyGoal->13],12]
k31 =k /. N[FindRoot[-(Sqrt[k^2+u0](sbesseljPrime[3,Sqrt[k^2+u0]])sBesselj[3,k]
-sBesselj[3,Sqrt[k^2+u0]] k sbesseljPrime[3,k]),{k,5.4},AccuracyGoal->13],12]
k32 =k /. N[FindRoot[-(Sqrt[k^2+u0](sbesseljPrime[3,Sqrt[k^2+u0]])sBesselj[3,k]
-sBesselj[3,Sqrt[k^2+u0]] k sbesseljPrime[3,k]),{k,6.1},AccuracyGoal->13],12]

fd0[k_]=If[k<k0,d2[0,k]+2 Pi,d2[0,k]]
fd1[k_]=If[k<k1,d2[1,k]+2 Pi,d2[1,k]]
fd2[k_]=If[k<k2,d2[2,k]+2 Pi,d2[2,k]]
fd3[k_]=If[(k-k31)*(k-k32)<0,d2[3,k]+2 Pi,d2[3,k]]

Plot[{d2[0,k],d2[1,k],d2[2,k],d2[3,k],d2[4,k],d2[5,k]},{k,.1,30},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]}}]

Plot[{fd0[k],fd1[k],fd2[k],fd3[k],d2[4,k],d2[5,k]},{k,.1,30},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]}}]

Plot[{fd0[k],fd1[k],fd2[k],fd3[k],d2[4,k],d2[5,k]},{k,.1,10},PlotStyle->
{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]},{RGBColor[1,0,1]}}]


----
f[theta_,k_,lmax_]=1/(2 k I) (Sum[(Exp[2 I d2[l,k]]-1) (2l+1)
LegendreP[l, Cos[theta]],{l,0,lmax}])

Plot[{Abs[f[theta,.1,1]]^2,Abs[f[theta,.1,2]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,.3,1]]^2,Abs[f[theta,.3,2]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,1.,3]]^2,Abs[f[theta,1.,6]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,3.,5]]^2,Abs[f[theta,3.,10]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,10.,15]]^2,Abs[f[theta,10.,30]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}}]

Plot[{Abs[f[theta,.1,2]]^2,Abs[f[theta,.3,2]]^2,Abs[f[theta,1.,3]]^2,
Abs[f[theta,3.,6]]^2,Abs[f[theta,10.,15]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[1,0,0]},{RGBColor[1,1,0]},{RGBColor[0,1,0]},
{RGBColor[0,1,1]},{RGBColor[0,0,1]}},PlotRange->{0,1.5}]


Plot[{Abs[f[theta,30.,45]]^2,Abs[f[theta,30.,60]]^2},{theta,0,Pi},
PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]}}]
----
sigma1[k_]:=(4 Pi/k^2)(Sum[(2l+1)Sin[d2[l,k]]^2,{l,0,15}] )
Plot[sigma1[k],{k,.01,10},PlotPoints->100,PlotRange->{0,20},Ticks->{{0,2,4,6,8,10},{0,5,10,15,20}}]

sigma[k_]:=(4 Pi/k^2)(Sum[(2l+1)Sin[d2[l,k]]^2,{l,0,35}] )
Plot[sigma[k],{k,10,30},PlotRange->{0,20},Ticks->{{10,15,20,25,30},{0,5,10,15,20}}]

----
For[i=0,i<21,i++,dd[i]=d2[i,10.]]
For[i=0,i<21,i++,n[i]=.5((Exp[2 I dd[i]]+1)sBesselj[i,10.] + 
I (Exp[2 I dd[i]]-1)sBessely[i,10.] )/sBesselj[i,Sqrt[150.]] ]

p1[r_,theta_]:= Sum[ I^l (2l+1)n[l] (sBesselj[l,Sqrt[150.] r] )LegendreP[l, Cos[theta]],{l,0,20}]

p0[r_,theta_]=Exp[I 10. r Cos[theta]]+
Sum[I^l (l+.5)(Exp[2 I dd[l]]-1) (sBesselj[l,10. r] + I sBessely[l,10. r])LegendreP[l, Cos[theta]] ,{l,0,15}]

p2[x_,z_]=If[x^2+z^2<1,Abs[p1[Sqrt[x^2+z^2],ArcTan[z,x]]]^2,Abs[p0[Sqrt[x^2+z^2],ArcTan[z,x]]]^2]

ContourPlot[p2[x,z],{z,-3,3},{x,-3,3},PlotPoints->100,ContourLines->False]

Plot[p2[x,0],{x,-3,10}]
Plot[p2[x,.5],{x,-3,3}]
Plot[p2[0,z],{z,-3,10}]
----
For[i=0,i<21,i++,dd[i]=d2[i,1.]]
For[i=0,i<21,i++,n[i]=.5((Exp[2 I dd[i]]+1)sBesselj[i,1.] + 
I (Exp[2 I dd[i]]-1)sBessely[i,1.] )/sBesselj[i,Sqrt[51.]] ]

p1[r_,theta_]:= Sum[ I^l (2l+1)n[l] (sBesselj[l,Sqrt[51.] r] )LegendreP[l, Cos[theta]],{l,0,5}]

p0[r_,theta_]=Exp[I 1. r Cos[theta]]+
Sum[I^l (l+.5)(Exp[2 I dd[l]]-1) (sBesselj[l,1. r] + I sBessely[l,1. r])LegendreP[l, Cos[theta]] ,{l,0,5}]

p2[x_,z_]=If[x^2+z^2<1,Abs[p1[Sqrt[x^2+z^2],ArcTan[z,x]]]^2,Abs[p0[Sqrt[x^2+z^2],ArcTan[z,x]]]^2]

ContourPlot[p2[x,z],{z,-10,10},{x,-10,10},PlotPoints->100,ContourLines->False]

Plot[p2[x,0],{x,-3,10}]
Plot[p2[x,.5],{x,-3,3}]
Plot[p2[0,z],{z,-10,10}]
-
---
n=Sqrt[1.5]
tn[b_]= (-2ArcSin[b])+(2ArcSin[b/n]) 


b=0.95
theta1=ArcSin[b]
alpha1=ArcSin[b/n]
t1=tn[b]
g1=Graphics[{{RGBColor[1,0,0],Circle[{0,0},1]},Line[{{2.5,b},{Cos[theta1],Sin[theta1]},
{Cos[theta1+Pi-2 alpha1],Sin[theta1+Pi-2 alpha1]},
{Cos[theta1+Pi-2 alpha1]-2 Cos[t1],Sin[theta1+Pi-2 alpha1]+2 Sin[t1]}}], 
Line[{{2.5,-b},{Cos[theta1],-Sin[theta1]},
{Cos[theta1+Pi-2 alpha1],-Sin[theta1+Pi-2 alpha1]},
{Cos[theta1+Pi-2 alpha1]-2 Cos[t1],-Sin[theta1+Pi-2 alpha1]-2 Sin[t1]}}]} ]
b=0.75
theta2=ArcSin[b]
alpha2=ArcSin[b/n]
t2=tn[b]
g2=Graphics[{Line[{{2.5,b},{Cos[theta2],Sin[theta2]},
{Cos[theta2+Pi-2 alpha2],Sin[theta2+Pi-2 alpha2]},
{Cos[theta2+Pi-2 alpha2]-2 Cos[t2],Sin[theta2+Pi-2 alpha2]+2 Sin[t2]}}], 
Line[{{2.5,-b},{Cos[theta2],-Sin[theta2]},
{Cos[theta2+Pi-2 alpha2],-Sin[theta2+Pi-2 alpha2]},
{Cos[theta2+Pi-2 alpha2]-2 Cos[t2],-Sin[theta2+Pi-2 alpha2]-2 Sin[t2]}}]} ]
b=0.55
theta3=ArcSin[b]
alpha3=ArcSin[b/n]
t3=tn[b]
g3=Graphics[{Line[{{2.5,b},{Cos[theta3],Sin[theta3]},
{Cos[theta3+Pi-2 alpha3],Sin[theta3+Pi-2 alpha3]},
{Cos[theta3+Pi-2 alpha3]-2 Cos[t3],Sin[theta3+Pi-2 alpha3]+2 Sin[t3]}}], 
Line[{{2.5,-b},{Cos[theta3],-Sin[theta3]},
{Cos[theta3+Pi-2 alpha3],-Sin[theta3+Pi-2 alpha3]},
{Cos[theta3+Pi-2 alpha3]-2 Cos[t3],-Sin[theta3+Pi-2 alpha3]-2 Sin[t3]}}]} ]
b=0.35
theta4=ArcSin[b]
alpha4=ArcSin[b/n]
t4=tn[b]
g4=Graphics[{Line[{{2.5,b},{Cos[theta4],Sin[theta4]},
{Cos[theta4+Pi-2 alpha4],Sin[theta4+Pi-2 alpha4]},
{Cos[theta4+Pi-2 alpha4]-2 Cos[t4],Sin[theta4+Pi-2 alpha4]+2 Sin[t4]}}], 
Line[{{2.5,-b},{Cos[theta4],-Sin[theta4]},
{Cos[theta4+Pi-2 alpha4],-Sin[theta4+Pi-2 alpha4]},
{Cos[theta4+Pi-2 alpha4]-2 Cos[t4],-Sin[theta4+Pi-2 alpha4]-2 Sin[t4]}}]} ]
b=0.15
theta5=ArcSin[b]
alpha5=ArcSin[b/n]
t5=tn[b]
g5=Graphics[{Line[{{2.5,b},{Cos[theta5],Sin[theta5]},
{Cos[theta5+Pi-2 alpha5],Sin[theta5+Pi-2 alpha5]},
{Cos[theta5+Pi-2 alpha5]-2 Cos[t5],Sin[theta5+Pi-2 alpha5]+2 Sin[t5]}}], 
Line[{{2.5,-b},{Cos[theta5],-Sin[theta5]},
{Cos[theta5+Pi-2 alpha5],-Sin[theta5+Pi-2 alpha5]},
{Cos[theta5+Pi-2 alpha5]-2 Cos[t5],-Sin[theta5+Pi-2 alpha5]-2 Sin[t5]}}]} ]

Show[g1,g2,g3,g4,g5,AspectRatio->Automatic]

----
ClearAll[b]
n=Sqrt[51.]
tn[b_]= (-2ArcSin[b])+(2ArcSin[b/n]) 


b=0.95
theta1=ArcSin[b]
alpha1=ArcSin[b/n]
t1=tn[b]
g1=Graphics[{{RGBColor[1,0,0],Circle[{0,0},1]},Line[{{2.5,b},{Cos[theta1],Sin[theta1]},
{Cos[theta1+Pi-2 alpha1],Sin[theta1+Pi-2 alpha1]},
{Cos[theta1+Pi-2 alpha1]-2 Cos[t1],Sin[theta1+Pi-2 alpha1]+2 Sin[t1]}}], 
Line[{{2.5,-b},{Cos[theta1],-Sin[theta1]},
{Cos[theta1+Pi-2 alpha1],-Sin[theta1+Pi-2 alpha1]},
{Cos[theta1+Pi-2 alpha1]-2 Cos[t1],-Sin[theta1+Pi-2 alpha1]-2 Sin[t1]}}]} ]
b=0.75
theta2=ArcSin[b]
alpha2=ArcSin[b/n]
t2=tn[b]
g2=Graphics[{Line[{{2.5,b},{Cos[theta2],Sin[theta2]},
{Cos[theta2+Pi-2 alpha2],Sin[theta2+Pi-2 alpha2]},
{Cos[theta2+Pi-2 alpha2]-2 Cos[t2],Sin[theta2+Pi-2 alpha2]+2 Sin[t2]}}], 
Line[{{2.5,-b},{Cos[theta2],-Sin[theta2]},
{Cos[theta2+Pi-2 alpha2],-Sin[theta2+Pi-2 alpha2]},
{Cos[theta2+Pi-2 alpha2]-2 Cos[t2],-Sin[theta2+Pi-2 alpha2]-2 Sin[t2]}}]} ]
b=0.55
theta3=ArcSin[b]
alpha3=ArcSin[b/n]
t3=tn[b]
g3=Graphics[{Line[{{2.5,b},{Cos[theta3],Sin[theta3]},
{Cos[theta3+Pi-2 alpha3],Sin[theta3+Pi-2 alpha3]},
{Cos[theta3+Pi-2 alpha3]-2 Cos[t3],Sin[theta3+Pi-2 alpha3]+2 Sin[t3]}}], 
Line[{{2.5,-b},{Cos[theta3],-Sin[theta3]},
{Cos[theta3+Pi-2 alpha3],-Sin[theta3+Pi-2 alpha3]},
{Cos[theta3+Pi-2 alpha3]-2 Cos[t3],-Sin[theta3+Pi-2 alpha3]-2 Sin[t3]}}]} ]
b=0.35
theta4=ArcSin[b]
alpha4=ArcSin[b/n]
t4=tn[b]
g4=Graphics[{Line[{{2.5,b},{Cos[theta4],Sin[theta4]},
{Cos[theta4+Pi-2 alpha4],Sin[theta4+Pi-2 alpha4]},
{Cos[theta4+Pi-2 alpha4]-2 Cos[t4],Sin[theta4+Pi-2 alpha4]+2 Sin[t4]}}], 
Line[{{2.5,-b},{Cos[theta4],-Sin[theta4]},
{Cos[theta4+Pi-2 alpha4],-Sin[theta4+Pi-2 alpha4]},
{Cos[theta4+Pi-2 alpha4]-2 Cos[t4],-Sin[theta4+Pi-2 alpha4]-2 Sin[t4]}}]} ]
b=0.15
theta5=ArcSin[b]
alpha5=ArcSin[b/n]
t5=tn[b]
g5=Graphics[{Line[{{2.5,b},{Cos[theta5],Sin[theta5]},
{Cos[theta5+Pi-2 alpha5],Sin[theta5+Pi-2 alpha5]},
{Cos[theta5+Pi-2 alpha5]-2 Cos[t5],Sin[theta5+Pi-2 alpha5]+2 Sin[t5]}}], 
Line[{{2.5,-b},{Cos[theta5],-Sin[theta5]},
{Cos[theta5+Pi-2 alpha5],-Sin[theta5+Pi-2 alpha5]},
{Cos[theta5+Pi-2 alpha5]-2 Cos[t5],-Sin[theta5+Pi-2 alpha5]-2 Sin[t5]}}]} ]

Show[g1,g2,g3,g4,g5,AspectRatio->Automatic]
--
sBesselj[l_,x_]=Sqrt[Pi/(2 x)] BesselJ[l+(1/2),x]
sBessely[l_,x_]=Sqrt[Pi/(2 x)] BesselY[l+(1/2),x]
sbesseljPrime[l_,k_]=D[sBesselj[l,k],k]
sbesselyPrime[l_,k_]=D[sBessely[l,k],k]
u0=30^2(1.75^2-1)

d2[l_,k_]=ArcTan[
-(Sqrt[k^2+u0](sbesseljPrime[l,Sqrt[k^2+u0]])sBessely[l,k]
-sBesselj[l,Sqrt[k^2+u0]] k sbesselyPrime[l,k]),
-(Sqrt[k^2+u0](sbesseljPrime[l,Sqrt[k^2+u0]])sBesselj[l,k]
-sBesselj[l,Sqrt[k^2+u0]] k sbesseljPrime[l,k])]

For[i=0,i<51,i++,dd[i]=d2[i,30.]]

f[theta_]:=1/(2 30. I) (Sum[(Exp[2 I dd[l]]-1) (2l+1)LegendreP[l, Cos[theta]],{l,0,50}])

Plot[Abs[f[theta]]^2,{theta,0,Pi},PlotPoints->50]
---
<<NumericalMath`BesselZeros`
For[i=0,i<31,i++,z[i]=BesselJZeros[i,12-Floor[i/3]]; Print[Last[z[i]]^2]]
sqrtE=Sqrt[3000.]
n=0
For[j=1,j<12,j++,If[z[0][[j]]<sqrtE,n=n+1]]
For[i=1,i<31,i++,For[j=1,j<12-Floor[i/3],j++,If[z[i][[j]]<sqrtE,n=n+2] ] ] 
Print[n]
ListPlot[{{30,5},{100,21},{300,67},{1000,232}}]

For[i=0,i<55,i++,z[i]=BesselJZeros[i,18-Floor[i/3]]; Print[Last[z[i]]^2]]
sqrtE=Sqrt[3000.]
n=0
For[j=1,j<19,j++,If[z[0][[j]]<sqrtE,n=n+1]]
For[i=1,i<54,i++,For[j=1,j<18-Floor[i/3],j++,If[z[i][[j]]<sqrtE,n=n+2] ] ] 
Print[n]
ListPlot[{{30,5},{100,21},{300,67},{1000,232},{3000,725}}]