montecarlomaplet.mw

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

 </Font>    <Font size="9">  RETURN(di);fi;
od;end:</Font>
checkR:=proc(radius)<Font size="9">
local ri,difference,ddr;
ddr:=evalf(dr);
for ri from 0 to Nr-1 do    
      difference := abs(radius - evalf(RTd[ri]));              
      if difference &lt; ddr then
 </Font>    <Font size="9">  RETURN(ri);  fi;
od;end:</Font></Font><Font size="9">


<Font foreground="[255,51,102]">
</Font>
</Font><Font size="20" foreground="[0,0,0]">#13. Interact (update step and absorb)</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">interact := proc()</Font><Font size="9">
  global x, y, z, sleft, dw, A, w:

  <Font foreground="[0,0,0]"># Update position</Font>
      x := x + ux*s;
      y := y + uy*s;                   <Font italic="true" foreground="[0,0,0]"># New Cartesian Coordinates of photon</Font>
      z := z + uz*s;
      sleft := 0:                      <Font italic="true" foreground="[0,0,0]"># No remaining step size</Font>
  <Font foreground="[0,0,0]"># Absorb</Font>
      dw := evalf(w*mua[layer]/mut[layer]);    <Font italic="true" foreground="[0,0,0]"># dw - photon weight absorbed by interaction  </Font>                      
      A := A + dw;            <Font italic="true" foreground="[0,0,0]"># Deposited weight   </Font>           
      w := w - dw;            <Font italic="true" foreground="[0,0,0]"># New weight is original minus</Font><Font italic="true"> <Font foreground="[0,0,0]">absorbed weight</Font></Font>
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#14. Depth Resolved absorption</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">DRabsorb := proc()</Font><Font size="9">
  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 &lt; smallz then
        smallz := difference;
        smallzi := zi;
      fi;
     od:
     DRAv[smallzi] := DRAv[smallzi] + dw:        #  print(&quot;DRabsorb&quot;,smallzi,DRAv[smallzi]);
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#15. Scatter and update directions</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">scatter := proc()</Font><Font size="9">
    local theta, psi, ux1, uy1, uz1:
    global ux, uy, uz:

    <Font foreground="[0,0,0]"># Scatter -- new directions</Font>
      if (g[layer] &lt;&gt; 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)));
                                                                                  <Font italic="true"> <Font foreground="[0,0,0]"># transmitted angle</Font></Font>
         psi := evalf(2*Pi*xi());                                                  <Font italic="true" foreground="[0,0,0]"># azimuthal angle [rad]</Font>
      else
         theta := evalf(arccos(2*xi()-1));
         psi := evalf(2*Pi*xi());
      fi;

      if abs(uz) &gt; 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;                            <Font italic="true" foreground="[0,0,0]"># New directions</Font>
      uz := uz1;
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#16. Photon roulette</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">roulette := proc()</Font><Font size="9">
  global w, wmin, m;
    if w &lt; wmin then
       if xi() &lt;= 1/m then
          w := m*w;
       else
          w := 0;
       fi;
    fi;
  <Font foreground="[153,0,153]">end:</Font>

</Font><Font size="20" foreground="[0,0,0]">#17. Load tissue properties</Font><Font size="9">
</Font><Font size="18" foreground="[153,0,153]">tissues := proc()</Font><Font size="9">
  global layers, n, mua, mus, mut, g, d, z0, z1:
  local i, j:
                      <Font italic="true" foreground="[0,0,0]"># ambient indices of refraction</Font>
  n[0] := 1: # upper ambience
  n[layers+1] := 1: # lower ambience
  d[0] := 0:<Font background="[255,255,102]" opaque="true">

</Font>                          <Font italic="true" foreground="[0,0,0]"># Find top &amp; bottom boundaries of each layer w/rt top of tissue</Font>
  for i from 0 to layers do 
      if i &lt;&gt; 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:

</Font><Font size="20" foreground="[0,0,0]">#18. New step size when crossing layers</Font><Font size="9">
  </Font><Font size="18" foreground="[153,0,153]">newstep := proc()</Font><Font size="9">
  global s;
    if uz &gt; 0 then
       s := s*mut[layer]/mut[layer+1]:
    elif uz &lt; 0 then
       s := s*mut[layer]/mut[layer-1]:
    fi:
  <Font foreground="[153,0,153]">end:
</Font></Font></Text-field>
</Input>
</Group>
</Section>
<Section collapsed="false" MultipleChoiceAnswerIndex="-1" MultipleChoiceRandomizeChoices="false" TrueFalseAnswerIndex="-1" EssayAnswerRows="5" EssayAnswerColumns="60"><Title>
<Text-field style="Heading 1" layout="Heading 1">3. Procedures to plot graphs </Text-field></Title>
<Group labelreference="L113" drawlabel="true">
<Input>
<Text-field prompt="&gt; " style="Maple Input" layout="Normal"><Equation executable="true" style="2D Input" input-equation="Typesetting:-mrow(Typesetting:-mi(&quot;&quot;))">LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnLyUnb3BhcXVlR0Y3LyUrZXhlY3V0YWJsZUdGOi8lKXJlYWRvbmx5R0Y3LyUpY29tcG9zZWRHRjcvJSpjb252ZXJ0ZWRHRjcvJStpbXNlbGVjdGVkR0Y3LyUscGxhY2Vob2xkZXJHRjcvJTBmb250X3N0eWxlX25hbWVHUSkyRH5JbnB1dEYnLyUqbWF0aGNvbG9yR0ZDLyUvbWF0aGJhY2tncm91bmRHRkYvJStmb250ZmFtaWx5R0YxLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy8lKW1hdGhzaXplR0Y0</Equation><Font size="9" foreground="[0,0,0]">## 1. Depth Resolved Absorption</Font><Font size="9">
<Font background="[255,255,51]" opaque="true" foreground="[153,0,153]">absorbgraph</Font><Font foreground="[153,0,153]"> := proc(photons)</Font>
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=[&quot;depth(cm)&quot;, Absorption], axes=boxed):
ppa:=[seq([DRA[n],temp[n]],n=0..Nz-1)];
save ppa,`Abs_graph_points.m`:
RETURN(display[plots](PP));
<Font foreground="[153,0,153]">end:</Font>
<Font foreground="[0,0,0]">## 2. Radially Resolved Reflectance</Font>
<Font background="[255,255,51]" opaque="true" foreground="[153,0,153]">radRgraph</Font><Font foreground="[153,0,153]"> := proc(photons)</Font>
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=[&quot;radius (cm)&quot;, Reflectance], axes=boxed):
pprR:=[seq([RRd[n],temp[n]],n=1..Nr-1)];
save pprR,`Ref_rad_graph_points.m`:
RETURN(display[plots](PP));
<Font foreground="[153,0,153]">end:</Font>
<Font foreground="[0,0,0]">## 3. Angularly Resolved Reflectance</Font>
<Font background="[255,255,51]" opaque="true" foreground="[153,0,153]">angRgraph</Font><Font foreground="[153,0,153]"> := proc(photons)</Font>
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);               <Font italic="true" foreground="[0,0,0]">#solid angle</Font>
    temp[i] := evalf(ARdv[i]/(photons*sa));
      od:
  PP:= pointplot({seq([ARd[n],temp[n]],n=1..Na-1)}, labels=[&quot;Angle (rad)&quot;, Reflectance], axes=boxed):
ppaR:=[seq([ARd[n],temp[n]],n=1..Na-1)];
save ppaR,`Ref_ang_graph_points.m`:
RETURN(display[plots](PP));
<Font foreground="[153,0,153]">end:</Font><Font foreground="[0,0,0]">

</Font>
<Font foreground="[0,0,0]">## 4. Angularly Resolved Transmittance</Font>
<Font background="[255,255,51]" opaque="true" foreground="[153,0,153]">angTgraph</Font><Font foreground="[153,0,153]"> := proc(photons)</Font>
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);           <Font italic="true" foreground="[0,0,0]">#solid angle</Font>
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=[&quot;Angle (rad)&quot;, Transmittance], axes=boxed,colour=black):
#PPa:= pointplot({seq([ATd[n],tempa[n]],n=1..Na-1)}, labels=[&quot;Angle (rad)&quot;, 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));
<Font foreground="[153,0,153]">end:</Font>

<Font foreground="[0,0,0]">## 5. Radially Resolved Transmittance</Font>
<Font background="[255,255,51]" opaque="true" foreground="[153,0,153]">radTgraph</Font><Font foreground="[153,0,153]"> := proc(photons)</Font>
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(&quot;radT&quot;, evalf(area));fi;
 tempa[i] := evalf(RTdv[i]/photons);
  temp[i] := evalf(RTdv[i]/(photons*area)); 
      od: