(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialization Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

(*      CircleGreensFunctions.m

	Author S.Arridge
	Date: September 1992

	last update : Nov 3rd 1995

	This file contains various Mathematica functions for computing
	Integrated intensity and mean time on boundary of a circle

	Note : extended to include extrapolated boundary condition
	This is now stable, although there is a scaling problem
	
	
    August 28th 2002:
             extension to Robin boundary conditions


*)
(*=========================================================================*)
BeginPackage["CircleGreensFunctions`"]

CircleGreensFunctions::usage =
"CircleGreensFunctions supplies various functions for calculating Green's
functions of the diffusion equation on an homogeneous circle\n
\n
Functions supplied (use ? for help) :\n
\tCircleIntegIntensityAndMT\n
\tCircleComplexData\n
\tExtrapolatedCircleIntegIntensityAndMT\n
\tExtrapolatedCircleComplexData\n
\tNeumannCircleIntegIntensityAndMT\n
\tListAngles\n
\tListIntegIntensity\n
\tListMeanTime\n
\tListComplexIntensity\n
\tListPhase\n
\tListModulation\n
\tDisplayIntegIntensity\n
\tDisplayMeanTime\n
\tDisplayPhase\n
\tDisplayModulation\n
\tComparePhaseAndMT\n
"
CircleIntegIntensityAndMT::usage =
"CircleIntegIntensityAndMT[filename, a, s, r, ref, nang, P] calculates integrated 
intensity and mean time on the surface of a circle.\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 circle 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=CircleIntegIntensityAndMT[\"circ.results\",\n
\t\t\t0.025, 2, 25, 1.4, 16, 50];
"

CircleIntegIntensityAndMT/: CircleIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	circcalc1[out,NTerm, nang, P]
    ]

CircleComplexData::usage =
"CircleComplexData[filename, a, s, r, ref, nang, mhz,P] calculates complex intensity on the surface of a circle.\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 circle 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=CircleComplexData[\"circ.results\",\n
	0.025, 2, 25, 1.4, 16, 200, 50];
"

CircleComplexData/: CircleComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	circcalc[out,NTerm, mhz, nang, P]
    ]

ExtrapolatedCircleIntegIntensityAndMT::usage =
"ExtrapolatedCircleIntegIntensityAndMT[filename, a, s, r, ref, nang, P] calculates \n
integrated intensity and mean time on the surface of a circle, 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 circle 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=ExtrapolatedCircleIntegIntensityAndMT[\"circ.results\",\n
\t\t\t0.025, 2, 25, 1.4, 16, 50];
"

ExtrapolatedCircleIntegIntensityAndMT/: ExtrapolatedCircleIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	extrapolatedcirccalc1[out,NTerm, nang, P]
    ]

ExtrapolatedCircleComplexData::usage =
"ExtrapolatedCircleComplexData[filename, a, s, r, ref, nang, mhz,P] calculates complex intensity on the surface of a circle.\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 circle 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 = ExtrapolatedCircleComplexData[\"circ.results\",\n
\t\t\t0.025, 2, 25, 1.4, 16, 200, 50];
"

ExtrapolatedCircleComplexData/: ExtrapolatedCircleComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	extrapolatedcirccalc[out,NTerm, mhz, nang, P]
    ]
NeumannCircleIntegIntensityAndMT::usage =
"NeumannCircleIntegIntensityAndMT[filename, a, s, r, ref, nang, P] calculates integrated 
intensity and mean time on the surface of a circle, 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 circle 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=NeumannCircleIntegIntensityAndMT[\"circ.results\",\n
\t\t\t0.025, 2, 25, 1.4, 16, 50];
"
RobinCircleComplexData::usage =
"RobinCircleComplexData[filename, a, s, r, ref, nang, mhz,P] calculates complex intensity on the surface of a circle.\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 = RobinCircleComplexData[\"robcircle.results\",\n
\t\t\t0.025, 2, 25, 1.4, 16, 200, 50];
"

