montecarlomaplet.mw

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

<Output>
<Text-field style="Warning" layout="Warning">Warning, the name changecoords has been redefined</Text-field>
</Output>
</Group>
</Section>
<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title>
<Text-field style="Heading 1" layout="Heading 1">2. Procedures for the simulation</Text-field></Title>
<Group labelreference="L112" drawlabel="true">
<Input>
<Text-field prompt="&gt; " style="Maple Input" layout="Normal"><Font size="9">
   </Font><Font size="20" foreground="[0,0,0]">#1. Load Simulation properties </Font><Font size="9">
   </Font><Font size="18" foreground="[153,0,153]">simprop := proc()</Font><Font size="9">
   local i:
   global ARdv, ARd, ATdv, ATd, RRdv, RRd, RTdv, RTd, DRAv, DRA, Na, Nr, dr,  Nz, dz, wmin, m, da:
<Font foreground="[0,0,0]">   # Parameters for Angular Resolution. Can be changed in INPUT</Font>

<Font background="[255,255,204]" opaque="true">      Na := 30:                  <Font italic="true" foreground="[0,0,0]"># INPUT number of grid elements for exit angle, alphat</Font>
      da := Pi/(2*Na):     <Font italic="true" foreground="[0,0,0]"># delta_alpha, angle resolution    </Font>                   </Font>
      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;             <Font italic="true" foreground="[0,0,0]"># Rd values initialized to 0</Font>
        ARd[i] := (i+0.5)*da; #(1/2)*(i*Pi/(2*Na)+(i+1)*Pi/(2*Na));
          ATdv[i] := 0;           <Font italic="true" foreground="[0,0,0]"># Td values initialized to 0</Font>
          ATd[i] := (i+0.5)*da;#(1/2)*(i*Pi/(2*Na)+(i+1)*Pi/(2*Na));
      od:

   <Font foreground="[0,0,0]"># Parameters for Radial Resolution. Can be changed in INPUT</Font>
<Font background="[255,255,204]" opaque="true">      Nr := 200:                    <Font italic="true" foreground="[0,0,0]"># INPUT number of radial grid elements</Font>
      dr := 0.05:            <Font italic="true" foreground="[0,0,0]"># INPUT radial resolution [cm]  </Font><Font foreground="[0,0,0]">   </Font>   </Font>
      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;              <Font italic="true" foreground="[0,0,0]"># Rd values initialized to 0</Font>
        RRd[i] := (i+0.5)*dr;
        RTdv[i] := 0;              <Font italic="true" foreground="[0,0,0]"># Td values initialized to 0</Font>
        RTd[i] := (i+0.5)*dr;
      od:
   <Font foreground="[0,0,0]"># Parameters for Depth Absorption Resolution</Font>
   <Font background="[102,255,102]" opaque="true">
</Font><Font background="[255,255,204]" opaque="true">      Nz := 20;    <Font italic="true" foreground="[0,0,0]"># *INPUT* number of depth grid elements</Font><Font italic="true"> </Font> 
      dz := totald/Nz;     <Font italic="true" foreground="[0,0,0]"># *INPUT* depth resolution where totald is the sum of the thickness of all layers</Font>        </Font></Font><Font background="[255,255,204]" opaque="true" size="10">  </Font><Font size="9">
      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:
   <Font foreground="[0,0,0]"># Photon Properties</Font>
      <Font background="[255,255,204]" opaque="true">wmin := 0.0001: m := 10:</Font>
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#2.  Load random number generator</Font><Font size="9">
<Font opaque="true">  </Font></Font><Font opaque="true" size="18" foreground="[153,0,153]">xi := proc()</Font><Font opaque="true" size="9">
    evalf(rand()/999999999999);
  <Font foreground="[153,0,153]">end:</Font>
</Font><Font size="9">
</Font><Font size="20" foreground="[0,0,0]">#3.  Load initial conditions</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">initialcond := proc()</Font><Font size="9">
      global s, sleft, w,  x, y, z, Rsp: <Font italic="true" foreground="[0,0,0]">
# step size, remaining step, photon weight, x,y,z) co-ordinates, specular reflectance
</Font>
         x := 0: y := 0: z := 0:       <Font italic="true" foreground="[0,0,0]"> </Font>    
             s := 0: 
         sleft := 0: 
             w := 1: 
         Rsp := 0: 
          
  <Font foreground="[153,0,153]">end:</Font>
</Font><Font size="20" foreground="[0,0,0]">#4. Initial directions for launch angle </Font><Font size="9"> 
  </Font><Font size="18" foreground="[153,0,153]">initialangle := proc(langle)</Font><Font size="9">
  local temp;
  global ux, uy, uz:  <Font italic="true" foreground="[0,0,0]"># (ux, uy, uz) angle cossines</Font>
    temp := 90 - langle: <Font italic="true" foreground="[0,0,0]"># angle relative to surface's normal</Font>
    temp := evalf(temp*Pi/180): <Font italic="true" foreground="[0,0,0]"># convert to radians</Font>
    uz := cos(temp):
    ux := sqrt((1-uz*uz)/2):
    uy := ux:	<Font italic="true" foreground="[0,0,0]">#  ux=uy since ux, uy are assumed arbitrary</Font></Font>
<Font size="9">end:

</Font><Font opaque="true" size="20" foreground="[0,0,0]">#5. Specular Reflection</Font><Font opaque="true" size="9" foreground="[255,255,255]">
  </Font><Font opaque="true" size="18" foreground="[153,0,153]">specular := proc()</Font><Font opaque="true" size="9" foreground="[255,255,255]">
  </Font><Font opaque="true" size="9">  global Rsp, TRsp, w: 
    Rsp := ((n[0]-n[1])^2)/((n[0]+n[1])^2): 
    TRsp := TRsp + Rsp: 
    w := 1-Rsp:<Font foreground="[255,255,255]">
 </Font><Font foreground="[153,0,153]"> end:</Font><Font foreground="[255,255,255]">
</Font></Font><Font size="9">
</Font><Font size="20" foreground="[0,0,0]">#6. Step size calculation</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">stepsize := proc()</Font><Font size="9">
    global s, sleft:
    if sleft = 0 then
        s := -ln(xi())/mut[layer]:
    else
        s := sleft/mut[layer]:
        sleft := 0:
    fi:
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font opaque="true" size="20" foreground="[0,0,0]">#7. Find distance calculation</Font><Font opaque="true" size="9">
  </Font><Font opaque="true" size="18" foreground="[153,0,153]">finddistance := proc()</Font><Font opaque="true" size="9">
    global db, r:
      if uz &gt; 0 then
         db := (z1[layer]-z)/uz:      <Font italic="true" foreground="[0,0,0]"># from lower layer</Font>
      elif uz &lt; 0 then
         db := (z0[layer]-z)/uz:      <Font italic="true" foreground="[0,0,0]"># from upper layer</Font>
      fi:
        r := sqrt(x*x+y*y):           <Font italic="true" foreground="[0,0,0]"># radial position</Font>
  <Font foreground="[153,0,153]">end:</Font></Font><Font size="9">

</Font><Font size="20" foreground="[0,0,0]">#8. Move photon to boundary</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">movetoboundary := proc()</Font><Font size="9">
    global x, y, z, w, layer;
      if uz &gt; 0 then  	<Font italic="true" foreground="[0,0,0]"># going down</Font>
         x := x + ux*db;
         y := y + uy*db;
         z := z1[layer];
      elif uz &lt; 0 then	<Font italic="true" foreground="[0,0,0]"># going up</Font>
         x := x + ux*db;
         y := y + uy*db;
         z := z0[layer];
      else		  <Font italic="true" foreground="[0,0,0]"># total internal reflection</Font>
         w := 0;
      fi:
 end:

</Font><Font size="20" foreground="[0,0,0]">#9. Photon angles</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">angles := proc()</Font><Font size="9">
    global alphai, alphat, alphac, ni, nt;
    <Font italic="true" foreground="[0,0,0]"># alphai, t, c used for transmit out or reflect 
    # alphat, r used for angularly resolved diffuse</Font>

      ni := n[layer]:

      if uz &gt; 0 then
         nt := n[layer+1]:  <Font italic="true" foreground="[0,0,0]"># If heading down, transmitted layer is next one</Font>
      elif uz &lt; 0 then
         nt := n[layer-1]:  <Font italic="true" foreground="[0,0,0]"># If heading up, transmitted layer is previous one</Font><Font background="[255,255,102]" opaque="true">
</Font>      fi:

      alphai := arccos(abs(uz));              <Font italic="true" foreground="[0,0,0]"># alphai, incident angle</Font>
      alphat := arcsin(ni/nt*sin(alphai));    <Font italic="true" foreground="[0,0,0]"># alphat, transmitted angle</Font>
      alphac := arcsin(nt/ni);                <Font italic="true" foreground="[0,0,0]"># alphac = critical angle</Font>
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#10. Calculate the Fresnel Reflectance</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">fresnel := proc()</Font><Font size="9">
    global R;
      if evalf(alphai) &gt; evalf(alphac) then     <Font italic="true" foreground="[0,0,0]"># total reflection</Font>
          R := 1;
      else                                      <Font italic="true" foreground="[0,0,0]"># total transmission</Font>
          if alphai &lt;&gt; 0 then
               R := 0.5*((sin(alphai-alphat))^2/(sin(alphai+alphat))^2+(tan(alphai-alphat))^2/(tan(alphai+alphat))^2);
          else                                     <Font italic="true" foreground="[0,0,0]"># normal incidence</Font>
               R := (ni-nt)^2/(ni+nt)^2:
          fi;
      fi;
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#11. Angular/radially ly Resolved transmittance
</Font><Equation executable="true" style="2D Input" input-equation="">NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNictSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEyRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHUS5bMjU1LDI1NSwyNTVdRigvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y7LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKTJEfklucHV0RigvJSptYXRoY29sb3JHRkQvJS9tYXRoYmFja2dyb3VuZEdGRy8lK2ZvbnRmYW1pbHlHRjIvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YoLyUpbWF0aHNpemVHRjUtSSdtc3BhY2VHRiU2Ji8lJ2hlaWdodEdRJzAuMH5leEYoLyUmd2lkdGhHUScwLjB+ZW1GKC8lJmRlcHRoR0Zjby8lKmxpbmVicmVha0dRKG5ld2xpbmVGKEYsRl5vRiw3I0Yo</Equation><Font size="18" foreground="[153,0,153]">transmit:=proc()</Font><Font size="9">
 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] &lt; dr*(Nr) then ri:=checkR(TTR[k]);
             RTdv[ri] := RTdv[ri] + TT[k]: fi;
  od;<Font background="[255,255,102]" opaque="true">
</Font><Font foreground="[255,51,102]">  </Font><Font foreground="[153,0,153]"> end:</Font><Font foreground="[255,51,102]">
</Font></Font><Font size="20" foreground="[0,0,0]">#12. Angular/radially Resolved Diffuse reflectance</Font><Font size="9" foreground="[255,51,102]">
</Font><Font size="18" foreground="[153,0,153]">reflect:=proc()</Font><Font size="9" foreground="[255,51,102]">
</Font><Font size="9"> 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] &lt; dr*(Nr)  then  ri:=checkR(RTT[k]);
            RRdv[ri] := RRdv[ri] + RR[k]: fi;
  od;<Font foreground="[255,51,102]">
  </Font><Font foreground="[153,0,153]"> end:</Font>
</Font><Font foreground="[153,0,153]">checkA:=proc(alpha)<Font size="9">
local di,difference,dda;
dda:=evalf(da);
for di from 0 to Na-1 do      
      difference := abs(alpha - evalf(ATd[di]));              
      if difference &lt; dda then