montecarlomaplet.maplet

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

      sleft := 0:                      # No remaining step size
  # Absorb
      dw := evalf(w*mua[layer]/mut[layer]);    # dw - photon weight absorbed by interaction                        
      A := A + dw;            # Deposited weight              
      w := w - dw;            # New weight is original minus absorbed weight
  end:

#14. Depth Resolved absorption
  DRabsorb := proc()
  global smallz, smallzi, DRAv:
  local zi, difference:

     smallz := 99:
     for zi from 0 to Nz-1 do
     difference := evalf(abs(z-DRA[zi]));
      if difference < smallz then
        smallz := difference;
        smallzi := zi;
      fi;
     od:
     DRAv[smallzi] := DRAv[smallzi] + dw:        #  print("DRabsorb",smallzi,DRAv[smallzi]);
  end:

#15. Scatter and update directions
  scatter := proc()
    local theta, psi, ux1, uy1, uz1:
    global ux, uy, uz:

    # Scatter -- new directions
      if (g[layer] <> 0) then
         #cos(theta) := 1/(2*g)*(1+g^2-((1-g^2)/(1-g+2*g*xi()))^2);
         theta := evalf(arccos(1/(2*g[layer])*(1+g[layer]*g[layer]-((1-g[layer]*g[layer])/(1-g[layer]+2*g[layer]*xi()))^2)));
                                                                                   # transmitted angle
         psi := evalf(2*Pi*xi());                                                  # azimuthal angle [rad]
      else
         theta := evalf(arccos(2*xi()-1));
         psi := evalf(2*Pi*xi());
      fi;

      if abs(uz) > 0.99999 then
         ux1 := sin(theta)*cos(psi);
         uy1 := sin(theta)*sin(psi);
         uz1 := signum(uz)*cos(theta);
      else
         ux1 := sin(theta)*(ux*uz*cos(psi)-uy*sin(psi))/sqrt(1-uz*uz)+ux*cos(theta);
         uy1 := sin(theta)*(uy*uz*cos(psi)+ux*sin(psi))/sqrt(1-uz*uz)+uy*cos(theta);
         uz1 := -sin(theta)*cos(psi)*sqrt(1-uz*uz)+uz*cos(theta);
      fi;

      ux := ux1;
      uy := uy1;                            # New directions
      uz := uz1;
  end:

#16. Photon roulette
  roulette := proc()
  global w, wmin, m;
    if w < wmin then
       if xi() <= 1/m then
          w := m*w;
       else
          w := 0;
       fi;
    fi;
  end:

#17. Load tissue properties
tissues := proc()
  global layers, n, mua, mus, mut, g, d, z0, z1:
  local i, j:
                      # ambient indices of refraction
  n[0] := 1: # upper ambience
  n[layers+1] := 1: # lower ambience
  d[0] := 0:

                          # Find top & bottom boundaries of each layer w/rt top of tissue
  for i from 0 to layers do 
      if i <> 0 then
         mut[i] := mua[i] + mus[i]:
      fi:
      for j from 1 to i do
          z0[i] := z0[i-1] + d[j-1]:
          z1[i] := z1[i] + d[j]:
      od:
  od:
end:

#18. New step size when crossing layers
  newstep := proc()
  global s;
    if uz > 0 then
       s := s*mut[layer]/mut[layer+1]:
    elif uz < 0 then
       s := s*mut[layer]/mut[layer-1]:
    fi:
  end:

# 3. Procedures to plot graphs 
## 1. Depth Resolved Absorption
absorbgraph := proc(photons)
local i, temp,PP,ppa:
temp := array(0..Nz-1):
  for i from 0 to Nz-1 do
    temp[i] := evalf(DRAv[i]/(photons*dz));
      od:
PP:= pointplot({seq([DRA[n],temp[n]],n=0..Nz-1)}, labels=["depth(cm)", Absorption], axes=boxed):
ppa:=[seq([DRA[n],temp[n]],n=0..Nz-1)];
save ppa,`Abs_graph_points.m`:
RETURN(display[plots](PP));
end:
## 2. Radially Resolved Reflectance
radRgraph := proc(photons)
local i, temp,PP,pprR,area:
temp := array(0..Nr-1):
  for i from 0 to Nr-1 do
