montecarlomaplet.maplet

A Monte-Carlo Maplet for the Study of the Optical Properties of Biological Tissues

# Source code for MonteCarloMaplet
# Authors: Man Ho Yip and M.J. Carvalho
# Ms. Ref. No.:  CPC-D-07-00087
# I - Monte Carlo Simulation
# 1. The main simulation procedure
restart;with(plots):
Digits := 20: MonteCarloInMaple:=proc(graphs,geom,photons,entryang)
global layers,n,mua,mus,mut,g,d,z0,z1,w,A,Rd,Td,TRsp,layer,sleft,uz,tempRd,tempTd,rmax,ymax,xmax,RTT,ATT,TTR,TTA,TT,RR,DW;
local i,TR,TD,RD,AA,PPa,PPrT,PPrR,PPaT,PPaR;
TTR:=[NULL];TTA:=[NULL];TT:=[NULL];
RTT:=[NULL];ATT:=[NULL];RR:=[NULL];


simprop():
 
  for i from 0 to layers do # Initialize to be 0 in order to setup bounds
      z0[i] := 0:
      z1[i] := 0:
  od:

tissues():

A := 0: Rd := 0: Td := 0: TRsp := 0:
 
# Launching  photon
  for i from 1 to photons do
      
initialcond():
initialangle(entryang):


      layer := 1:       # currently in the first layer
      specular():
      while w > wmin do
            stepsize():
            finddistance():
if geom=cylindrical then 
                if r > rmax then w := 0 fi:
elif geom=prismatic then          
             if x > xmax or y > ymax then w := 0 fi:
fi;
            if s >= db then
               sleft := (s-db)*mut[layer]:
               angles(): fresnel():
               if xi() <= R then
                # reflect
                  uz := -1*uz:
               else
                # transmit 
                  movetoboundary():
                  if z = z0[0] then                     ### leaving through the top
                     # diffuse reflectance
                       # Some w reflected (reduces deviations)
                         fresnel():
                     tempRd := (1-R)*w: Rd := Rd+tempRd: w := w*R: uz := -1*uz;
RTT:=[op(RTT),r];
ATT:=[op(ATT),alphat];                 

# print("alpha=0",i,w,tempRd);
RR:=[op(RR),tempRd];
                             
tempRd := 0: tempTd := 0:
                  elif z = z1[layers] then                  ####  end of layers, exiting
                     # diffuse transmittance
                       # Some w transmitted
                         fresnel():
                         tempTd := (1-R)*w; Td := Td+tempTd; 
TTR:=[op(TTR),r];
TTA:=[op(TTA),alphat];

TT:=[op(TT),tempTd];#print(TT);
                         w := w*R; uz := -1*uz;
                     tempRd := 0: tempTd := 0:
                  else                                 ####  transmitted into next layer
                          newstep():
                       if uz > 0 then
                         layer := layer + 1:
                       elif uz < 0 then
                         layer := layer - 1:
                       fi:
                  fi:
                   # Reset temp values
               fi:
            else
               interact(): # Updates directions and weights

               DRabsorb(): scatter():
            fi:
            if w <> 0 then
             
               roulette():   # Photon Weight Roulette
            fi:
      od:
  od:
TR:=TRsp/photons;RD:=Rd/photons;AA:=A/photons;TD:=Td/photons;
save TR,RD,AA,TD, `MCresults.m`:
reflect();
transmit();

          if graphs=1 then   ##### compute graph points

PPa:=absorbgraph(photons):
save PPa,`MCgraphA.m`:

PPrT:=radTgraph(photons):  PPrR:=radRgraph(photons):
save PPrT,PPrR, `MCgraphRad.m`:

PPaT:=angTgraph(photons):  PPaR:=angRgraph(photons):
save PPaT,PPaR, `MCgraphAng.m`:

          fi;
end:

# 2. Procedures for the simulation

   #1. Load Simulation properties 
   simprop := proc()
   local i:
   global ARdv, ARd, ATdv, ATd, RRdv, RRd, RTdv, RTd, DRAv, DRA, Na, Nr, dr,  Nz, dz, wmin, m, da:
   # Parameters for Angular Resolution. Can be changed in INPUT

      Na := 30:                  # INPUT number of grid elements for exit angle, alphat
      da := Pi/(2*Na):     # delta_alpha, angle resolution                       
      ARdv := array(0..Na-1): ARd := array(0..Na-1):
      ATdv := array(0..Na-1): ATd := array(0..Na-1):
      for i from 0 to Na-1 do
        ARdv[i] := 0;             # Rd values initialized to 0
        ARd[i] := (i+0.5)*da; #(1/2)*(i*Pi/(2*Na)+(i+1)*Pi/(2*Na));
          ATdv[i] := 0;           # Td values initialized to 0
          ATd[i] := (i+0.5)*da;#(1/2)*(i*Pi/(2*Na)+(i+1)*Pi/(2*Na));
      od:

   # Parameters for Radial Resolution. Can be changed in INPUT
      Nr := 200:                    # INPUT number of radial grid elements
      dr := 0.05:            # INPUT radial resolution [cm]        
      RRdv := array(0..Nr-1): RRd := array(0..Nr-1):
      RTdv := array(0..Nr-1): RTd := array(0..Nr-1):
      for i from 0 to Nr-1 do
        RRdv[i] := 0;              # Rd values initialized to 0
        RRd[i] := (i+0.5)*dr;
        RTdv[i] := 0;              # Td values initialized to 0
        RTd[i] := (i+0.5)*dr;
      od:
   # Parameters for Depth Absorption Resolution
   
      Nz := 20;    # *INPUT* number of depth grid elements  
      dz := totald/Nz;     # *INPUT* depth resolution where totald is the sum of the thickness of all layers          
      DRAv := array(0..Nz-1): DRA := array(0..Nz-1):
      for i from 0 to Nz-1 do
        DRA[i] := (i+0.5)*dz:
        DRAv[i] := 0:
      od:
   # Photon Properties
      wmin := 0.0001: m := 10:
  end:

#2.  Load random number generator
  xi := proc()
    evalf(rand()/999999999999);
  end:

#3.  Load initial conditions
  initialcond := proc()
      global s, sleft, w,  x, y, z, Rsp: 
# step size, remaining step, photon weight, x,y,z) co-ordinates, specular reflectance

         x := 0: y := 0: z := 0:            
             s := 0: 
         sleft := 0: 
             w := 1: 
         Rsp := 0: 
          
  end:
#4. Initial directions for launch angle  
  initialangle := proc(langle)
  local temp;
  global ux, uy, uz:  # (ux, uy, uz) angle cossines
    temp := 90 - langle: # angle relative to surface's normal
    temp := evalf(temp*Pi/180): # convert to radians
    uz := cos(temp):
    ux := sqrt((1-uz*uz)/2):
    uy := ux:	#  ux=uy since ux, uy are assumed arbitrary
end:

#5. Specular Reflection
  specular := proc()
    global Rsp, TRsp, w: 
    Rsp := ((n[0]-n[1])^2)/((n[0]+n[1])^2): 
    TRsp := TRsp + Rsp: 
    w := 1-Rsp:
  end:

#6. Step size calculation
  stepsize := proc()
    global s, sleft:
    if sleft = 0 then
        s := -ln(xi())/mut[layer]:
    else
        s := sleft/mut[layer]:
        sleft := 0:
    fi:
  end:

#7. Find distance calculation
  finddistance := proc()
    global db, r:
      if uz > 0 then
         db := (z1[layer]-z)/uz:      # from lower layer
      elif uz < 0 then
         db := (z0[layer]-z)/uz:      # from upper layer
      fi:
        r := sqrt(x*x+y*y):           # radial position
  end:

#8. Move photon to boundary
  movetoboundary := proc()
    global x, y, z, w, layer;
      if uz > 0 then  	# going down
         x := x + ux*db;
         y := y + uy*db;
         z := z1[layer];
      elif uz < 0 then	# going up
         x := x + ux*db;
         y := y + uy*db;
         z := z0[layer];
      else		  # total internal reflection
         w := 0;
      fi:
 end:

#9. Photon angles
  angles := proc()
    global alphai, alphat, alphac, ni, nt;
    # alphai, t, c used for transmit out or reflect 
    # alphat, r used for angularly resolved diffuse

      ni := n[layer]:

      if uz > 0 then
         nt := n[layer+1]:  # If heading down, transmitted layer is next one
      elif uz < 0 then
         nt := n[layer-1]:  # If heading up, transmitted layer is previous one
      fi:

      alphai := arccos(abs(uz));              # alphai, incident angle
      alphat := arcsin(ni/nt*sin(alphai));    # alphat, transmitted angle
      alphac := arcsin(nt/ni);                # alphac = critical angle
  end:

#10. Calculate the Fresnel Reflectance
  fresnel := proc()
    global R;
      if evalf(alphai) > evalf(alphac) then     # total reflection
          R := 1;
      else                                      # total transmission
          if alphai <> 0 then
               R := 0.5*((sin(alphai-alphat))^2/(sin(alphai+alphat))^2+(tan(alphai-alphat))^2/(tan(alphai+alphat))^2);
          else                                     # normal incidence
               R := (ni-nt)^2/(ni+nt)^2:
          fi;
      fi;
  end:

#11. Angular/radially ly Resolved transmittance
  transmit:=proc()
 global  RTdv, ATdv, TTA, TTR, TT;
    local di,k,ri,i,j;
for i from 0 to Na-1 do 
ATdv[i]:=0;od;
for j from 0 to Nr-1 do
RTdv[j]:=0;od;
for k to nops(TT) do
di:=checkA(TTA[k]); ATdv[di] := ATdv[di] + TT[k]:
if TTR[k] < dr*(Nr) then ri:=checkR(TTR[k]);
             RTdv[ri] := RTdv[ri] + TT[k]: fi;
  od;
   end:
#12. Angular/radially Resolved Diffuse reflectance
reflect:=proc()
 global  RRdv, ARdv, ATT, RTT, RR;
    local di,k,ri,i,j;
for i from 0 to Na-1 do 
ARdv[i]:=0;od;
for j from 0 to Nr-1 do
RRdv[j]:=0;od;
for k to nops(RR) do
di:=checkA(ATT[k]); ARdv[di] := ARdv[di] + RR[k]:
if RTT[k] < dr*(Nr)  then  ri:=checkR(RTT[k]);
            RRdv[ri] := RRdv[ri] + RR[k]: fi;
  od;
   end:
checkA:=proc(alpha)
local di,difference,dda;
dda:=evalf(da);
for di from 0 to Na-1 do      
      difference := abs(alpha - evalf(ATd[di]));              
      if difference < dda then
       RETURN(di);fi;
od;end:
checkR:=proc(radius)
local ri,difference,ddr;
ddr:=evalf(dr);
for ri from 0 to Nr-1 do    
      difference := abs(radius - evalf(RTd[ri]));              
      if difference < ddr then
       RETURN(ri);  fi;
od;end:




#13. Interact (update step and absorb)
  interact := proc()
  global x, y, z, sleft, dw, A, w:

  # Update position
      x := x + ux*s;
      y := y + uy*s;                   # New Cartesian Coordinates of photon
      z := z + uz*s;