Cell["\<\
(*      SphereGreensFunctions.m

\tAuthor S.Arridge
\tDate: September 1992

\tlast update : Feb 14th 1996

\tThis file contains various Mathematica functions for computing
\tIntegrated intensity and mean time on boundary of a sphere

\tNote : extended to include extrapolated boundary condition
        July 15th 2002!:
             extension to Robin boundary conditions

*)
(*=========================================================================*)
BeginPackage[\"SphereGreensFunctions`\"]

SphereGreensFunctions::usage =
\"SphereGreensFunctions supplies various functions for calculating Green's
functions of the diffusion equation on an homogeneous sphere\\n
\\n
Functions supplied (use ? for help) :\\n
\\tSphereIntegIntensityAndMT\\n
\\tSphereComplexData\\n
\\tExtrapolatedSphereIntegIntensityAndMT\\n
\\tExtrapolatedSphereComplexData\\n
\\tNeumannSphereIntegIntensityAndMT\\n
\\tListAngles\\n
\\tListIntegIntensity\\n
\\tListMeanTime\\n
\\tListComplexIntensity\\n
\\tListPhase\\n
\\tListModulation\\n
\\tDisplayIntegIntensity\\n
\\tDisplayMeanTime\\n
\\tDisplayPhase\\n
\\tDisplayModulation\\n
\\tComparePhaseAndMT\\n
\"
SphereIntegIntensityAndMT::usage =
\"SphereIntegIntensityAndMT[filename, a, s, r, ref, nang, P] calculates \
integrated 
intensity and mean time on the surface of a sphere.\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n
\\tReturned Value : {angle (list),intensity (list), \\n
\\t\\t\\ttime weighted intensity (list), mean time (list)}\\n
\\n
\\tExample:\\n
  soln=SphereIntegIntensityAndMT[\\\"sphere.results\\\",\\n
\\t\\t\\t0.025, 2, 25, 1.4, 16, 50];
\"

SphereIntegIntensityAndMT/: \
SphereIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\tspherecalc1[out,NTerm, nang, P]
    ]

SphereComplexData::usage =
\"SphereComplexData[filename, a, s, r, ref, nang, mhz,P] calculates complex \
intensity on the surface of a sphere.\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tmhz is frequency in mhz\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n
\\n
\\tReturned Value : {angle (list),complex_intensity (list),frequency \
(number)}
\\n
\\tExample:\\n
  soln=SphereComplexData[\\\"sphere.results\\\",\\n
\t0.025, 2, 25, 1.4, 16, 200, 50];
\"

SphereComplexData/: SphereComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\tspherecalc[out,NTerm, mhz, nang, P]
    ]

ExtrapolatedSphereIntegIntensityAndMT::usage =
\"ExtrapolatedSphereIntegIntensityAndMT[filename, a, s, r, ref, nang, P] \
calculates \\n
integrated intensity and mean time on the surface of a sphere, using\\n
extrapolated boundary condition.\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n
\\n
\\tReturned Value : {angle (list),intensity (list), \\n
\\t\\t\\ttime weighted intensity (list), mean time (list)}\\n
\\n
\\tExample:\\n
  soln=ExtrapolatedSphereIntegIntensityAndMT[\\\"sphere.results\\\",\\n
\\t\\t\\t0.025, 2, 25, 1.4, 16, 50];
\"

ExtrapolatedSphereIntegIntensityAndMT/: \
ExtrapolatedSphereIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\textrapolatedspherecalc1[out,NTerm, nang, P]
    ]

ExtrapolatedSphereComplexData::usage =
\"ExtrapolatedSphereComplexData[filename, a, s, r, ref, nang, mhz,P] \
calculates complex intensity on the surface of a sphere.\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tmhz is frequency in mhz\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n

\\tReturned Value : {angle (list),complex_intensity (list),frequency \
(number)}
\\n
\\tExample:\\n
  soln = ExtrapolatedSphereComplexData[\\\"sphere.results\\\",\\n
\\t\\t\\t0.025, 2, 25, 1.4, 16, 200, 50];
\"

ExtrapolatedSphereComplexData/: \
ExtrapolatedSphereComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\textrapolatedspherecalc[out,NTerm, mhz, nang, P]
    ]



RobinSphereComplexData::usage =
\"RobinSphereComplexData[filename, a, s, r, ref, nang, mhz,P] calculates \
complex intensity on the surface of a sphere.\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tmhz is frequency in mhz\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n

\\tReturned Value : {angle (list),complex_intensity (list),frequency \
(number)}
\\n
\\tExample:\\n
  soln = RobinSphereComplexData[\\\"sphere.results\\\",\\n
\\t\\t\\t0.025, 2, 25, 1.4, 16, 200, 50];
\"

RobinSphereComplexData/: \
RobinSphereComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\trobinspherecalc[out,NTerm, mhz, nang, P]
    ]

