(* SlabReflectGreensFunctions.m Author S.Arridge Date: September 1992 last update : May 11th 1993 This file contains various Mathematica functions for computing Integrated intensity and mean time on reflection boundary of a slab *) (*=========================================================================*) BeginPackage["SlabReflectGreensFunctions`"] SlabReflectGreensFunctions::usage = "SlabReflectGreensFunctions supplies various functions for calculating Green's functions of the diffusion equation on reflection side of homogeneous slab\n \n Functions supplied (use ? for help) :\n \tSlabBothData\n \tListPoints\n \tListIntegIntensity\n \tListMeanTime\n \tDisplayIntegIntensity\n \tDisplayMeanTime\n " SlabBothData::usage = "SlabBothData[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 reflection 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\tSlabBothData[\"rslab.results\",0.025, 2, 25, 1.4, 100, 21, 50] " SlabBothData/: SlabBothData[out_,a_,s_,r_,ref_,maxd_, np_, P_] := Module[{NTerm}, SetParameters[a,s,r,ref,P]; NTerm = Floor[(5*s)]; rslabcalc1[out,NTerm, maxd, np, P] ] ListPoints::usage = "ListPoints prints the points at which functions calculated" ListPoints/: ListPoints := SlabReflectGreensFunctions`Private`radtab ListIntegIntensity::usage = "ListIntegIntensity prints the integrated intensity at each point" ListIntegIntensity := SlabReflectGreensFunctions`Private`if5 ListMeanTime::usage = "ListMeanTime prints the mean time at each point" ListMeanTime := SlabReflectGreensFunctions`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] (*=========================================================================*) 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 gamma/: gamma[P_] := N[Sqrt[c/(3*(mua + musd))],P] gammaSym/: gammaSym := Sqrt[c/(3*(mua + musd))] ic = 3/10 refind = 14/10 c = ic/refind 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 )ds + z0 zm/: zm[n_] := (2 n )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] 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] ) gfslab/: gfslab[xi_,omega_,M_,P_] := z0(1 + alpha[omega,P] rhop[xi,0,P]) Exp[ -alpha[omega,P] rhop[xi,0,P] ]/ rhop[xi,0,P]^3 + Sum[pterm[xi,j,omega,P],{j,1,M}] ifslab/: ifslab[xi_,M_,P_] := z0(1 + a0 rhop[xi,0,P]) Exp[ -a0 rhop[xi,0,P] ]/ rhop[xi,0,P]^3 + Sum[bterm[xi,j,P],{j,1,M}] tfslab/: tfslab[xi_,M_,P_] := z0 Exp[ -a0 rhop[xi,0,P] ]/rhop[xi,0,P] + Sum[dterm[xi,j,P],{j,1,M}] rad/: rad[maxrad_,Ma_,P_] := Table[ N[(j-1)maxrad/(Ma-1), P],{j,Ma}] radtab = rad[100,21,50] (*==========================================================================*) rslabcalc/: rslabcalc[M_, P_, fq_] := { 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]; pf5 = initab[21]; 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, 21} ]; nf5 = tf5/(if5 * 2 gamma[P]^2); af5 = -Arg[pf5]; Save["rslab6a.out",if5 , tf5, pf5, nf5, af5 ]; pf0 = initab[21]; 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["rslab6b.out",wtab,pf0,af0,mf0]; } (*==========================================================================*) rslabcalc1::usage = "rslabcalc1[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 reflection 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 " rslabcalc1/: rslabcalc1[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 ]; } (*==========================================================================*) rslabcalc2/: rslabcalc2[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["rslab6a2.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["rslab6.moments",nmoment]; },{j,mnts} ]; } End[] (* end of Private context *) EndPackage[]