RobinCircleComplexData/: RobinCircleComplexData[out_,a_,s_,r_,ref_,nang_,mhz_,  P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	robincirclecalc[out,NTerm, mhz, nang, P]
    ]
    
NeumannCircleIntegIntensityAndMT/: NeumannCircleIntegIntensityAndMT[out_,a_,s_,r_,ref_,nang_, P_] :=
    Module[{NTerm},
	SetParameters[a,s,r,ref,P];
	NTerm = Floor[(100*s*a*r*r/nang)]*nang;
	If[NTerm < 1536, NTerm = 1536];
	NeumannCircCalc1[out,NTerm, nang, P]
    ]

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

ListAngles/: ListAngles := CircleGreensFunctions`Private`angtab

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

ListIntegIntensity := CircleGreensFunctions`Private`if5

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

ListMeanTime := CircleGreensFunctions`Private`nf5

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

ListComplexIntensity := CircleGreensFunctions`Private`pf5

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

ListPhase := CircleGreensFunctions`Private`af5

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

ListModulation := Abs[CircleGreensFunctions`Private`pf5]

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

DisplayIntegIntensity :=
	Module[{pairs},
	pairs = Transpose[ {N[180/Pi ListAngles],Log[10,ListIntegIntensity]} ];
	ListPlot[pairs, PlotJoined->True, PlotLabel->"Log Intensity on circle",
		AxesLabel->{"Degrees","Log_10 I"}]
	]

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

DisplayMeanTime :=
	Module[{pairs},
	pairs = Transpose[ {N[180/Pi ListAngles],ListMeanTime} ];
	ListPlot[pairs, PlotJoined->True,PlotLabel->"Mean Time on circle",
		AxesLabel->{"Degrees","<t> picoseconds"}]
	]

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

DisplayPhase :=
	Module[{pairs},
	pairs = Transpose[ {N[180/Pi ListAngles],ListPhase} ];
	ListPlot[pairs, PlotJoined->True,PlotLabel->"Phase on circle",
		AxesLabel->{"Degrees","Arg[G] radians"}]
	]

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

DisplayModulation :=
	Module[{pairs},
	pairs = Transpose[ {N[180/Pi ListAngles],
			Log[10,Abs[ListComplexIntensity]]} ];
	ListPlot[pairs, PlotJoined->True,PlotLabel->"Modulation on circle",
		AxesLabel->{"Degrees","Abs[G] "}]
	]

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_] :=
	Module[{pairs1,pairs2,w,g1,g2},
	w = soln1[[3]];         (* frequency *)
	pairs1 = Transpose[ {soln1[[1]],-Arg[soln1[[2]]]/w} ];
	g1=ListPlot[pairs1, PlotJoined->True,PlotStyle->{RGBColor[1,0,0]},
		DisplayFunction->Identity];
	pairs2 = Transpose[ {soln2[[1]],soln2[[4]]} ];
	g2=ListPlot[pairs2, PlotJoined->True,PlotStyle->{RGBColor[0,1,0]},
		DisplayFunction->Identity];
	Show[g1,g2,PlotLabel->"Phase and MeanTime on circle",
		AxesLabel->{"Degrees","Phase,<t>"},
		DisplayFunction:>$DisplayFunction]
	]
(*=========================================================================*)
Begin["`Private`"]

SetParameters::usage = 
"SetParameters[a, s, r, ref, P] defines mua to be a, reduced mus to be s
a circle 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[{},
	mua  = SetPrecision[a,P];
	musd = SetPrecision[s,P];
	sa = SetPrecision[r,P];
	refind = SetPrecision[ref,P];
	zeta = sa - 1/musd;
	ic = 3/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);
	sab = sa + ebP[P];
	rA  = 2; (* Robin boundary condition, set to 2 for now *)
	Print["Setting mua : ",a," mus ",s," radius ",r," rind ",ref];
	a0 = 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] * 
	ArcTanh[ 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 ,alpha[omega,P]*zeta]*N[Cos[n* theta],P])/
  BesselI[n ,alpha[omega,P]*sa]

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

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

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

(*      here I am attempting to define things by replacement rules ! *)
(*
pkdef = (BesselI[n, alpha[omega,P]*zeta])*alpha[omega,P]*(D[BesselK[n, y],y] - 
	D[BesselI[n, y],y] * BesselK[n, alpha[omega,P]*sab]/
		BesselI[n, alpha[omega,P]*sab])/.{y->alpha[omega,P]*sa};
bessRules1 = 
{BesselI[n, alpha[omega,P]*zeta]->BItabzeta[[n + 1]],
 BesselI[n, alpha[omega,P]*sab]->BItabsab[[n+1]],
 BesselK[n, alpha[omega,P]*sab]->BKtabsab[[n+1]],
 D[BesselK[n, y],y]-> -(BKtabsa[[n]] + BKtabsa[[n+2]])/2,
 D[BesselI[n, y],y]->  (BItabsa[[n]] + BItabsa[[n+2]])/2,
};

btermK/: btermK[m_, w_, p_] := pkdef/.{n->m,omega->w,P->p}

fastbtermK/: fastbtermK[m_,p_] :=
   If[ m > 0,
    pkdef/.bessRules1/.{n->m,omega->w,P->p},
    btermK[m,0,p] 
   ]

*)

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

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

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

fastdbtermI/: fastdbtermI[n_, P_ ] :=
   If[ n > 0, 
   -( 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] 
    ]

(* This function got by differentiating btermK, w.r.t alpha. 
*)
dbtermK/: dbtermK[n_,P_] :=
(I*D[alpha[omega,P],omega]/.omega->0) * (BesselI[n, a0*zeta]*
     (-((BesselI[-1 + n, a0*sa] + BesselI[1 + n, a0*sa])*
	   BesselK[n, a0*sab])/(2*BesselI[n, a0*sab]) + 
       (-BesselK[-1 + n, a0*sa] - BesselK[1 + n, a0*sa])/2) + 
    (a0*zeta*(BesselI[-1 + n, a0*zeta] + BesselI[1 + n, a0*zeta])*
       (-((BesselI[-1 + n, a0*sa] + BesselI[1 + n, a0*sa])*
	     BesselK[n, a0*sab])/(2*BesselI[n, a0*sab]) + 
	 (-BesselK[-1 + n, a0*sa] - BesselK[1 + n, a0*sa])/2))/2 + 
    a0*BesselI[n, a0*zeta]*
     ((sab*(BesselI[-1 + n, a0*sa] + BesselI[1 + n, a0*sa])*
	  (BesselI[-1 + n, a0*sab] + BesselI[1 + n, a0*sab])*
	  BesselK[n, a0*sab])/(4*BesselI[n, a0*sab]^2) - 
       (((sa*(BesselI[-2 + n, a0*sa] + BesselI[n, a0*sa]))/2 + 
	    (sa*(BesselI[n, a0*sa] + BesselI[2 + n, a0*sa]))/2)*
	  BesselK[n, a0*sab])/(2*BesselI[n, a0*sab]) - 
       (sab*(BesselI[-1 + n, a0*sa] + BesselI[1 + n, a0*sa])*
	  (-BesselK[-1 + n, a0*sab] - BesselK[1 + n, a0*sab]))/
	(4*BesselI[n, a0*sab]) + 
       (-(sa*(-BesselK[-2 + n, a0*sa] - BesselK[n, a0*sa]))/2 - 
	  (sa*(-BesselK[n, a0*sa] - BesselK[2 + n, a0*sa]))/2)/2) )

fastdbtermK/: fastdbtermK[n_,P_] :=
(*
	just need to substitute in the tabulated Bessel functions
*)
If[ n > 1, 
 (I*D[alpha[omega,P],omega]/.omega->0) * (BItabzeta[[1+n]]*
     (-((BItabsa[[n]] + BItabsa[[2 + n]])*
	   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]])*
	     BKtabsab[[1+n]])/(2*BItabsab[[1+n]]) + 
	 (-BKtabsa[[n]] - BKtabsa[[2 + n]])/2))/2 + 
    a0*BItabzeta[[1+n]]*
     ((sab*(BItabsa[[n]] + BItabsa[[2 + n]])*
	  (BItabsab[[n]] + BItabsab[[2 + n]])*
	  BKtabsab[[1+n]])/(4*BItabsab[[1+n]]^2) - 
       (((sa*(BItabsa[[-1 + n]] + BItabsa[[1+n]]))/2 + 
	    (sa*(BItabsa[[1+n]] + BItabsa[[3 + n]]))/2)*
	  BKtabsab[[1+n]])/(2*BItabsab[[1+n]]) - 
       (sab*(BItabsa[[n]] + BItabsa[[2 + n]])*
	  (-BKtabsab[[ n]] - BKtabsab[[2 + n]]))/
	(4*BItabsab[[1+n]]) + 
       (-(sa*(-BKtabsa[[-1 + n]] - BKtabsa[[1+n]]))/2 - 
	  (sa*(-BKtabsa[[1+n]] - BKtabsa[[3 + n]]))/2)/2) ),
    dbtermK[n,P] 
]
 
dnbtermI/: dnbtermI[n_,omega_,P_,d_] := 
	(I^d)*D[btermI[n,w,P],{w,d}]/.w->omega


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

(*
NeumannBterm/: NeumannBterm[n_, omega_, P_] :=
	alpha[omega,P] *(mdbi[n , alpha[omega,P]*sa]/
		BesselI[n , alpha[omega,P]*sa] )

fastNeumannBterm/: fastNeumannBterm[n_, P_] :=
	If[ n > 0, 
		a0*(BItabsa[[n]]+BItabsa[[n + 2]])/(2 BItabsa[[n+1]]) ,
		NeumannBterm[n,0,P] 
	]
*)

NeumannBterm/: NeumannBterm[n_, omega_, P_] :=
	alpha[omega,P] *(mdbi[n , alpha[omega,P]*sa]/
		BesselI[n , alpha[omega,P]*sa] - n/alpha[omega,P]*sa)

fastNeumannBterm/: fastNeumannBterm[n_, P_] :=
	If[ n > 0, 
		a0*( (BItabsa[[n]]+BItabsa[[n + 2]])/(2 BItabsa[[n+1]]) - n/(a0*sa)),
		NeumannBterm[n,0,P] 
	]

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_] :=
   2*Sum[btab[[i + 1]]*N[Cos[i theta],P], {i, 1, M}] + btab[[1]]

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

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_] :=
(*
	Evaluate BesselI[n,x] for n = 0..M+1, by downward recursion
*)
       Module[ {},
	tab = initab[M+1];
	tab[[M+1]] = BesselI[M , N[x, P] ];
	tab[[M ]]  = BesselI[M-1, N[x, P] ];
	Do[
	   {k = M-j; 
	    tab[[k]] = (tab[[k+2]] + (2 k )tab[[k+1]]/x) ;
	    tab[[k]] = SetPrecision[tab[[k]],P]} ,
	   {j,1,M-1}
	];
       ]
 
SetBesselTablesK/: SetBesselTablesK[ tab_, M_, x_, P_] :=
(*
	Evaluate BesselK[n,x] for n = 0..M+1, by upward recursion
*)
       Module[{},
	tab = initab[M+1];
	tab[[1]] = BesselK[0, N[x, P] ];
	tab[[2]] = BesselK[1, N[x, P] ];
	Do[
	   {k = j+2; 
	    tab[[k]] = (tab[[j]] + (2 j )tab[[j+1]]/x);
	    tab[[k]] = SetPrecision[tab[[k]],P]} ,
	   {j,1,M-1}
	];
       ]
 
 
 SetDBesselTablesI/: SetDBesselTablesI[ tab_, M_, x_, P_] :=
(*
	Evaluate DBesselI[n,x] for n = 0..M+1,
	First set up Bessel Tables, then use recurrence relation
*)
       Module[ {bbtab},
       	(*bbtab = initab[M+2];*)
        SetBesselTablesI[ bbtab, M+2, x, P];
        tab = 1/2(RotateLeft[bbtab,1]+RotateRight[bbtab,1]);
        tab =Join[{DBessI[0,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[x]},Drop[Drop[tab,1],-1]]
    ]

 
SetCosTables/: SetCosTables[ tab_, M_, x_, P_] :=
       Module[{stab},
	tab = initab[M+1];
	stab = initab[M+1];
	tab[[1]]  = 1;
	tab[[2]]  = N[Cos[x], P];
	stab[[1]]  = 0;
	stab[[2]] = N[Sin[x], P];
	Do[
	   {
	    tab[[k+1]] = tab[[k]]*tab[[2]] - stab[[k]]*stab[[2]];
	    tab[[k+1]] = SetPrecision[tab[[k+1]],P];
	    stab[[k+1]] = stab[[k]]*tab[[2]] + tab[[k]]*stab[[2]];
	    stab[[k+1]] = SetPrecision[stab[[k+1]],P] },
	   {k,2,M}
	];
       ]  
(*=========================================================================*)
circcalc::usage =
"circcalc[filename, M, fq, nang, P] calculates complex intensity\n
on the surface of a circle 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.
"

(*=========================================================================*)
circcalc/: circcalc[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[{
	  SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
	  pf5[[nang*2 + 2 - k]] = pf5[[k]] = fastsbterm[btab,Ltab,M];
	 Print[k, "\t",N[180*angtab[[k]]/Pi],"\t",N[pf5[[k]]]  ];
	  Clear[Ltab];
	  }, {k, nang + 1}
       ];
       af5 = -Arg[pf5];
       pf5 /= (2 Pi sa);
       ff5 = Abs[pf5];
       nangdeg = N[180*angtab/Pi, P];
       Save[out, nangdeg,  pf5, ff5, af5];
       {nangdeg,pf5,wpN}    (* return values *)
   ] (* End Module - circcalc *)

(*=========================================================================*)
circcalc1::usage =
"circcalc1[filename, M, nang, P] calculates integrated intensity and mean time\n
on the surface of a circle 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.
"

circcalc1/: circcalc1[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[btermI[j - 1, 0, P], {j, M+1}]; *)
       btab = Table[fastbtermI[j - 1], {j, M+1}];
       Print["Finished btab"];
 (*      dbtab = Table[dbtermI[j - 1, P], {j, M+1}]; *)
       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[{
	  SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
	  if5[[nang*2 + 2 - k]] = if5[[k]] = fastsbterm[btab,Ltab,M];
	  tf5[[nang*2 + 2 - k]] = tf5[[k]] = fastsbterm[dbtab,Ltab,M];
	 Print[k, " ",N[180*angtab[[k]]/Pi]," ",N[if5[[k]]], " ",N[tf5[[k]]] ];
	  Clear[Ltab];
	  }, {k, nang + 1}
       ];
       nf5 = tf5/(if5 * 2 *a0 * gamma[P]^2);
       if5 /= (2 Pi sa);
       tf5 /= (2 Pi sa);
	   nangdeg = N[180*angtab/Pi, P];
       Save[out, nangdeg, if5 , tf5, nf5];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - circcalc1 *)

(*=========================================================================*)
extrapolatedcirccalc::usage =
"extrapolatedcirccalc1[filename, M, fq, nang, P] calculates integrated intensity and mean time\n
on the surface of a circle 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.
"

