On the use of IsotropicSource in MCCL

[Karthik Sasihithlu]

Hi,

I recently started using the MCCL code to model light propation in multilayered structures.

In particular, I am trying to model the propagation of radiation from an isotropic source embedded in a multilayered structure, and monitor how the reflectance and transmittance gets modified for different configurations. 

However, I notice that when I add up the output from  RDiffuse,  TDiffuse and Absorption detectors, they are not adding up to 1.  I am not able to figure out what I am entering wrong.  Can you please help me out ?

For, the same configuration if I use a DirectionalPoint source with incidence from the top,  then RDiffuse + RSpecular+ TDiffuse+ Absorption does indeed add up to 1.  But not so for isotropic source.

I have attached the infile I am using.  I am monitoring the total reflectance, total transmittance and total absorption in the structure.

Thanks!

Karthik

[Carole Hayakawa]

Hi Karthik,

Thanks for trying our software!  I see from your infile that you have a slab that is 3mm thick with optical properties (OPs):

```

          "Mua": 0.0,
          "Mus": 1E-10,
          "G": 1.0,
          "N": 2,

          "Musp": 1E-10

```

This is essentially air with a higher refractive index.   So your isotropic source photons are just going to transport in straight rays until they exit the top or bottom of the slab (the slab is infinite laterally).   Here are a few things that I did:

1) set TallySecondMoment = true for all detectors.  This will allow you to determine the error in the Mean values (load_results_script) will give you access to the standard deviation of the estimates. I have modified your infile and load_results_script.m to show you how you have access to the error.  The standard deviation divided by the mean should be < 5% if possible.  If it isn't, then increase the number of photons launched.

2) In "Options" section of infile, "AbsorptionWeightingType" is set to "Discrete".  This type of tally accounts for absorption by deweighting the photon at each collision and tallying to ATotal at that collision.  Since the slab has OPs of air no collisions will occur within the slab, so ATotal will be 0.  Modify "Discrete" to "Continuous".  This will account for absorption by deweighting continuously along the photon path length and will result in an non-zero ATotal.  I just tried this and ATotal is still 0.  This is a bug!

Let me debug this and post again.

Best,

Carole

[Carole Hayakawa]

Hi Karthik,

I stepped through the code and ATotal being 0 makes sense to me now.  Its because you have Mua=0 in the slab.  (I should have realized that earlier).  I get the following results from load_results_script.m using 10^5 photons:

Total absorption captured by ATotal detector: 0
Standard Deviation captured by ATotal detector: 0
Total reflectance captured by RDiffuse detector: 0.06661
Standard Deviation captured by RDiffuse detector: 0.0007885
Total reflectance captured by RSpecular detector: 0
Total transmittance captured by TDiffuse detector: 0.06808
Standard Deviation captured by TDiffuse detector: 0.00079652

My relative error [SD/Mean] 1.2% for RDiffuse and 1.2% for TDiffuse, however RDiffuse + TDiffuse = 0.1347, not 1.

I forgot to attached my edited version of your infile and also my edits to load_results_script.m to see the error in the estimates.

In Options, I set "TrackStatistics" to true.   With this set to true, a statistics.txt file is created after running the simulation.  In that file it shows me that 865531 photons were killed because their path length is greater than the maximum path length allowable 60,000mm.  

```

  "NumberOfPhotonsOutTopOfTissue": 6661,
  "NumberOfPhotonsOutBottomOfTissue": 6808,
  "NumberOfPhotonsAbsorbed": 0,
  "NumberOfPhotonsSpecularReflected": 0,
  "NumberOfPhotonsKilledOverMaximumPathLength": 86531,
  "NumberOfPhotonsKilledOverMaximumCollisions": 0,
  "NumberOfPhotonsKilledByRussianRoulette": 0

```

This accounts for the loss.  This amounts to a weight loss of 0.8651, so now 0.8651+0.06661+0.06808=0.9998.  These are the photons that are emanating from the isotropic source moving laterally until the maximum path length allowable.

Let me know if you have any further questions.

Best,

Carole

[Karthik Sasihithlu]

Thanks a lot!  Makes sense now.  

So these are being lost due to 'waveguiding' effect.  These are trapped in the layers due to total internal reflection and never reach the reflection/transmission detectors. 

Karthik

[Karthik Sasihithlu]

Dear Carole,

You mentioned that the maximum path length allowable 60,000mm.  Is there a quick way to custom set this to a lower value, or would I have to make changes in the source code for this ?  

Thanks!

Karthik

[Carole Hayakawa]

Hi Karthik,

You would have to make changes in the source code.

I took a stab at what might work for you.  We have a tissue type that surrounds the tissue layer with a cylinder.  If you make that an infinite absorber, we have a detector that will tally all photons that enter the absorber (see the attached infile).  I have defined a *bounding* cylinder the height of the 3mm slab with radius=5mm.  I specified that it is an infinite absorber and have the N of the slab so that any photons that hit the cylinder, directly enter and get absorbed.  I added a detector ATotalBoundingVolume to the infile.  There is a little nit,  I modified the slab layer Mus and Musp to be 1E-9 instead of 1E-10.  We have code that checks that the "tissue" layer is not air and checks if these values are <=1E-10.  I might remove that check after seeing your post here and what you are trying to do.

When I run this infile and run load_results_script.m (I added standard deviation output for added detector and attached below), I get the following output:

Total absorption captured by ATotal detector: 0
Standard Deviation captured by ATotal detector: 0

Total reflectance captured by RDiffuse detector: 0.06712
Standard Deviation captured by RDiffuse detector: 0.0007913

Total reflectance captured by RSpecular detector: 0

Total transmittance captured by TDiffuse detector: 0.06626
Standard Deviation captured by TDiffuse detector: 0.00078657
Total absorption captured by ATotalBoundingVolume detector: 0.86662
Standard Deviation captured by ATotalBoundingVolume detector: 0.0010751

This detector captures that lateral lost light.  Adding 0.06712+0.06626+0.86662 gives back 1.

Let me know if this is helpful or you had something else in mind.

Best,

Carole