(* SlabTransGreensFunctions.m Author S.Arridge Date: September 1992 last update : May 3rd 1995 This file contains various Mathematica functions for computing Integrated intensity and mean time on transmission boundary of a slab Notes : 1) Now includes option to use Neumann source. The mean time works beautifully. Integrated Intensity includes a strange scaling that I can't get to the bottom of. *) (*=========================================================================*) BeginPackage["SlabTransGreensFunctions`"] SlabTransGreensFunctions::usage = "SlabTransGreensFunctions supplies various functions for calculating Green's functions of the diffusion equation on transmission side of homogeneous slab\n \n Functions supplied (use ? for help) :\n \tSlabIntegIntensityAndMTData\n \tSlabComplexData\n \tNeumannSlabIntegIntensityAndMTData\n \tListPoints\n \tListIntegIntensity\n \tListMeanTime\n \tListComplexIntensity\n \tListPhase\n \tDisplayIntegIntensity\n \tDisplayMeanTime\n \tDisplayPhase\n " SlabIntegIntensityAndMTData::usage = "SlabIntegIntensityAndMTData[filename, a, s, d, ref, maxd, np, P] calculates integrated intensity and mean time on the surface of a slab.\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 \td is thickness of slab in mm\n \tref is refractive index\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tP is working precision. \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 \tExample:\n \t\tSlabIntegIntensityAndMTData[\"slab.results\",0.025, 2, 25, 1.4, 100, 21, 16] " SlabIntegIntensityAndMTData/: SlabIntegIntensityAndMTData[out_,a_,s_,r_,ref_,maxd_, np_, P_] := Module[{NTerm}, SetParameters[a,s,r,ref,P]; NTerm = Floor[(5*s)]; slabcalc1[out,NTerm, maxd, np, P] ] SlabComplexData::usage = "SlabComplexData[filename, a, s, d, ref, maxd, np,mhz, P] calculates complex intensity on the surface of a slab.\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 \td is thickness of slab in mm\n \tref is refractive index\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tmhz is frequency in mhz\n \tP is working precision. \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 \tExample:\n \t\tSlabComplexData[\"slab.results\",0.025, 2, 25, 1.4, 100, 21,200,16] " SlabComplexData/: SlabComplexData[out_,a_,s_,r_,ref_,maxd_,np_,mhz_, P_] := Module[{NTerm}, SetParameters[a,s,r,ref,P]; NTerm = Floor[(5*s)]; slabcalc[out,NTerm,mhz, maxd, np, P] ] NeumannSlabIntegIntensityAndMTData::usage = "NeumannSlabIntegIntensityAndMTData[filename, a, s, d, ref, maxd, np, P] calculates integrated intensity and mean time on the surface of a slab, assuming 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 \td is thickness of slab in mm\n \tref is refractive index\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tP is working precision. \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 \tExample:\n \t\tNeumannSlabIntegIntensityAndMTData[\"neumann_slab.results\",0.025, 2, 25, 1.4, 100, 21,50] " NeumannSlabIntegIntensityAndMTData/: NeumannSlabIntegIntensityAndMTData[out_,a_,s_,r_,ref_,maxd_, np_, P_] := Module[{NTerm}, SetParameters[a,s,r,ref,P]; NTerm = Floor[(5*s)]; NeumannSlabCalc1[out,NTerm, maxd, np, P] ] ListPoints::usage = "ListPoints prints the points at which functions calculated" ListPoints/: ListPoints := SlabTransGreensFunctions`Private`radtab ListIntegIntensity::usage = "ListIntegIntensity prints the integrated intensity at each point" ListIntegIntensity := SlabTransGreensFunctions`Private`if5 ListComplexIntensity::usage = "ListComplexIntensity prints the complex intensity at each point" ListComplexIntensity := SlabTransGreensFunctions`Private`pf5 ListPhase::usage = "ListPhase prints the phase at each point" ListPhase := SlabTransGreensFunctions`Private`af5 ListMeanTime::usage = "ListMeanTime prints the mean time at each point" ListMeanTime := SlabTransGreensFunctions`Private`nf5 DisplayIntegIntensity::usage = "DisplayIntegIntensity plots the log of integrated intensity at each point" DisplayIntegIntensity := ListPlot[Log[10,ListIntegIntensity], PlotJoined->True] DisplayMeanTime::usage = "DisplayMeanTime plots the mean time at each point" DisplayMeanTime := ListPlot[ListMeanTime, PlotJoined->True] DisplayPhase::usage = "DisplayPhase plots the phase at each point" DisplayPhase := ListPlot[ListPhase, PlotJoined->True] (*=========================================================================*) Begin["`Private`"] SetParameters::usage = "SetParameters[a, s, d, ref, P] defines mua to be a, reduced mus to be s a slab of width ds 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]; ds = SetPrecision[r,P]; refind = SetPrecision[ref,P]; z0 = 1/musd; c = ic/refind; a0 = a0P[P]; ] (*=========================================================================*) mua = 4/100; musd = 9; ds = 15; z0 = 1/musd; ic = 3/10; refind = 14/10; c = ic/refind; gamma/: gamma[P_] := N[Sqrt[c/(3*(mua + musd))],P] gammaSym/: gammaSym := Sqrt[c/(3*(mua + musd))] kappa/: kappa := c/(3*(mua + musd)) 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] a0 = a0P[40]; initab/: initab[N_] := Table[0, {i, N}] wpN = wp[200,40]; wp/: wp[mhz_,P_] := N[ 2 Pi mhz / 1000000, P] (*=========================================================================*) (* Following are specific to a slab *) zp/: zp[n_] := (2 n + 1)ds + z0 zm/: zm[n_] := (2 n + 1)ds - z0 rhop/: rhop[xi_,n_,P_] := N[Sqrt[xi^2 + zp[n]^2],P] rhom/: rhom[xi_,n_,P_] := N[Sqrt[xi^2 + zm[n]^2],P] (* repeat with source as a parameter *) zp/: zp[zz_,n_] := (2 n + 1)ds + zz zm/: zm[zz_,n_] := (2 n + 1)ds - zz rhop/: rhop[zz_,xi_,n_,P_] := N[Sqrt[xi^2 + zp[zz,n]^2],P] rhom/: rhom[zz_,xi_,n_,P_] := N[Sqrt[xi^2 + zm[zz,n]^2],P] bterm/: bterm[zz_,xi_,n_,P_] := -( zp[zz,n](1 + a0 rhop[zz,xi,n,P])Exp[ -a0 rhop[zz,xi,n,P] ]/ rhop[zz,xi,n,P]^3 - zm[zz,n](1 + a0 rhom[zz,xi,n,P])Exp[ -a0 rhom[zz,xi,n,P] ]/ rhom[zz,xi,n,P]^3 ) dterm/: dterm[zz_, xi_,n_,P_] := -( zp[zz,n](Exp[ -a0 rhop[zz,xi,n,P] ] )/rhop[zz,xi,n,P] - zm[zz,n](Exp[ -a0 rhom[zz,xi,n,P] ] )/rhom[zz,xi,n,P] ) (* repeat as analytical expressions *) Arhop/: Arhop[zz_,xi_,n_] := Sqrt[xi^2 + zp[zz,n]^2] Arhom/: Arhom[zz_,xi_,n_] := Sqrt[xi^2 + zm[zz,n]^2] Abterm/: Abterm[zz_,xi_,n_] := -( zp[zz,n](1 + a0 Arhop[zz,xi,n])Exp[ -a0 Arhop[zz,xi,n] ]/ Arhop[zz,xi,n]^3 - zm[zz,n](1 + a0 Arhom[zz,xi,n])Exp[ -a0 Arhom[zz,xi,n] ]/ Arhom[zz,xi,n]^3 ) Adterm/: Adterm[zz_, xi_,n_] := -( zp[zz,n](Exp[ -a0 Arhop[zz,xi,n] ] )/Arhop[zz,xi,n] - zm[zz,n](Exp[ -a0 Arhom[zz,xi,n] ] )/Arhom[zz,xi,n] ) bterm/: bterm[xi_,n_,P_] := -( zp[n](1 + a0 rhop[xi,n,P])Exp[ -a0 rhop[xi,n,P] ]/rhop[xi,n,P]^3 - zm[n](1 + a0 rhom[xi,n,P])Exp[ -a0 rhom[xi,n,P] ]/rhom[xi,n,P]^3 ) pterm/: pterm[xi_,n_,omega_,P_] := -( zp[n](1 + alpha[omega,P] rhop[xi,n,P]) * Exp[ -alpha[omega,P] rhop[xi,n,P] ] /rhop[xi,n,P]^3 - zm[n](1 + alpha[omega,P] rhom[xi,n,P] ) * Exp[ -alpha[omega,P] rhom[xi,n,P] ] /rhom[xi,n,P]^3 ) dterm/: dterm[xi_,n_,P_] := -( zp[n](Exp[ -a0 rhop[xi,n,P] ] )/rhop[xi,n,P] - zm[n](Exp[ -a0 rhom[xi,n,P] ] )/rhom[xi,n,P] ) dpterm/: dpterm[xi_,n_,omega_,P_,d_] := (I^d)*D[pterm[xi,n,w,P],{w,d}]/.w->omega gfslab/: gfslab[xi_,omega_,M_,P_] := Sum[pterm[xi,j,omega,P],{j,0,M}] ifslab/: ifslab[xi_,M_,P_] := Sum[bterm[xi,j,P],{j,0,M}] tfslab/: tfslab[xi_,M_,P_] := Sum[dterm[xi,j,P],{j,0,M}] dgfslab/: dgfslab[xi_,omega_,M_,P_,d_] := Sum[dpterm[xi,j,omega,P,d],{j,0,M}] (* Neumann versions of things *) dsn/: dsn[n_] := ds*(2 n +1) NeumRhoSq/: NeumRhoSq[xi_,n_] := (dsn[n]^2 + xi^2) NeumRho/: NeumRho[xi_,n_,P_] := N[Sqrt[dsn[n]^2 + xi^2],P] NeumannBterm/: NeumannBterm[xi_,n_] := Exp[-a0*NeumRhoSq[xi,n]^(1/2)]*kappa* ( (-2*a0*dsn[n]^2)/NeumRhoSq[xi,n]^2 + (6*dsn[n]^2*(1 + a0*NeumRhoSq[xi,n]^(1/2)))/NeumRhoSq[xi,n]^(5/2) + (2*a0*dsn[n]^2*(1 + a0*NeumRhoSq[xi,n]^(1/2)))/NeumRhoSq[xi,n]^2 - (2*(1 + a0*NeumRhoSq[xi,n]^(1/2)))/NeumRhoSq[xi,n]^(3/2) ) NeumannBterm/: NeumannBterm[xi_,n_,P_] := Exp[-a0*NeumRho[xi,n,P]]*gamma[P]^2* ( (-2*a0*dsn[n]^2)/(NeumRhoSq[xi,n]^2) + (6*dsn[n]^2*(1 + a0*NeumRho[xi,n,P]))/(NeumRho[xi,n,P]^(5)) + (2*a0*dsn[n]^2*(1 + a0*NeumRho[xi,n,P]))/(NeumRhoSq[xi,n]^2) - (2*(1 + a0*NeumRho[xi,n,P]))/(NeumRho[xi,n,P]^(3)) ) NeumannDterm/: NeumannDterm[xi_,n_] := Exp[-(a0*NeumRhoSq[xi,n]^(1/2))]*kappa* ( (2*dsn[n]^2)/(NeumRhoSq[xi,n]^(3/2)) + (2*a0*dsn[n]^2)/(NeumRhoSq[xi,n]) - 2/(NeumRhoSq[xi,n]^(1/2)) ) NeumannDterm/: NeumannDterm[xi_,n_,P_] := Exp[-a0*NeumRho[xi,n,P]]*gamma[P]^2* ( (2*dsn[n]^2)/(NeumRho[xi,n,P]^(3)) + (2*a0*dsn[n]^2)/(NeumRhoSq[xi,n]) - 2/(NeumRho[xi,n,P]) ) NeumannPterm/: NeumannPterm[xi_,n_,omega_, P_] := Exp[-(alpha[omega,P]*NeumRho[xi,n,P]) ]*gamma[P]^2 * ( (-2*alpha[omega,P]*dsn[n]^2)/(NeumRhoSq[xi,n]^2) + (6*dsn[n]^2*(1 + alpha[omega,P]*( dsn[n]^2 + xi^2)^(1/2)))/ (NeumRho[xi,n,P]^(5)) + (2*alpha[omega,P]*dsn[n]^2*(1 + alpha[omega,P]*NeumRho[xi,n,P]))/ (NeumRhoSq[xi,n]^2) - (2*(1 + alpha[omega,P]*NeumRho[xi,n,P]))/(NeumRho[xi,n,P]^(3)) ) NeumannGfslab/: NeumannGfslab[xi_,omega_,M_,P_] := Sum[NeumannPterm[xi,j,omega,P],{j,0,M}] NeumannIfslab/: NeumannIfslab[xi_,M_] := Sum[NeumannBterm[xi,j],{j,0,M}] NeumannIfslab/: NeumannIfslab[xi_,M_,P_] := Sum[NeumannBterm[xi,j,P],{j,0,M}] NeumannTfslab/: NeumannTfslab[xi_,M_] := Sum[NeumannDterm[xi,j],{j,0,M}] NeumannTfslab/: NeumannTfslab[xi_,M_,P_] := Sum[NeumannDterm[xi,j,P],{j,0,M}] rad/: rad[maxrad_,Ma_,P_] := Table[ N[(j-1)maxrad/(Ma-1), P],{j,Ma}] radtab = rad[100,21,50]; (*=========================================================================*) slabcalc::usage = "slabcalc[filename, M, fq,maxd, np, P] calculates integrated intensity and \n mean time and phase on the surface of a slab using parameters set by SetParameters.\n \tM is the number of terms used in the summation of the Bessel functions.\n \tfq is frequency in Mhz\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tP is working precision. \n \tfilename is an output file - NOTE filename must be supplied with\n \t\tquotes, so that Mathematica interprets it as a string.\n " slabcalc/: slabcalc[out_,M_, fq_,maxd_,np_,P_] := { wpN = wp[fq,P]; a0 = a0P[P]; radtab = rad[maxd,np,P]; Print["Finding I phi: ",M," terms ",P," precision ",fq," MHz"]; if5 = initab[np]; tf5 = initab[np]; pf5 = initab[np]; Print["Iteration integrated intensity FT of TPSF "]; Do[{ if5[[k]] = ifslab[radtab[[k]],M,P]; tf5[[k]] = tfslab[radtab[[k]],M,P]; pf5[[k]] = gfslab[radtab[[k]],wpN,M,P]; Print[k, " ", if5[[k]], " ", tf5[[k]], " ", pf5[[k]]]}, {k, np} ]; nf5 = tf5/(if5 * 2 gamma[P]^2); af5 = -Arg[pf5]; Save[out,if5 , tf5, pf5, nf5, af5 ]; pf0 = initab[np]; wtab = Table[ wp[50*i,P],{i,20}]; Do[{ pf0[[k+1]] = gfslab[0,wtab[[k]],M,P]; Print[wtab[[k]]," ",pf0[[k+1]] ]},{k,20} ]; pf0[[1]] = if5[[1]]; af0 = -Arg[pf0]; mf0 = Abs[pf0]/pf0[[1]]; Save[out,wtab,pf0,af0,mf0]; } (*==========================================================================*) slabcalc1::usage = "slabcalc1[filename, M, maxd, np, P] calculates integrated intensity and \n mean time on the surface of a slab using parameters set by SetParameters.\n \tM is the number of terms used in the summation of the Bessel functions.\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tP is working precision. \n \tfilename is an output file - NOTE filename must be supplied with\n \t\tquotes, so that Mathematica interprets it as a string.\n " slabcalc1/: slabcalc1[out_, M_, maxd_, np_, P_] := { wpN = wp[fq,P]; a0 = a0P[P]; radtab = rad[maxd,np,P]; Print["Finding I : ",M," terms ",P," precision ",np," points"]; if5 = initab[np]; tf5 = initab[np]; pf5 = initab[np]; Print["Iteration point integrated intensity time-weighted"]; Do[{ if5[[k]] = ifslab[radtab[[k]],M,P]; tf5[[k]] = tfslab[radtab[[k]],M,P]; Print[k, " ", radtab[[k]], " ", if5[[k]], " ", tf5[[k]] ]}, {k, np} ]; if5 /= (2 Pi); tf5 /= (2 Pi); nf5 = tf5/(if5 * 2 gamma[P]^2); Save[out,if5 ,tf5, nf5 ]; } (*==========================================================================*) NeumannSlabCalc1::usage = "NeumannSlabCalc1[filename, M, maxd, np, P] calculates integrated intensity\n and mean time on the surface of a slab using parameters set by SetParameters.\n \tM is the number of terms used in the summation of the Bessel functions.\n \tmaxd is maximum radius from source that is required, \n \tnp is number of points sampled on transmission side, \n \tP is working precision. \n \tfilename is an output file - NOTE filename must be supplied with\n \t\tquotes, so that Mathematica interprets it as a string.\n " NeumannSlabCalc1/: NeumannSlabCalc1[out_, M_, maxd_, np_, P_] := { wpN = wp[fq,P]; a0 = a0P[P]; radtab = rad[maxd,np,P]; Print["Finding I : ",M," terms ",P," precision ",np," points"]; if5 = initab[np]; tf5 = initab[np]; pf5 = initab[np]; Print["Iteration point integrated intensity time-weighted"]; Do[{ if5[[k]] = NeumannIfslab[radtab[[k]],M,P]; tf5[[k]] = NeumannTfslab[radtab[[k]],M,P]; Print[k, " ", radtab[[k]], " ", if5[[k]], " ", tf5[[k]] ]}, {k, np} ]; if5 /= (2 Pi); tf5 /= (2 Pi); nf5 = tf5/(if5 * 2 gamma[P]^2); Save[out,if5 ,tf5, nf5 ]; } (*==========================================================================*) slabcalc2/: slabcalc2[M_, P_, fq_, mnts_] := { wpN = wp[fq,P]; a0 = a0P[P]; radtab = rad[100,21,P]; Print["Finding I phi: ",M," terms ",P," precision ",fq," MHz"]; if5 = initab[21]; tf5 = initab[21]; Print["Iteration integrated intensity FT of TPSF "]; Do[{ if5[[k]] = ifslab[radtab[[k]],M,P]; tf5[[k]] = tfslab[radtab[[k]],M,P]; Print[k, " ", if5[[k]], " ", tf5[[k]] ]}, {k, 21} ]; nf5 = tf5/(if5 * 2 gamma[P]^2); Save["slab6a2.out",if5 , tf5, nf5 ]; Do[ { moments = initab[21]; Do[{ moments[[k]] = dgfslab[radtab[[k]],0,M,P,j]; Print[k, " ", if5[[k]], " ", moments[[k]]]}, {k, 21} ]; nmoment = moments/if5; Save["slab6.moments",nmoment]; },{j,mnts} ]; } End[] (* end of Private context *) EndPackage[]