extrapolatedcirccalc/: extrapolatedcirccalc[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}];    
(*       btab = Table[btermK[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[{
	  SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
	  pf5[[nang*2 +2 -k]] =pf5[[k]]= -1/(2 Pi) * fastsbterm[btab,Ltab,M];
	 Print[k, "\t",N[180*angtab[[k]]/Pi],"\t",N[pf5[[k]]] ];
	  Clear[Ltab];
	  }, {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 - extrapolatedcirccalc *)
(*=========================================================================*)
extrapolatedcirccalc1::usage =
"extrapolatedcirccalc1[filename, M, nang, P] calculates integrated intensity and mean time\n
on the surface of a circle 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.
"

extrapolatedcirccalc1/: extrapolatedcirccalc1[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}];     
(*       btab = Table[btermK[j-1, 0, P], {j, M+1}]; *)
       Print["Finished btab"];
       dbtab = Table[fastdbtermK[j-1, P], {j, M+1}];
(*       dbtab = Table[dnbtermK[j -1, 0, P,1], {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[{
	  SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
	  if5[[nang*2 +2 -k]] =if5[[k]]= -1/(2 Pi)* fastsbterm[btab,Ltab,M];
	  tf5[[nang*2 +2 -k]] =tf5[[k]]= -1/(2 Pi)* fastsbterm[dbtab,Ltab,M];
	 Print[k,"\t",N[180*angtab[[k]]/Pi],"\t",N[if5[[k]]],"\t",N[tf5[[k]]] ];
	  Clear[Ltab];
	  }, {k, nang + 1}
       ];
       nf5 = tf5/if5;           (* scale factor already taken care of *) 
       nangdeg = N[180*angtab/Pi, P];
       Save[out, nangdeg, if5 , tf5, nf5];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - extrapolatedcircalc1 *)


(*=========================================================================*)
robincirclecalc::usage =
"robincirclecalc[filename, M, fq, nang, P] calculates complex intensity\n
on the surface of a circle 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.
"



(*=========================================================================*)
robincirclecalc/: robincirclecalc[out_,M_, fq_, nang_,P_] :=
       Module[ {wpN, a0, btab, dbtab, ptab,Ltab,rcons},
       wpN = wp[fq,P];
       a0 = a0P[P];
       rcons = rA*N[alpha[wpN]]*c/(3*(mua + musd));
       angtab = ang[2*nang+1,P];
       Print["Circle 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"];    
        
       SetDBesselTablesI[ BDItabsa, M +1, N[alpha[wpN]*sa,P] ,P];
       Print["Finished BDItabsa"];
      
       btab = Table[fastbtermIR[j - 1,rcons], {j, M+1}];
       
       Print["Finished btab"];
      
       Clear[BItabzeta];
       Clear[BItabsa];
   
       Clear[BDItabsa];
                
       pf5 = initab[nang*2 + 1];
       Print["Iteration  \t angle \t complex intensity "];
       Do[{               
			SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
            pf5[[nang*2 + 2 - k]] = pf5[[k]] = fastsbterm[btab,Ltab,M];
            Print[k, "\t",N[180*angtab[[k]]/Pi],"\t",N[pf5[[k]]]  ];
            Clear[Ltab];
          }, {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];
       rpf5N = N[pf5]; rff5N = N[ff5]; raf5N = N[af5];
       Save[out, nangdeg,  rpf5N, rff5N, raf5N];
       {nangdeg,pf5,wpN}    (* return values *)
   ] (* End Module - robinspherecalc *)


(*=========================================================================*)
NeumannCircCalc1::usage =
"NeumannCircCalc1[filename, M, nang, P] calculates integrated intensity and \n
mean time on the surface of a circle 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.
"

NeumannCircCalc1/: NeumannCircCalc1[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[{
	  SetCosTables[ Ltab, M, N[angtab[[k]] ,P],P]; 
	  if5[[nang*2 + 2 - k]] = if5[[k]] = fastsbterm[btab,Ltab,M];
(*          tf5[[nang*2 + 2 - k]] = tf5[[k]] = fastsbterm[dbtab,Ltab,M]; *)
	 Print[k, " ",N[180*angtab[[k]]/Pi]," ",N[if5[[k]]], " ",N[tf5[[k]]] ];
	  Clear[Ltab];
	  }, {k, nang + 1}
       ];
       nf5 = tf5/(if5 * 2 *a0 * gamma[P]^2);
       if5 *= gamma[P]^2/(2 Pi sa);
       tf5 /= (2 Pi sa);
	   nangdeg = N[180*angtab/Pi, P];
       Save[out, nangdeg, if5 , tf5, nf5];
       {nangdeg,if5,tf5,nf5}    (* return values *)
   ] (* End Module - NeumannCircCalc1 *)

End[]   (* end of Private context *)
EndPackage[]