NeumannSphereIntegIntensityAndMT::usage =
\"NeumannSphereIntegIntensityAndMT[filename, a, s, r, ref, nang, P] \
calculates integrated 
intensity and mean time on the surface of a sphere, usung a Neumann source\\n
\\ta is value of mua in mm^-1\\n
\\ts is value of reduced mus (i.e. (1-g)mus) in mm^-1\\n
\\tr is radius of sphere in mm\\n
\\tref is refractive index\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tP is working precision. <I suggest 50>\\n
\\tfilename is an output file - \\n
\\t\\tNOTE : filename must be supplied with quotes, \\n
\\t\\tso that Mathematica interprets it as a string.\\n
\\n
\\tReturned Value : {angle (list),intensity (list), \\n
\\t\\t\\ttime weighted intensity (list), mean time (list)}\\n
\\n
\\n
\\tExample:\\n
  soln=NeumannSphereIntegIntensityAndMT[\\\"sphere.results\\\",\\n
\\t\\t\\t0.025, 2, 25, 1.4, 16, 50];
\"

NeumannSphereIntegIntensityAndMT/: \
NeumannSphereIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
\tSetParameters[a,s,r,ref,P];
\tNTerm = Floor[(100*s*a*r*r/nang)]*nang;
\tIf[NTerm < 1536, NTerm = 1536];
\tNeumannSphereCalc1[out,NTerm, nang, P]
    ]

ListAngles::usage =
\"ListAngles prints the angles (in radians) at which functions calculated\"

ListAngles/: ListAngles := SphereGreensFunctions`Private`angtab

ListIntegIntensity::usage =
\"ListIntegIntensity prints the integrated intensity at each angle\"

ListIntegIntensity := SphereGreensFunctions`Private`if5

ListMeanTime::usage =
\"ListMeanTime prints the mean time at each angle\"

ListMeanTime := SphereGreensFunctions`Private`nf5

ListComplexIntensity::usage =
\"ListComplexIntensity prints the complex intensity at each point\"

ListComplexIntensity := SphereGreensFunctions`Private`pf5

ListPhase::usage =
\"ListPhase prints the phase at each point\"

ListPhase := SphereGreensFunctions`Private`af5

ListModulation::usage =
\"ListPhase prints the phase at each point\"

ListModulation := Abs[SphereGreensFunctions`Private`pf5]

DisplayIntegIntensity::usage =
\"DisplayIntegIntensity plots the log of integrated intensity at each angle\"

