Variable phase calculation for electromagnetic scattering.

All contents of this package are copyright Noah Graham and Aden Forrow, 2012, all rights reserved.  This package may be used and modified for noncommercial research purposes, provided that citation to

A. Forrow and N. Graham, "Variable Phase S-Matrix Calculations for Asymmetric Potentials and Dielectrics," arXiv:1210.0777, Phys. Rev. A 86 (2012) 062715

is included in all publications or other products in which the program or any programs derived from it were used.  This program and all accompanying materials are provided "AS IS," WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.

DOperator[A_,expr_]/;(FreeQ[A,#1]):=A expr
DOperator[L1_+L2_,expr_]:=DOperator[L1,expr]+DOperator[L2,expr]
DOperator[A_ L_,expr_]/;(FreeQ[A,#1]):=A DOperator[L,expr]
DOperator[A:HoldPattern[D[__]&],expr_]:=A[expr]
DOperator[L1_^n_Integer,expr_]/;n>1:=Nest[Expand[DOperator[L1,#1]]&,expr,n]
DOperator[A_,expr_]/;(!FreeQ[A,#1]):=A[expr]

Dr=(∂_r)#1 &;
indextoℓp[index_]:=Floor[Sqrt[index-(1/2)]]
indextomp[index_]:=-indextoℓp[index]+index-(((Floor[Sqrt[index-1]]))^2)-1
maxj:=2
maxIndexp:=(((maxj+1))^2)

indextoj[index_]:=Floor[Sqrt[Floor[(((index+4))/3)]-(1/2)]]
indextoℓ[index_]:=indextoj[index]+index-3Floor[(((index+1))/3)]
indextoδ[index_]:=indextoj[index]-indextoℓ[index]
indextom[index_]:=-indextoj[index]+Floor[(((index+4))/3)]-(((Floor[Sqrt[Floor[((index+4)/3)]-1]]))^2)-1
maxIndex:=3(((maxj+1))^2)-2
gradminus[l_]:=(Sqrt[(l/(2l+1))])(((l+1)/r)+Dr)
gradplus[l_]:=(Sqrt[((l+1)/(2l+1))])(-(l/r)+Dr)
 
gradOperator[L1_+L2_]:=gradOperator[L1]+gradOperator[L2]

(*This method will work as long as it never replaces "sphericalHarmonicY[...]" with the expression for the spherical harmonic*)
gradOperator[(L1_:1)*sphericalHarmonicY[l_,m_]]/;(FreeQ[L1,vectorSphericalHarmonic]&&FreeQ[L1,sphericalHarmonic]):=DOperator[gradminus[l],L1 ]*vectorSphericalHarmonicY[l,l-1,m]-DOperator[gradplus[l],L1 ]*vectorSphericalHarmonicY[l,l+1,m]

gradOperator[expr_]/;(FreeQ[expr,sphericalHarmonicY]&&FreeQ[expr,r]):=0
(*For VSH with l=j+1*)
divplus[j_]:=-(Sqrt[((j+1)/(2j+1))])(Dr+(((j+2))/r))
(*For VSH with l=j-1*)
divminus[j_]:=(Sqrt[(j/(2j+1))])(Dr-((j-1)/r))
divOperator[L1_+L2_]:=divOperator[L1]+divOperator[L2]
divOperator[A_(B_+C_)]:=divOperator[A B]+divOperator[A C]

divOperator[(L1_:1)*vectorSphericalHarmonicY[j_,l_,m_]]/;(FreeQ[L1,vectorSphericalHarmonic]&&FreeQ[L1,sphericalHarmonic]):=If[l===j+1,DOperator[divplus[j],L1]*sphericalHarmonicY[j,m],If[l===j-1,DOperator[divminus[j],L1]*sphericalHarmonicY[j,m],0]]
divOperator[L1_]/;(FreeQ[L1,r]&&FreeQ[L1,vectorSphericalHarmonicY]):=0


(*For curls of VSH with l=j+1*)
pluscurl[j_]:=I (Sqrt[(j/(2j+1))])(Dr + ((j+2)/r))

(*For curls of VSH with l=j*)
zerocurlplus[j_]:=I (Sqrt[(j/(2j+1))])(Dr-(j/r))
zerocurlminus[j_]:=I(Sqrt[((j+1)/(2j+1))])(Dr+(((j+1))/r))

(*For curls of VSH with l=j-1*)
minuscurl[j_]:=I(Sqrt[((j+1)/(2j+1))])(Dr-(((j-1))/r))
curlOperator[L1_+L2_]:=curlOperator[L1]+curlOperator[L2]

curlOperator[(L1_:1)*vectorSphericalHarmonicY[j_,l_,m_]]/;(FreeQ[L1,vectorSphericalHarmonic]&&FreeQ[L1,sphericalHarmonic]):=If[l===j+1,DOperator[pluscurl[j],L1]*vectorSphericalHarmonicY[j,j,m],If[l===j-1,DOperator[minuscurl[j],L1]*vectorSphericalHarmonicY[j,j,m],If[l===j,DOperator[zerocurlplus[j],L1]*vectorSphericalHarmonicY[j,j+1,m]+DOperator[zerocurlminus[j],L1]*vectorSphericalHarmonicY[j,j-1,m],0]]]

(*This piece is necessary mostly so that curlOperator[0]=0. Any other argument should be some constant or function of r times a vector spherical harmonic.*)
curlOperator[L1_]/;(FreeQ[L1,r]&&FreeQ[L1,vectorSphericalHarmonicY]):=0

(*This operator is equivalent to multiplying (dot product) by the complex conjugate of VSH[j,l,m] and integrating over solid angle, if the function has been expanded in a basis of VSH*)

angleIntegrate[L1_+L2_,j_,l_,m_]:=angleIntegrate[L1,j,l,m]+angleIntegrate[L2,j,l,m]

(*If the angular dependent part is not a VSH, angleIntegrate gives zero. This is not a problem if everything is expanded as VSH*)
angleIntegrate[A_*B_,j_,l_,m_]/;FreeQ[A,vectorSphericalHarmonicY]:=A angleIntegrate[B,j,l,m]
angleIntegrate[vectorSphericalHarmonicY[j_,l_,m_],j_,l_,m_]:=1

angleIntegrate[A_,j_,l_,m_]/;FreeQ[A,vectorSphericalHarmonicY[j,l,m]]:=0

Off[ClebschGordan::tri]
Off[ClebschGordan::phy]
Off[SixJSymbol::tri]

dotproduct[vectorSphericalHarmonicY[j1_,l1_,m1_],vectorSphericalHarmonicY[j2_,l2_,m2_]]:=Sum[Simplify[(((-1))^(j2 +l1+indextoℓp[i]))(Sqrt[(((2j1+1)(2j2+1)(2l1+1)(2l2+1))/(4π(2indextoℓp[i] +1)))])SixJSymbol[{l1,l2,indextoℓp[i]},{j2,j1,1}]ClebschGordan[{l1,0},{l2,0},{indextoℓp[i],0}]ClebschGordan[{j1,m1},{j2,m2},{indextoℓp[i],indextomp[i]}]]*sphericalHarmonicY[indextoℓp[i],indextomp[i]],{i,1,maxIndexp}]

dotproduct[A1_,A_ A2_]/;FreeQ[A,vectorSphericalHarmonicY]:=A dotproduct[A1,A2]
dotproduct[A_ A1_,A2_]/;FreeQ[A,vectorSphericalHarmonicY]:=A dotproduct[A1,A2]
dotproduct[A1_,A2_+A3_]:=dotproduct[A1,A2]+dotproduct[A1,A3]
dotproduct[A1_+A2_,A3_]:=dotproduct[A1,A3]+dotproduct[A2,A3]

(*We need to define a new "product" function to multiply scalar and vector spherical harmonics together.*)

(*Scalar harmonics are put in first*)
product[sphericalHarmonicY[l_,m_],vectorSphericalHarmonicY[j2_,l2_,m2_]]:=Sum[(((-1))^(indextoj[i]+indextoℓ[i]+indextom[i]+1))(((-1))^(j2+indextoℓ[i]+l)) (((-1))^m) (Sqrt[(((2 indextoj[i]+1) (2 j2+1) (2 indextoℓ[i]+1) (2 l2+1))/(4 π (2 l+1)))])SixJSymbol[{indextoℓ[i],l2,l},{j2,indextoj[i],1}] ClebschGordan[{indextoℓ[i],0},{l2,0},{l,0}] ClebschGordan[{indextoj[i],-indextom[i]},{j2,m2},{l,-m}]*vectorSphericalHarmonicY[indextoj[i],indextoℓ[i],indextom[i]],{i,1,maxIndex}]

product[L1_+L2_,L3_]:=product[L1,L3]+product[L2,L3]
product[L1_,L2_+L3_]:=product[L1,L2]+product[L1,L3]
product[A_ L1_,L2_]/;(FreeQ[A,sphericalHarmonicY]):=A product[L1,L2]
product[L1_,A_ L2_]/;(FreeQ[A,vectorSphericalHarmonicY]):=A product[L1,L2]

evector[r_]:=Sum[(r^(-1))(Subscript[e, {indextoj[i],indextoℓ[i],indextom[i]}]) [r]*vectorSphericalHarmonicY[indextoj[i],indextoℓ[i],indextom[i]],{i,1,maxIndex}]
myδ[a_,b_,c_,ip_]:=KroneckerDelta[a,indextoj[ip]]KroneckerDelta[b,indextoℓ[ip]]KroneckerDelta[c,indextom[ip]]

depth=4/(Sqrt[4 π])
steepness=4096
width=1
depth2=(1/4)
ϵ[ℓ_,m_]/;(m≠0 ||ℓ>1):=Function[0]
ϵtotal[r_]:=Sum[ϵ[indextoℓp[i],indextomp[i]][r]*sphericalHarmonicY[indextoℓp[i],indextomp[i]],{i,1,maxIndexp}]
$RecursionLimit=∞
$Assumptions={r>0}
Module[{curlcurl,ϵevec,ϵdivϵ,eqn},
curlcurl[k_]=curlOperator[curlOperator[evector[r]]];
ϵevec[k_]=product[ϵtotal[r],evector[r]];
ϵdivϵ[k_]=product[ϵtotal[r],gradOperator[divOperator[ϵevec[k]]]];
eqn[k_]=Expand[curlcurl[k]-ϵdivϵ[k]-(k^2)ϵevec[k]];
system[k_]=Simplify[Table[angleIntegrate[eqn[k],indextoj[i],indextoℓ[i],indextom[i]],{i,1,maxIndex}]];
]
d2[k_]=Table[-(r^3)system[k][[i]]/.{(Subscript[e, {a_,b_,c_}])[r]->0,(((Subscript[e, {a_,b_,c_}])))'[r]->0,(((Subscript[e, {a_,b_,c_}])))''[r]->myδ[a,b,c,ip]},{i,1,maxIndex},{ip,1,maxIndex}];
d1[k_]=Table[(r^3)system[k][[i]]/.{(Subscript[e, {a_,b_,c_}])[r]->0,(((Subscript[e, {a_,b_,c_}])))''[r]->0,(((Subscript[e, {a_,b_,c_}])))'[r]->myδ[a,b,c,ip]},{i,1,maxIndex},{ip,1,maxIndex}];
d[k_]=Table[(r^3)system[k][[i]]/.{(((Subscript[e, {a_,b_,c_}])))'[r]->0,(((Subscript[e, {a_,b_,c_}])))''[r]->0,(Subscript[e, {a_,b_,c_}])[r]->myδ[a,b,c,ip]},{i,1,maxIndex},{ip,1,maxIndex}];
d2inv[k_]=Simplify[Inverse[d2[k]]];

d12[k_]=Simplify[d2inv[k].d1[k]];
d02[k_]=Simplify[d2inv[k].d[k]+(k^2)IdentityMatrix[maxIndex]];
d12prime[k_]=Simplify[D[d12[k],r]];

ϵ[0,0]:=Function[r,(Sqrt[4 π])(1-(depth/2) (Tanh[steepness (r-width)]-1))]
ϵ[1,0]:=Function[r,depth2 Exp[((-(r^2))/(width^2))]]

radialFunction[ℓ_,x_]:=x SphericalHankelH1[Abs[ℓ],x]
logderivRadialFunction[ℓ_,x_]:=Simplify[FunctionExpand[((D[radialFunction[ℓ,xx],xx]/.{xx->x})/(radialFunction[ℓ,x]))]]
signmatrix=DiagonalMatrix[Table[(((-1))^(indextoℓ[i])),{i,1,maxIndex}]];
besselMatrix[x_]=DiagonalMatrix[Table[radialFunction[indextoℓ[i],x],{i,1,maxIndex}]];
besselMatrixInverse[x_]=DiagonalMatrix[Table[(1/(radialFunction[indextoℓ[i],x])),{i,1,maxIndex}]];
logDerivBesselMatrix[x_]=DiagonalMatrix[Table[logderivRadialFunction[indextoℓ[i],x],{i,1,maxIndex}]];
logDerivBesselMatrixPrime[x_]=Simplify[logDerivBesselMatrix'[x]];

lmatrix=DiagonalMatrix[Table[indextoℓ[i](indextoℓ[i]+1),{i,1,maxIndex}]];
Nvector[j_,m_]:={Table[KroneckerDelta[indextoj[i],j]KroneckerDelta[indextom[i],m](KroneckerDelta[indextoℓ[i],j-1](Sqrt[((j+1)/(2j+1))])-KroneckerDelta[indextoℓ[i],j+1](Sqrt[(j/(2j+1))])),{i,1,maxIndex}]}
projectN[j_,m_]:=Transpose[Nvector[j,m]].Nvector[j,m]
Lvector[j_,m_]:={Table[KroneckerDelta[indextoj[i],j]KroneckerDelta[indextom[i],m](KroneckerDelta[indextoℓ[i],j+1](Sqrt[((j+1)/(2j+1))])+KroneckerDelta[indextoℓ[i],j-1](Sqrt[(j/(2j+1))])),{i,1,maxIndex}]}
projectL[j_,m_]:=Transpose[Lvector[j,m]].Lvector[j,m]
transverseNmatrix=Sum[projectN[indextoℓp[i],indextomp[i]],{i,1,maxIndexp}];
transverseMmatrix=Table[KroneckerDelta[i,ip]KroneckerDelta[indextoj[i],indextoℓ[i]],{i,1,maxIndex},{ip,1,maxIndex}];
transverseprojection=transverseNmatrix+transverseMmatrix;
infinity=16
zero=(1/128)
matrixsystemg[k_]=-gmatrix''[r]+(d12[k].gmatrix[r]-2 gmatrix'[r]).logDerivBesselMatrix[k r]k+d12[k].gmatrix'[r]+d02[k].gmatrix[r]-gmatrix[r].(lmatrix/(r^2))==0;
boundary1g=gmatrix[infinity]==IdentityMatrix[maxIndex];
boundary2g=gmatrix'[infinity]==DiagonalMatrix[Table[0,{i,1,maxIndex}]];

matrixsystemh[k_]=-hmatrix''[r]+k logDerivBesselMatrix[k r].(2 hmatrix'[r]+hmatrix[r].d12[k])+2(k^2)logDerivBesselMatrixPrime[k r].hmatrix[r]-hmatrix'[r].d12[k]-hmatrix[r].d12prime[k]+hmatrix[r].d02[k]-(lmatrix/(r^2)).hmatrix[r]==0;
boundary1h=hmatrix'[zero]==IdentityMatrix[maxIndex];
boundary2h=hmatrix[zero]==DiagonalMatrix[Table[0,{i,1,maxIndex}]];

solg[k_,x_:zero]:=First[NDSolve[{matrixsystemg[k],boundary1g,boundary2g},gmatrix,{r,x,infinity},MaxSteps->100000]]
solh[k_,x_:infinity]:=First[NDSolve[{matrixsystemh[k],boundary1h,boundary2h},hmatrix,{r,zero,x},MaxSteps->100000]]

wronk[g_,h_,x_,k_]:=besselMatrixInverse[k x].(h[x].(k g[x].logDerivBesselMatrix[k x]+g'[x])-(h'[x]-k logDerivBesselMatrix[k x].h[x]).g[x] - h[x].(d12[k]/.{r->x}).g[x]).besselMatrix[k x]

s=Module[{k=1.5,x=0.5,g,h,gm,hm},g=gmatrix/.solg[k,x]; h=hmatrix/.solh[k,x];gm=gmatrix/.solg[-k,x]; hm=hmatrix/.solh[-k,x];Inverse[wronk[g,h,x,k]].besselMatrixInverse[k zero].besselMatrix[-k zero].wronk[gm,hm,x,-k].signmatrix]
s.Conjugate[Transpose[s]]//MatrixForm
Eigenvalues[s]//MatrixForm
Eigenvalues[transverseprojection.s.transverseprojection]//MatrixForm