area:=2*Pi*(i+0.5)*dr*dr;   
temp[i] := evalf(RRdv[i]/(photons*area));
      od:
PP:= pointplot({seq([RRd[n],temp[n]],n=1..Nr-1)}, labels=["radius (cm)", Reflectance], axes=boxed):
pprR:=[seq([RRd[n],temp[n]],n=1..Nr-1)];
save pprR,`Ref_rad_graph_points.m`:
RETURN(display[plots](PP));
end:
## 3. Angularly Resolved Reflectance
angRgraph := proc(photons)
local i, temp,PP,ppaR,sa:
temp := array(0..Na-1):
  for i from 0 to Na-1 do
sa:=4*Pi*sin((i+0.5)*da)*sin(da/2);               #solid angle
    temp[i] := evalf(ARdv[i]/(photons*sa));
      od:
  PP:= pointplot({seq([ARd[n],temp[n]],n=1..Na-1)}, labels=["Angle (rad)", Reflectance], axes=boxed):
ppaR:=[seq([ARd[n],temp[n]],n=1..Na-1)];
save ppaR,`Ref_ang_graph_points.m`:
RETURN(display[plots](PP));
end:


## 4. Angularly Resolved Transmittance
angTgraph := proc(photons)
local i, temp,tempa,PP,PPa,ppaT,ppaTa,sa:
temp := array(0..Na-1):tempa := array(0..Na-1):
  for i from 0 to Na-1 do
sa:=4*Pi*sin((i+0.5)*da)*sin(da/2);           #solid angle
tempa[i] := evalf(ATdv[i]/photons);
temp[i] := evalf(ATdv[i]/(photons*sa));
          od:
PP:= pointplot({seq([ATd[n],temp[n]],n=1..Na-1)}, labels=["Angle (rad)", Transmittance], axes=boxed,colour=black):
#PPa:= pointplot({seq([ATd[n],tempa[n]],n=1..Na-1)}, labels=["Angle (rad)", Transmittance], axes=boxed,colour=red):
ppaT:=[seq([ATd[n],temp[n]],n=1..Na-1)];
#ppaTa:=[seq([ATd[n],tempa[n]],n=1..Na-1)];
save ppaT, `Tra_ang_graph_points.m`:
RETURN(display[plots](PP));
end:

## 5. Radially Resolved Transmittance
radTgraph := proc(photons)
local i, temp,tempa,PP,PPa,pprT,pprTa,area:
temp := array(0..Nr-1):tempa := array(0..Nr-1):
  for i from 0 to Nr-1 do
area:=2*Pi*(i+0.5)*dr*dr;                  #####  if i=1 then print("radT", evalf(area));fi;
 tempa[i] := evalf(RTdv[i]/photons);
  temp[i] := evalf(RTdv[i]/(photons*area)); 
      od:
PP:= pointplot({seq([RTd[n],temp[n]],n=1..Nr-1)}, labels=["radius (cm)", Transmittance], axes=boxed,colour=black):
#PPa:= pointplot({seq([RTd[n],tempa[n]],n=1..Nr-1)}, labels=["radius (cm)", Transmittance], axes=boxed,colour=red):
pprT:=[seq([RTd[n],temp[n]],n=1..Nr-1)];
#pprTa:=[seq([RTd[n],tempa[n]],n=1..Nr-1)];
save pprT,`Tra_rad_graph_points.m`:
RETURN(display[plots](PP));
end:


# II - The Maplet
# 1. Construction of the Maplet 
use Maplets[Elements] in
with(Maplets[Elements]):
    MonteCarloMaplet:= Maplet('onstartup' = RunWindow( 'W0' ),
 Font['F00']( 'family' = "times",bold,'size'=12 ),
               Font['F0']( 'family' = "arial",bold,'size'=12 ),
               Font['F1']( 'family' = "arial",bold,'size'=14 ),
               Font['F2']( 'family' = "times",bold, 'size'=20 ),
               Font['F3']( 'family' = "times",bold, 'size'=14 ),

Window['W0']('title'= "MAIN WINDOW",
BoxRow(
BoxColumn( #### BoxColumn 1 starts
'hscroll'='always',
'vscroll'='always','inset'=0,'spacing'=15,'valign'='top',
       Label("INPUT PANEL      ",'font'=Font(Arial,bold,16)),
[Button("Explanation of INPUT ",'foreground'=black,'background'=white, 'font'='F1', RunWindow('Win'))],
[BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Number of photons  " ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left',"Enter a number >500  "), 
           GridCell('halign'='left',
              TextField['TF_2']('width'=6,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"   ")
))) #  end gridrow ; gridlayout ; BoxRow
,
BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Entry Angle  " ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left',"Enter a number between 0 and 90  "), 
           GridCell('halign'='left',
              TextField['TF_0']('width'=6,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"   ")
)))] #  end gridrow ; gridlayout ; BoxRow
,
[BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Number of tissue layers " ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left',"Enter an integer value  "), 
           GridCell('halign'='left',
              TextField[TF_1]('width'=3,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"    ")
))) #  end gridrow ; gridlayout ; boxrow
,

BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Layers' refraction indices" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left'," [n1, n2, ...]     "), 
           GridCell('halign'='left',
              TextField['TF_3']('width'=5,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"   ")
)))] #  end gridrow ; gridlayout ; boxrow
,

[BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Layers' scattering coefficients" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left'," [cs1, cs2, ...]      "), 
           GridCell('halign'='left',
              TextField['TF_4']('width'=5,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"  1/cm ")
))) #  end gridrow ; gridlayout ; BoxRow
,
BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Layers' absorption coefficients " ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left'," [ca1, ca2, ...]     "), 
           GridCell('halign'='left',
              TextField['TF_5']('width'=5,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"  1/cm ")
)))] #  end gridrow ; gridlayout ; BoxRow
,
[BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Layers' anisotropy coefficients" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left',"[g1, g2, ...]     "), 
           GridCell('halign'='left',
              TextField['TF_6']('width'=5,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"   ")
))) #  end gridrow ; gridlayout ; BoxRow
,
BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Layers' depths" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left'," [d1, d2, ...]     "), 
           GridCell('halign'='left',
              TextField['TF_7']('width'=5,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"  cm ")
)))] #  end gridrow ; gridlayout ; BoxRow
),
[BoxColumn('hscroll'='always',
'vscroll'='always','inset'=0,'spacing'=15,'valign'='top',
       Label("   SELECT GEOMETRY TYPE     ",'font'=Font(Arial,bold,16)),

BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Cylindrical Layers" ,'font'='F1',     
GridRow('halign'='left',
           GridCell('halign'='left',"Input maximum radius  "), 
           GridCell('halign'='left',
              TextField['TF_r']('width'=3,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"  cm ")
))), #  end gridrow ; gridlayout ; BoxRow

CheckBox['ChB1']('caption'="For cylindrical layers check here", 'value' = 'false' ),   
#[Button("Proceed ",'background'=white,'foreground'=black, 'font'='F1', RunWindow('W1'))],
BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Prismatic Layers" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left',"Input maximum x,y dimensions [x,y]    "), 
           GridCell('halign'='left',
              TextField['TF_xy']('width'=6,'background'=grey, 'editable'='true', 'font'='F1')),  
GridCell('halign'='left',"  cm ")
))), #  end gridrow ; gridlayout ; BoxRow

CheckBox['ChB2']('caption'="For prismatic layers check here", 'value' = 'false' ),   

BoxRow(
GridLayout('border'='true', 'inset'=0, 
          'caption'="Infinitely Long Layers" ,'font'='F1',        
GridRow('halign'='left',
           GridCell('halign'='left'," Infinite radius assumed   ")#, 
           
))
) ,  #  end gridrow ; gridlayout ; BoxRow
CheckBox['ChB3']('caption'="For inifnite layers check here", 'value' = 'false' ), 

MathMLViewer['ML']('foreground'=black,'background'=grey,'height'=2,
       'width'=1), 
Button("CLEAR Run Status BOX", 
Action(SetOption('ML' = ""))),
CheckBox['ChB4']('caption'="If only numerical results are required, check here", 'value' = 'false' ),   
Button("Start Simulation",'foreground'=black, 'background'=white,'font'='F00',Evaluate('ML' = 'SimulationMC()')) 
 , 
[Button("Explanation of OUTPUT",'background'=white, 'foreground'=black,'font'='F1',RunWindow('Wout'))],