DisplayIntegIntensity :=
\tModule[{pairs},
\tpairs = Transpose[ {N[180/Pi ListAngles],Log[10,ListIntegIntensity]} ];
\tListPlot[pairs, PlotJoined->True, PlotLabel->\"Log Intensity on sphere\",
\t\tAxesLabel->{\"Degrees\",\"Log_10 I\"}]
\t]

DisplayMeanTime::usage =
\"DisplayMeanTime plots the mean time at each angle\"

DisplayMeanTime :=
\tModule[{pairs},
\tpairs = Transpose[ {N[180/Pi ListAngles],ListMeanTime} ];
\tListPlot[pairs, PlotJoined->True,PlotLabel->\"Mean Time on sphere\",
\t\tAxesLabel->{\"Degrees\",\"<t> picoseconds\"}]
\t]

DisplayPhase::usage =
\"DisplayPhase plots the phase at each angle\"

DisplayPhase :=
\tModule[{pairs},
\tpairs = Transpose[ {N[180/Pi ListAngles],ListPhase} ];
\tListPlot[pairs, PlotJoined->True,PlotLabel->\"Phase on sphere\",
\t\tAxesLabel->{\"Degrees\",\"Arg[G] radians\"}]
\t]

DisplayModulation::usage =
\"DisplayModulation plots the amplitude of the complex intensity at each \
angle\"

DisplayModulation :=
\tModule[{pairs},
\tpairs = Transpose[ {N[180/Pi ListAngles],
\t\t\tLog[10,Abs[ListComplexIntensity]]} ];
\tListPlot[pairs, PlotJoined->True,PlotLabel->\"Modulation on sphere\",
\t\tAxesLabel->{\"Degrees\",\"Abs[G] \"}]
\t]

ComparePhaseAndMT::usage =
\"ComparePhaseAndMT[complex_solution,mean_time_solution] \\n
Displays a joint graph of the phase and mean time*frequency. Input are the
expexted solution lists from a function of time <name>ComplexIntensity and
<name>IntegIntensityAndMeanTime respectively\"

ComparePhaseAndMT[soln1_,soln2_] :=
\tModule[{pairs1,pairs2,w,g1,g2},
\tw = soln1[[3]];         (* frequency *)
\tpairs1 = Transpose[ {soln1[[1]],-Arg[soln1[[2]]]/w} ];
\tg1=ListPlot[pairs1, PlotJoined->True,PlotStyle->{RGBColor[1,0,0]},
\t\tDisplayFunction->Identity];
\tpairs2 = Transpose[ {soln2[[1]],soln2[[4]]} ];
\tg2=ListPlot[pairs2, PlotJoined->True,PlotStyle->{RGBColor[0,1,0]},
\t\tDisplayFunction->Identity];
\tShow[g1,g2,PlotLabel->\"Phase and MeanTime on sphere\",
\t\tAxesLabel->{\"Degrees\",\"Phase,<t>\"},
\t\tDisplayFunction:>$DisplayFunction]
\t]
(*=========================================================================*)
Begin[\"`Private`\"]

SetParameters::usage = 
\"SetParameters[a, s, r, ref, P] defines mua to be a, reduced mus to be s
a sphere of radius r to be used, and a refractive index of ref. All at a
precision of P\";

SetParameters/: SetParameters[a_, s_, r_, ref_, P_] :=
    Module[{},
\tmua  = SetPrecision[a,P];
\tmusd = SetPrecision[s,P];
\tsa = SetPrecision[r,P];
\trefind = SetPrecision[ref,P];
\tzeta = sa - 1/musd;
\tic = 3/10;
\tc = ic/refind;
\tr0 = (refind-1)^2 / (refind+1)^2;
\tthetac = ArcSin[1/refind];
\treffac = (2/(1-r0) - 1 + Abs[Cos[thetac]]^3) / (1 - Abs[Cos[thetac]]^2);
\tsab = sa + ebP[P];
\trA  = 2; (* Robin boundary condition, set to 2 for now *)
\tPrint[\"Setting mua : \",a,\" mus \",s,\" radius \",r,\" rind \",ref];
\ta0 = a0P[P];
    ] 
(*=========================================================================*)
(*
mua = 25/100;
musd = 2;
sa = 25;
zeta = sa - 1/musd;
ic = 3/10;
refind = 14/10;
c = ic/refind;
r0 = (refind-1)^2 / (refind+1)^2;
thetac = ArcSin[1/refind];
reffac = (2/(1-r0) - 1 + Abs[Cos[thetac]]^3) / (1 - Abs[Cos[thetac]]^2);
a0 = a0P[40];
*)
gamma/: gamma[P_] := N[Sqrt[c/(3*(mua + musd))],P]
gammaSym/: gammaSym := Sqrt[c/(3*(mua + musd))]
ebP/: ebP[P_] := N[Sqrt[1/(3*mua*(mua + musd))],P] * 
\tArcTanh[ 2 N[reffac Sqrt[ mua / (3*(mua + musd)) ],P]]
alpha/: alpha[omega_,P_] := N[Sqrt[mua*c + I*omega],P]/gamma[P]
alpha/: alpha[omega_] := Sqrt[mua*c + I*omega]/gammaSym
a0P/: a0P[P_] := alpha[0,P]
(*=========================================================================*)
DBessI/:DBessI[n_,x_]:=D[BesselI[n,xx],xx]/.{xx->x}

DBessK/:DBessK[n_,x_]:=D[BesselK[n,xx],xx]/.{xx->x}

sterm/: sterm[n_, omega_, theta_, P_] :=
 (BesselI[n + 1/2,alpha[omega,P]*zeta]*(2*n + \
1)*LegendreP[n,N[Cos[theta],P]])/
  BesselI[n + 1/2, alpha[omega,P]*sa]


btermI/: btermI[n_, omega_, P_] :=
   (BesselI[n + 1/2, alpha[omega,P]*zeta]*(2*n + 1))/
    BesselI[n + 1/2, alpha[omega,P]*sa]

fastbtermI/: fastbtermI[n_] :=
   (BItabzeta[[n + 1]]*(2*n + 1))/BItabsa[[n+1]]

fastbtermIR/: fastbtermIR[n_,rbc_] :=
   (BItabzeta[[n + 1]]*(2*n + 1))/(BItabsa[[n+1]]+rbc*BDItabsa[[n+1]])

dbtermI/: dbtermI[n_,P_] :=  (*  careful, what is a0 ? *)
   -(2 n + 1) (zeta*mdbi[n + 1/2, a0*zeta]*BesselI[n + 1/2, a0*sa] -
                sa*mdbi[n + 1/2, a0*sa]*BesselI[n + 1/2, a0*zeta])/
        BesselI[n + 1/2, a0*sa]^2

fastdbtermI/: fastdbtermI[n_, P_ ] :=
   If[ n > 0, 
   -(2 n + 1)( zeta*(BItabzeta[[n]]+BItabzeta[[n + 2]])*BItabsa[[n+1]] -
                sa*(BItabsa[[n]]+BItabsa[[n + 2]])*BItabzeta[[n + 1]] )/
        (2 BItabsa[[n+1]]^2),
    dbtermI[n,P] 
    ]

btermK/: btermK[n_, omega_, P_] := 
    (2 n +1)*BesselI[n+1/2, alpha[omega,P]*zeta]*(alpha[omega,P]*
     (D[BesselK[n+1/2, y],y] - 
\tD[BesselI[n+1/2, y],y] * BesselK[n+1/2, alpha[omega,P]*sab]/
\t\tBesselI[n+1/2, alpha[omega,P]*sab])  -
    1/(2 sa)* (BesselK[n+1/2, y] - 
\tBesselI[n+1/2, y] * BesselK[n+1/2, alpha[omega,P]*sab]/
\t\tBesselI[n+1/2, alpha[omega,P]*sab]))/.{y->alpha[omega,P]*sa} 

fastbtermK/: fastbtermK[n_,omega_,P_] :=
   If[ n > 0, 
   (2 n +1)*BItabzeta[[n + 1]]*(alpha[omega,P]*
\t( -(BKtabsa[[n]] + BKtabsa[[n+2]])/2 -   
\t(BItabsa[[n]]+BItabsa[[n+2]])*BKtabsab[[n+1]]/(2*BItabsab[[n+1]]) )  -
   1/(2 sa)*
\t(BKtabsa[[n+1]] - BItabsa[[n+1]]*BKtabsab[[n+1]]/BItabsab[[n+1]] ) ) ,
    btermK[n,omega,P] 
   ]

(* This function got by differentiating btermK, w.r.t alpha.

\tNot correct for extrapolated case...
*)
dbtermK/: dbtermK[n_,P_] :=
((2 n +1)*I*D[alpha[omega,P],omega]/.omega->0) * (BesselI[n+1/2, a0*zeta]*
     (-((BesselI[-1 + n+1/2, a0*sa] + BesselI[1 + n+1/2, a0*sa])*
\t   BesselK[n+1/2, a0*sab])/(2*BesselI[n+1/2, a0*sab]) + 
       (-BesselK[-1 + n+1/2, a0*sa] - BesselK[1 + n+1/2, a0*sa])/2) + 
    (a0*zeta*(BesselI[-1 + n+1/2, a0*zeta] + BesselI[1 + n+1/2, a0*zeta])*
       (-((BesselI[-1 + n+1/2, a0*sa] + BesselI[1 + n+1/2, a0*sa])*
\t     BesselK[n+1/2, a0*sab])/(2*BesselI[n+1/2, a0*sab]) + 
\t (-BesselK[-1 + n+1/2, a0*sa] - BesselK[1 + n+1/2, a0*sa])/2))/2 + 
    a0*BesselI[n+1/2, a0*zeta]*
     ((sab*(BesselI[-1 + n+1/2, a0*sa] + BesselI[1 + n+1/2, a0*sa])*
\t  (BesselI[-1 + n+1/2, a0*sab] + BesselI[1 + n+1/2, a0*sab])*
\t  BesselK[n+1/2, a0*sab])/(4*BesselI[n+1/2, a0*sab]^2) - 
       (((sa*(BesselI[-2 + n+1/2, a0*sa] + BesselI[n+1/2, a0*sa]))/2 + 
\t    (sa*(BesselI[n+1/2, a0*sa] + BesselI[2 + n+1/2, a0*sa]))/2)*
\t  BesselK[n+1/2, a0*sab])/(2*BesselI[n+1/2, a0*sab]) - 
       (sab*(BesselI[-1 + n+1/2, a0*sa] + BesselI[1 + n+1/2, a0*sa])*
\t  (-BesselK[-1 + n+1/2, a0*sab] - BesselK[1 + n+1/2, a0*sab]))/
\t(4*BesselI[n+1/2, a0*sab]) + 
       (-(sa*(-BesselK[-2 + n+1/2, a0*sa] - BesselK[n+1/2, a0*sa]))/2 - 
\t  (sa*(-BesselK[n+1/2, a0*sa] - BesselK[2 + n+1/2, a0*sa]))/2)/2) )

fastdbtermK/: fastdbtermK[n_,P_] :=
(*
\tjust need to substitute in the tabulated Bessel functions
\tNot correct for extrapolated case...
*)
If[ n > 1, 
 ((2 n +1)I*D[alpha[omega,P],omega]/.omega->0) * (BItabzeta[[1+n]]*
     (-((BItabsa[[n]] + BItabsa[[2 + n]])*
\t   BKtabsab[[1+n]])/(2*BItabsab[[1+n]]) + 
       (-BKtabsa[[n]] - BKtabsa[[2 + n]])/2) + 
    (a0*zeta*(BItabzeta[[n]] + BItabzeta[[2 + n]])*
       (-((BItabsa[[n]] + BItabsa[[2 + n]])*
\t     BKtabsab[[1+n]])/(2*BItabsab[[1+n]]) + 
\t (-BKtabsa[[n]] - BKtabsa[[2 + n]])/2))/2 + 
    a0*BItabzeta[[1+n]]*
     ((sab*(BItabsa[[n]] + BItabsa[[2 + n]])*
\t  (BItabsab[[n]] + BItabsab[[2 + n]])*
\t  BKtabsab[[1+n]])/(4*BItabsab[[1+n]]^2) - 
       (((sa*(BItabsa[[-1 + n]] + BItabsa[[1+n]]))/2 + 
\t    (sa*(BItabsa[[1+n]] + BItabsa[[3 + n]]))/2)*
\t  BKtabsab[[1+n]])/(2*BItabsab[[1+n]]) - 
       (sab*(BItabsa[[n]] + BItabsa[[2 + n]])*
\t  (-BKtabsab[[ n]] - BKtabsab[[2 + n]]))/
\t(4*BItabsab[[1+n]]) + 
       (-(sa*(-BKtabsa[[-1 + n]] - BKtabsa[[1+n]]))/2 - 
\t  (sa*(-BKtabsa[[1+n]] - BKtabsa[[3 + n]]))/2)/2) ),
    dbtermK[n,P] 
]
 
dnbtermI/: dnbtermI[n_,omega_,P_,d_] := 
\t(I^d)*D[btermI[n,w,P],{w,d}]/.w->omega


dnbtermK/: dnbtermK[n_,omega_,P_,d_] := 
\t(I^d)*D[btermK[n,w,P],{w,d}]/.w->omega


(*=========================================================================*)
dbi/: dbi[n_, x_] := D[BesselI[n, x], x]
dbif/: dbif[n_, y_] := dbi[n,x]/.x->y
mdbi/: mdbi[n_,x_] := (BesselI[n-1,x] + BesselI[n+1,x])/2

sbterm/: sbterm[btab_, M_, theta_,P_] :=
   Sum[btab[[i + 1]]*LegendreP[i, N[Cos[theta],P]], {i, 0, M}]

fastsbterm/: fastsbterm[btab_,Ltab_,M_] :=
   Sum[btab[[i + 1]]*Ltab[[i+1]], {i, 0, M}]

gf/: gf[tht_, omega_, N_, P_] := Sum[sterm[n, omega, tht, P], {n, 0, N}]

phase/: phase[tht_, omga_, N_, P_] :=
   ArcTan[Im[gf[tht, omga, N, P]]/Re[gf[tht, omga, N, P]]]

initab/: initab[N_] := Table[0, {i, N}]

wpN = wp[200,40];
wp/: wp[mhz_,P_] := N[ 2 Pi mhz / 1000000, P]
ang/: ang[Ma_,P_] := Table[ N[2 Pi (j-1)/(Ma-1), P],{j,Ma}]
angtab = ang[33,40];

(*=========================================================================*)
SetBesselTablesI/: SetBesselTablesI[ tab_, M_, x_, P_] :=
(*
\tEvaluate BesselI[n,x] for n = 0..M+1, by downward recursion
*)
       Module[ {},
       \ttab = initab[M+1];
       \ttab[[M+1]] = BesselI[M + 1/2, N[x, P] ];
\ttab[[M ]]  = BesselI[M - 1/2, N[x, P] ];
\tDo[
\t   {k = M-j; 
\t    tab[[k]] = (tab[[k+2]] + (2 k + 1)tab[[k+1]]/x) ;
\t    tab[[k]] = SetPrecision[tab[[k]],P]} ,
\t   {j,1,M-1}
\t];
       ]
 
SetBesselTablesK/: SetBesselTablesK[ tab_, M_, x_, P_] :=
(*
\tEvaluate BesselK[n,x] for n = 0..M+1, by upward recursion
*)
       Module[{},
\ttab = initab[M+1];
\ttab[[1]] = BesselK[1/2, N[x, P] ];
\ttab[[2]] = BesselK[3/2, N[x, P] ];
\tDo[
\t   {k = j+2; 
\t    tab[[k]] = (tab[[j]] + (2 j+1 )tab[[j+1]]/x);
\t    tab[[k]] = SetPrecision[tab[[k]],P]} ,
\t   {j,1,M-1}
\t];
       ]
 
 SetDBesselTablesI/: SetDBesselTablesI[ tab_, M_, x_, P_] :=
(*
\tEvaluate DBesselI[n,x] for n = 0..M+1,
\tFirst set up Bessel Tables, then use recurrence relation
*)
       Module[ {bbtab},
       \t(*bbtab = initab[M+2];*)
        SetBesselTablesI[ bbtab, M+2, x, P];
        tab = 1/2(RotateLeft[bbtab,1]+RotateRight[bbtab,1]);
        tab =Join[{DBessI[1/2,N[x,P]]}, Drop[Drop[tab,1],-1]]
       ]
       
  SetDBesselTablesK/:SetDBesselTablesK[tab_,M_,x_,P_]:=
  (*
  Evaluate DBesselI[n,x] for n=0..M+1,
  First set up Bessel Tables,then use recurrence relation
  *)
    Module[{bbtab},
    SetBesselTablesK[bbtab,M+2,x,P];
    tab=-1/2(RotateLeft[bbtab,1]+RotateRight[bbtab,1]);
    tab=Join[{DBessK[1/2,x]},Drop[Drop[tab,1],-1]]
    ]

SetCosTables/: SetCosTables[ tab_, M_, x_, P_] :=
       Module[{stab},
\ttab = initab[M+1];
\tstab = initab[M+1];
\ttab[[1]]  = 1;
\ttab[[2]]  = N[Cos[x], P];
\tstab[[1]]  = 0;
\tstab[[2]] = N[Sin[x], P];
\tDo[
\t   {
\t    tab[[k+1]] = tab[[k]]*tab[[2]] - stab[[k]]*stab[[2]];
\t    tab[[k+1]] = SetPrecision[tab[[k+1]],P];
\t    stab[[k+1]] = stab[[k]]*tab[[2]] + tab[[k]]*stab[[2]];
\t    stab[[k+1]] = SetPrecision[stab[[k+1]],P] },
\t   {k,2,M}
\t];
       ]  
SetLegendreTables/: SetLegendreTables[ tab_, M_, x_, P_] :=
       Module[{},
       \ttab = initab[M+1];
       \ttab[[1]]  = 1;
\ttab[[2]]  = N[x, P];
\tDo[
\t   {tab[[k+1]] = ((2 k -1) x tab[[k]] - (k -1)tab[[k-1]])/k ;
\t    tab[[k+1]] = SetPrecision[tab[[k+1]],P] },
\t   {k,2,M}
\t];
       ]
  
NextLegendre/: NextLegendre[x_,k_,lk_,lkm1_] :=
SetPrecision[((2k+1)x lk -k lkm1)/(k+1),Precision[x]]

NextLegendreStep/: NextLegendreStep[{x_,k_,lk_,lkm1_}] :=
{x,k+1,NextLegendre[x,k,lk,lkm1],lk}

MakeLegendreTable/: MakeLegendreTable[x_,n_] := Transpose[
NestList[NextLegendreStep,{x,1,x,1},n]][[4]]
(*=========================================================================*)
spherecalc::usage =
\"spherecalc[filename, M, fq, nang, P] calculates complex intensity\\n
on the surface of a sphere using parameters set by SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tfq is frequency in Mhz\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"

(*=========================================================================*)
spherecalc/: spherecalc[out_,M_, fq_, nang_,P_] :=
       Module[ {wpN, a0, btab, dbtab, ptab,Ltab},
       wpN = wp[fq,P];
       a0 = a0P[P];
       angtab = ang[2*nang+1,P];
       Print[\"Calculating  psi : \",M,\" terms \",P,\" precision \",fq,\" \
MHz\"];

       SetBesselTablesI[ BItabzeta , M +1, N[alpha[wpN]*zeta,P], P ];
       SetBesselTablesI[ BItabsa, M +1, N[alpha[wpN]*sa,P] ,P];
       btab = Table[fastbtermI[j - 1], {j, M+1}];
       Print[\"Finished btab\"];
       Clear[BItabzeta];
       Clear[BItabsa];
       pf5 = initab[nang*2 + 1];
       Print[\"Iteration  \\t angle \\t complex intensity \"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)
\t  pf5[[nang*2 + 2 - k]] = pf5[[k]] = fastsbterm[btab,Ltab,M];
\t Print[k, \"\\t\",N[180*angtab[[k]]/Pi],\"\\t\",N[pf5[[k]]]  ];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       af5 = -Arg[pf5];
       pf5 /= (2 Pi sa);
       pf5 /= 2 Sqrt[(2 Pi sa zeta)];
       ff5 = Abs[pf5];
       nangdeg = N[180*angtab/Pi, P];
       pf5N = N[pf5]; ff5N = N[ff5]; af5N = N[af5];
       Save[out, nangdeg,  pf5N, ff5N, af5N];
       {nangdeg,pf5,wpN}    (* return values *)
   ] (* End Module - spherecalc *)

(*=========================================================================*)
spherecalc1::usage =
\"spherecalc1[filename, M, nang, P] calculates integrated intensity and mean \
time\\n
on the surface of a sphere using parameters set by SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"

spherecalc1/: spherecalc1[out_, M_, nang_, P_] :=
       Module[ {wpN, a0, btab, dbtab,Ltab},
       a0 = a0P[P];
       angtab = ang[nang*2 + 1,P];
       Print[\"Calculating I <t> : \",nang,\" angles \",M,\" terms \",P,\" \
precision \"];
       SetBesselTablesI[ BItabzeta , M +1, N[alpha[0]*zeta,P], P ];
       SetBesselTablesI[ BItabsa, M +1, N[alpha[0]*sa,P] ,P];
       btab = Table[fastbtermI[j - 1], {j, M+1}];
       Print[\"Finished btab\"];
       dbtab = Table[fastdbtermI[j-1, P], {j, M+1}];
       Print[\"Finished dbtab\"];
       Clear[BItabzeta];
       Clear[BItabsa];
       if5 = initab[nang*2 + 1];
       tf5 = initab[nang*2 + 1];
       Print[\"Iteration  \\t angle \\t integrated intensity \\t \
time-weighted\"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)

\t  if5[[nang*2 + 2 - k]] = if5[[k]] = fastsbterm[btab,Ltab,M];
\t  tf5[[nang*2 + 2 - k]] = tf5[[k]] = fastsbterm[dbtab,Ltab,M];
\t Print[k, \" \",N[180*angtab[[k]]/Pi],\" \",N[if5[[k]]], \" \",N[tf5[[k]]] \
];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       if5 /= (2 Pi sa);
       tf5 /= (2 Pi sa);
       if5 /= 2 Sqrt[(2 Pi sa zeta)];
       tf5 /= 2 Sqrt[(2 Pi sa zeta)];
       nf5 = tf5/(if5 * 2 *a0 * gamma[P]^2);
       nangdeg = N[180*angtab/Pi, P];
       if5N = N[if5]; tf5N = N[tf5]; nf5N = N[nf5];
       Save[out, nangdeg, if5N , tf5N, nf5N];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - spherecalc1 *)

(*=========================================================================*)
extrapolatedspherecalc::usage =
\"extrapolatedspherecalc1[filename, M, fq, nang, P] calculates integrated \
intensity and mean time\\n
on the surface of a sphere using parameters set by SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tfq is frequency in Mhz\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"

extrapolatedspherecalc/: extrapolatedspherecalc[out_, M_, fq_, nang_, P_] :=
       Module [{wpN, aw, btab, dbtab,Ltab},
       wpN = wp[fq,P];
       aw = alpha[wpN,P];
       angtab = ang[nang*2 + 1,P];
       Print[\"Calculating complex I: \",nang,\" angles \",M,\" terms \",P,\" \
precision \"];
       SetBesselTablesI[ BItabzeta , M +2, N[alpha[wpN]*zeta,P], P ];
       SetBesselTablesI[ BItabsa, M +2, N[alpha[wpN]*sa,P] ,P];
       SetBesselTablesK[ BKtabsa, M +2, N[alpha[wpN]*sa,P] ,P];
       SetBesselTablesI[ BItabsab, M +2, N[alpha[wpN]*sab,P] ,P];
       SetBesselTablesK[ BKtabsab, M +2, N[alpha[wpN]*sab,P] ,P];
       btab = Table[fastbtermK[j - 1,wpN,P], {j, M+1}];    
       Print[\"Finished btab\"];
       Clear[BItabzeta];
       Clear[BItabsa];
       Clear[BKtabsa];
       Clear[BItabsab];
       Clear[BKtabsab];
       pf5 = initab[nang*2 + 1];
       Print[\"Iter.\\tangle\\tcomplex intensity\"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)
\t  pf5[[nang*2 +2 -k]] =pf5[[k]]= -1/(2 Pi) * fastsbterm[btab,Ltab,M];
\t Print[k, \"\\t\",N[180*angtab[[k]]/Pi],\"\\t\",N[pf5[[k]]] ];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       pf5 /= 2 Sqrt[(2 Pi sa zeta)];
       af5 = -Arg[pf5];
       ff5 = Abs[pf5];
       nangdeg = N[180*angtab/Pi, P];
       nnangdeg = N[180*angtab/Pi];
       epf5N = N[pf5]; eff5N = N[ff5]; eaf5N = N[af5];
       Save[out, nnangdeg, epf5N, eff5N, eaf5N];
       {nangdeg,pf5,wpN}    (* return values *)
   ] (* End Module - extrapolatedspherecalc *)
(*=========================================================================*)
extrapolatedspherecalc1::usage =
\"extrapolatedspherecalc1[filename, M, nang, P] calculates integrated \
intensity and mean time\\n
on the surface of a sphere using parameters set by SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"

extrapolatedspherecalc1/: extrapolatedspherecalc1[out_, M_, nang_, P_] :=
       Module [{a0, btab, dbtab,Ltab},
       a0 = a0P[P];
       angtab = ang[nang*2 + 1,P];
       Print[\"Calculating I <t> : \",nang,\" angles \",M,\" terms \",P,\" \
precision \"];
       SetBesselTablesI[ BItabzeta , M +2, N[alpha[0]*zeta,P], P ];
       SetBesselTablesI[ BItabsa, M +2, N[alpha[0]*sa,P] ,P];
       SetBesselTablesK[ BKtabsa, M +2, N[alpha[0]*sa,P] ,P];
       SetBesselTablesI[ BItabsab, M +2, N[alpha[0]*sab,P] ,P];
       SetBesselTablesK[ BKtabsab, M +2, N[alpha[0]*sab,P] ,P];
       btab = Table[fastbtermK[j-1, 0, P], {j, M+1}];     
       Print[\"Finished btab\"];
       dbtab = Table[fastdbtermK[j-1, P], {j, M+1}];
       Print[\"Finished dbtab\"]; 
       Clear[BItabzeta];
       Clear[BItabsa];
       Clear[BKtabsa];
       Clear[BItabsab];
       Clear[BKtabsab];
       if5 = initab[nang*2 + 1];
       tf5 = initab[nang*2 + 1];
       Print[\"Iteration  \\t angle \\t integrated intensity \\t \
time-weighted\"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)
\t  if5[[nang*2 +2 -k]] =if5[[k]]= -1/(2 Pi)* fastsbterm[btab,Ltab,M];
\t  tf5[[nang*2 +2 -k]] =tf5[[k]]= -1/(2 Pi)* fastsbterm[dbtab,Ltab,M];
\t Print[k,\"\\t\",N[180*angtab[[k]]/Pi],\"\\t\",N[if5[[k]]],\"\\t\",N[tf5[[k]\
]] ];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       if5 /= 2 Sqrt[(2 Pi sa zeta)];
       tf5 /= 2 Sqrt[(2 Pi sa zeta)];
       nf5 = tf5/if5;           (* scale factor already taken care of *) 
       nangdeg = N[180*angtab/Pi, P];
       eif5N = N[if5]; etf5N = N[tf5]; enf5N = N[nf5];
       Save[out, nangdeg, eif5N, etf5N, enf5N];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - extrapolatedspherecalc1 *)
(*=========================================================================*)
robinspherecalc::usage =
\"robinspherecalc[filename, M, fq, nang, P] calculates complex intensity\\n
on the surface of a sphere using Robin boundary conditions parameters set by \
SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tfq is frequency in Mhz\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"


(*=========================================================================*)
robinspherecalc/: robinspherecalc[out_,M_, fq_, nang_,P_] :=
(* UNDER CONSTRUCTION !*)
       Module[ {wpN, a0, btab, dbtab, ptab,Ltab,rcons},
       wpN = wp[fq,P];
       a0 = a0P[P];
       rcons = rA*N[alpha[wpN]]/(3*(mua + musd));
       angtab = ang[2*nang+1,P];
       Print[\"Sphere with Robin boundary condition \",rcons,\"\\nCalculating \
 psi : \",M,\" terms \",P,\" precision \",fq,\" MHz\"];
       SetBesselTablesI[ BItabzeta , M +1, N[alpha[wpN]*zeta,P], P ];
        Print[\"Finished BItabzeta\"];
       SetBesselTablesI[ BItabsa, M +1, N[alpha[wpN]*sa,P] ,P];
        Print[\"Finished BItabsa\"];
    (*   SetBesselTablesK[ BKtabsa, M +1, N[alpha[wpN]*sa,P] ,P];
        Print[\"Finished BKtabsa\"];
        *)
        
       SetDBesselTablesI[ BDItabsa, M +1, N[alpha[wpN]*sa,P] ,P];
        Print[\"Finished BDItabsa\"];
      (* SetDBesselTablesK[ BDKtabsa, M +1, N[alpha[wpN]*sa,P] ,P];
        Print[\"Finished BDKtabsa\"];
        *)
      (* btab = Table[BItabzeta[[j]]*BKtabsa[[j]]
       \t(1  - \
(1+rcons*BDKtabsa[[j]]/BKtabsa[[j]])/(1+rcons*BDItabsa[[j]]/BItabsa[[j]])), \
{j, M+1}];
       \t*)
       \t(*
       \tbtab = rA/(sa *3*(mua + musd)) Table[(2*j - \
1)*BItabzeta[[j]]/(BItabsa[[j]]+rcons*BDItabsa[[j]]), {j, M+1}];
       \t*)
       btab = Table[fastbtermIR[j - 1,rcons], {j, M+1}];
       
       Print[\"Finished btab\"];
      
       Clear[BItabzeta];
       Clear[BItabsa];
       Clear[BKtabsa];
       Clear[BDItabsa];
       Clear[BDKtabsa];                
       pf5 = initab[nang*2 + 1];
       Print[\"Iteration  \\t angle \\t complex intensity \"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)
\t  pf5[[nang*2 + 2 - k]] = pf5[[k]] = fastsbterm[btab,Ltab,M];
\t Print[k, \"\\t\",N[180*angtab[[k]]/Pi],\"\\t\",N[pf5[[k]]]  ];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       af5 = -Arg[pf5];
       pf5 /= (2 Pi sa);
       pf5 /= 2 Sqrt[(2 Pi sa zeta)];
       ff5 = Abs[pf5];
       nangdeg = N[180*angtab/Pi, P];
       pf5N = N[pf5]; ff5N = N[ff5]; af5N = N[af5];
       Save[out, nangdeg,  pf5N, ff5N, af5N];
       {nangdeg,pf5,wpN}    (* return values *)
   ] (* End Module - robinspherecalc *)

(*=========================================================================*)
NeumannSphereCalc1::usage =
\"NeumannSphereCalc1[filename, M, nang, P] calculates integrated intensity \
and \\n
mean time on the surface of a sphere using parameters set by \
SetParameters.\\n
\\tM is the number of terms used in the summation of the Bessel functions.\\n
\\tP is working precision.\\n
\\tnang is number of angles in half the circumference e.g if nang is 16, \\n
\\t\\t32 angles are output, plus 0 degrees.\\n
\\tfilename is an output file\\n
\\t\\tNOTE filename must be supplied with quotes, so that Mathematica\\n
\\t\\tinterprets it as a string.
\"

NeumannSphereCalc1/: NeumannSphereCalc1[out_, M_, nang_, P_] :=
       Module [{wpN, a0, btab, dbtab,Ltab},
       a0 = a0P[P];
       angtab = ang[nang*2 + 1,P];
       Print[\"Calculating I <t> : \",nang,\" angles \",M,\" terms \",P, \" \
precision \"];
(*       SetBesselTablesI[ BItabzeta , M +1, N[alpha[0]*zeta,P], P ]; *)
       SetBesselTablesI[ BItabsa, M +1, N[alpha[0]*sa,P] ,P];
       btab = Table[fastNeumannBterm[j - 1,P], {j, M+1}];
       Print[\"Finished btab\"];
(*       dbtab = Table[fastdbtermI[j-1, P], {j, M+1}]; *)
(*       Print[\"Finished dbtab\"]; *)
(*       Clear[BItabzeta]; *)
       Clear[BItabsa];
       if5 = initab[nang*2 + 1];
       tf5 = initab[nang*2 + 1]; 
       Print[\"Iteration  \\t angle \\t integrated intensity \\t \
time-weighted\"];
       Do[{
\t  SetLegendreTables[ Ltab, M, N[Cos[angtab[[k]] ] ,P],P]; 
(*\t  Ltab=MakeLegendreTable[N[Cos[angtab[[k]] ] ,P],M]; *)
\t  if5[[nang*2 + 2 - k]] = if5[[k]] = fastsbterm[btab,Ltab,M];
(*          tf5[[nang*2 + 2 - k]] = tf5[[k]] = fastsbterm[dbtab,Ltab,M]; *)
\t Print[k, \" \",N[180*angtab[[k]]/Pi],\" \",N[if5[[k]]], \" \",N[tf5[[k]]] \
];
\t  Clear[Ltab];
\t  }, {k, nang + 1}
       ];
       nf5 = tf5/(if5 * 2 *a0 * gamma[P]^2);
       if5 *= gamma[P]^2/(2 Pi sa);
       tf5 /= (2 Pi sa);
\t   nangdeg = N[180*angtab/Pi, P];
       nif5N = N[if5]; ntf5N = N[tf5]; nnf5N = N[nf5];
       Save[out, nangdeg, nif5N, ntf5N, nnf5N];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - NeumannSphereCalc1 *)

End[]   (* end of Private context *)
EndPackage[]
\
\>", "Input",
  PageWidth->Infinity,
  InitializationCell->True,
  ShowSpecialCharacters->False]