Variable phase calculation for potential scattering, vector case.

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.

jmδToIndex[j_,m_,δ_]:=3((j^2)+j +m)+Switch[δ,0,3,-1,1,1,2]
reducei[i_]:=Floor[(i/3)]
indexToj[i_]:=Floor[Sqrt[reducei[i]]]
indexTom[i_]:=reducei[i]-reducei[jmδToIndex[indexToj[i],0,1]]
indexToδ[i_]:=Switch[Mod[i,3],0,0,1,-1,2,1]
indexToℓ[i_]:=indexToj[i]+indexToδ[i]
minindex=1  (* drop ℓ=j=0 channel; other ℓ=0 channels will be projected out when dropping longitudinal component *)
maxj=1
maxindex=3(((maxj+1))^2)-1
maxr = 32
minr=(1/1024)
width=2
depth=2
potential[ℓ_,m_,r_]:=If[ℓ==0||(ℓ==1 && m≠0),(I^ℓ)depth Exp[((-(r^2))/(width^2))],0]
lmatrix=DiagonalMatrix[Table[indexToℓ[i](indexToℓ[i]+1),{i,minindex,maxindex}]]
signmatrix=DiagonalMatrix[Table[(((-1))^(indexToℓ[i])),{i,minindex,maxindex}]]
wigner[minindex,minindex,0,0]:=(1/(Sqrt[4π]))
wigner[i1_,i2_,ℓ_,m_]:=Module[{j1=indexToj[i1],m1=indexTom[i1],j2=indexToj[i2],m2=indexTom[i2],ℓ1=indexToℓ[i1],ℓ2=indexToℓ[i2]},
(((-1))^(j2+ℓ2+m2+1))(((-1))^m)(((-1))^(j2+ℓ1+ℓ))(Sqrt[(((2j1+1)(2j2+1)(2ℓ1+1)(2ℓ2+1))/(4 π(2ℓ+1)))])SixJSymbol[{ℓ1,ℓ2,ℓ},{j2,j1,1}]ClebschGordan[{ℓ1,0},{ℓ2,0},{ℓ,0}]ClebschGordan[{j1,m1},{j2,-m2},{ℓ,-m}]]

vmatrix[r_]:=Table[Module[{m=indexTom[ipp]-indexTom[ip]},Sum[potential[ℓ,m,r]wigner[ip,ipp,ℓ,m],{ℓ,Max[Abs[m],Abs[indexToℓ[ip]-indexToℓ[ipp]],Abs[indexToj[ip]-indexToj[ipp]]],Min[indexToℓ[ip]+indexToℓ[ipp],indexToj[ip]+indexToj[ipp]]}]],{ip,minindex,maxindex},{ipp,minindex,maxindex}]

radialFunction[ℓ_,x_]:=x SphericalHankelH1[Abs[ℓ],x]
logderivBesselMatrix[x_]=FullSimplify[FunctionExpand[DiagonalMatrix[Table[((D[radialFunction[indexToℓ[i],xx],xx]/.{xx->x})/(radialFunction[indexToℓ[i],x])),{i,minindex,maxindex}]]]]
besselMatrix[x_]=DiagonalMatrix[Table[radialFunction[indexToℓ[i],x],{i,minindex,maxindex}]]
besselMatrixInverse[x_]=DiagonalMatrix[Table[(1/(radialFunction[indexToℓ[i],x])),{i,minindex,maxindex}]]

solg[k_,x_:minr]:=First[NDSolve[{-f''[r]-2k f'[r].logderivBesselMatrix[k r]+(1/(r^2))(lmatrix.f[r]-f[r].lmatrix)+vmatrix[r].f[r]==0,
f[maxr]==IdentityMatrix[maxindex-minindex+1],f'[maxr]==DiagonalMatrix[Table[0,{i,minindex,maxindex}]]},f,{r,x,maxr}]]
solh[k_,x_:maxr]:=First[NDSolve[{-f''[r]+2k logderivBesselMatrix[k r].f'[r]-(1/(r^2))(lmatrix.f[r]-f[r].lmatrix)+f[r].vmatrix[r]+2(k^2)logderivBesselMatrix'[k r].f[r]==0,
f'[minr]==IdentityMatrix[maxindex-minindex+1],f[minr]==DiagonalMatrix[Table[0,{i,minindex,maxindex}]]},f,{r,minr,x}]]
wronk[g_,h_,x_,k_]:=besselMatrixInverse[k x].h[x].(g[x].(k besselMatrix'[k x])+g'[x].besselMatrix[k x])-(besselMatrixInverse[k x].h'[x]+(k besselMatrixInverse'[k x]).h[x]).g[x].besselMatrix[k x]

s=Module[{k=1.5,x=2.,g,h,gm,hm},g=f/.solg[k,x]; h=f/.solh[k,x];gm=f/.solg[-k,x]; hm=f/.solh[-k,x];Inverse[wronk[g,h,x,k]].besselMatrixInverse[k minr].besselMatrix[-k minr].wronk[gm,hm,x,-k].signmatrix]
Chop[s,(10^(-6))]//MatrixForm
Chop[s.ConjugateTranspose[s],(10^(-6))]//MatrixForm