Validated refractive index mismatched scenarios using same example fr…

…om lewis's original paper on sMCRT

app/main.f90

@@ -8,7 +8,10 @@ program mcpolar
     num_args = command_argument_count()
     if(num_args == 0)then
         allocate(args(1))
-        args(1) = "scat_test.toml"
+        !args(1) = "scat_test.toml"
+        !args(1) = "validation1.toml"   !Validate hgg scattering, absorption and scattering
+        !args(1) = "validation2.toml"   !Validate refractive index mismatch is working
+        args(1) = "validation3.toml"   !Validate refractive index mismatch
     else
         allocate(args(num_args))
         do i = 1, num_args

res/scat_test.toml

@@ -1,8 +1,8 @@
 [source]
 name = "pencil"
-nphotons = 1000000
+nphotons = 10000000
 direction = "z"
-position = [0.0,0.0,-1.0]
+position = [0.0,0.0,-2.45]
 spectrum_type = "constant"
 wavelength = 500.0
 
@@ -12,34 +12,27 @@ nyg = 500
 nzg = 500
 xmax = 50.0
 ymax = 50.0
-zmax = 1.5
+zmax = 5.0
 
 [geometry]
 geom_name = "box"
-BoxDimensions = [100.0,100.0,2.0]
-boundingBox = [100.0,100.0,3.0]
-mus = [0.9]
-hgg = [0.75]
-mua = [0.1]
+BoxDimensions = [100.0,100.0,4.9]
+boundingBox = [100.0,100.0,5.0]
+position = [0.0,0.0,0.0]
+mus = [90]
+hgg = [0]
+mua = [10]
+n = [1.5]
 
 [[detectors]]
 type="circle"
-position=[0.0, 0.0, -1.001]
+position=[0.0, 0.0, -2.46]
 direction=[0.0, 0.0, -1.0]
 radius1=20
 nbins=100
 maxval = 1.0
 trackHistory=false
 
-[[detectors]]
-type="circle"
-position=[0.0, 0.0, +1.001]
-direction=[0.0, 0.0, +1.0]
-radius1=100
-nbins=100
-maxval = 1.0
-trackHistory=false
-
 [output]
 fluence = "fluence.nrrd"
 render_geometry_name = "geom_render.nrrd"

res/validation1.toml

@@ -0,0 +1,59 @@
+[source]
+name = "pencil"
+nphotons = 1000000
+direction = "z"
+position = [0.0,0.0,-0.01]
+spectrum_type = "constant"
+wavelength = 500.0
+
+[grid]
+nxg = 500
+nyg = 500
+nzg = 500
+xmax = 50.0
+ymax = 50.0
+zmax = 0.015
+
+[geometry]
+geom_name = "box"
+BoxDimensions = [100.0,100.0,0.02]
+boundingBox = [100.0,100.0,0.03]
+position = [0.0,0.0,0.0]
+mus = [90]
+hgg = [0.75]
+mua = [10]
+n = [1.0]
+
+[[detectors]]
+type="circle"
+position=[0.0, 0.0, -0.011]
+direction=[0.0, 0.0, -1.0]
+radius1=20
+nbins=100
+maxval = 1.0
+trackHistory=false
+
+[[detectors]]
+type="circle"
+position=[0.0, 0.0, 0.011]
+direction=[0.0, 0.0, 1.0]
+radius1=20
+nbins=100
+maxval = 1.0
+trackHistory=false
+
+[output]
+fluence = "fluence.nrrd"
+render_geometry_name = "geom_render.nrrd"
+render_geomerty = true
+render_source_name = "source_render.nrrd"
+render_source = true
+render_size = [200,200,200]
+overwrite = true
+
+[simulation]
+iseed = 123456789
+absorb=true
+load_checkpoint=false
+checkpoint_file="check.ckpt"
+checkpoint_every_n=10000

res/validation2.toml

@@ -0,0 +1,43 @@
+[source]
+name = "uniform"
+nphotons = 1000000
+direction = "-z"
+point1 = [ -5.0, -5.0, 2.0]
+point2 = [ 10.0, 0.0, 0.0]
+point3 = [ 0.0, 10.0, 0.0]
+spectrum_type = "constant"
+wavelength = 500.0
+
+[grid]
+nxg = 250
+nyg = 250
+nzg = 1000
+xmax = 50.0
+ymax = 50.0
+zmax = 2.0
+
+[geometry]
+geom_name = "box"
+BoxDimensions = [100.0,100.0,3.9]
+boundingBox = [100.0,100.0,4.0]
+position = [0.0,0.0,0.0]
+mus = [820]
+hgg = [0.9]
+mua = [1.8]
+n = [1.38]
+
+[output]
+fluence = "fluence.nrrd"
+render_geometry_name = "geom_render.nrrd"
+render_geomerty = true
+render_source_name = "source_render.nrrd"
+render_source = true
+render_size = [200,200,200]
+overwrite = true
+
+[simulation]
+iseed = 123456789
+absorb=true
+load_checkpoint=false
+checkpoint_file="check.ckpt"
+checkpoint_every_n=10000

res/validation3.toml

@@ -0,0 +1,43 @@
+[source]
+name = "uniform"
+nphotons = 1000000
+direction = "-z"
+point1 = [ -5.0, -5.0, 2.0]
+point2 = [ 10.0, 0.0, 0.0]
+point3 = [ 0.0, 10.0, 0.0]
+spectrum_type = "constant"
+wavelength = 500.0
+
+[grid]
+nxg = 250
+nyg = 250
+nzg = 1000
+xmax = 50.0
+ymax = 50.0
+zmax = 2.0
+
+[geometry]
+geom_name = "box"
+BoxDimensions = [100.0,100.0,3.9]
+boundingBox = [100.0,100.0,4.0]
+position = [0.0,0.0,0.0]
+mus = [210]
+hgg = [0.9]
+mua = [0.23]
+n = [1.38]
+
+[output]
+fluence = "fluence.nrrd"
+render_geometry_name = "geom_render.nrrd"
+render_geomerty = true
+render_source_name = "source_render.nrrd"
+render_source = true
+render_size = [200,200,200]
+overwrite = true
+
+[simulation]
+iseed = 123456789
+absorb=true
+load_checkpoint=false
+checkpoint_file="check.ckpt"
+checkpoint_every_n=10000

src/inttau2.f90

@@ -34,8 +34,9 @@ subroutine tauint2(grid, packet, sdfs_array, dects, history)
         type(history_stack_t), intent(inout) :: history
 
         real(kind=wp) :: tau, d_sdf, t_sdf, taurun, ds(size(sdfs_array)), dstmp(size(sdfs_array))
+        real(kind=wp) :: dsNew(size(sdfs_array))
         real(kind=wp) :: eps, dtot, old(size(sdfs_array)), new(size(sdfs_array)), n1, n2, Ri
-        integer       :: i, oldlayer, new_layer, smallStepLayer
+        integer       :: i, old_layer, new_layer, smallStepLayer, layer
         type(vector)  :: pos, dir, oldpos, N, smallStepPos
         logical       :: rflag
 
@@ -148,7 +149,7 @@ subroutine tauint2(grid, packet, sdfs_array, dects, history)
                 !print*, taurun>= tau, packet%tflag
                 exit
             end if
-            
+
             !move to the edge or till the packet has moved the full optical depth
             do while(d_sdf >= eps)
                 t_sdf = d_sdf * sdfs_array(packet%layer)%getkappa()
@@ -210,14 +211,19 @@ subroutine tauint2(grid, packet, sdfs_array, dects, history)
             !print*, "across boundary"
 
             !step a small bit into next sdf to get n2
+            ! choose a different step distance to ensure that you actually move into the next medium
+            !
+            !
+            !
             d_sdf = minval(abs(ds), dim=1) + 2.0_wp*eps
             smallStepPos = pos + d_sdf*dir
-            ds = 0._wp
+            dsNew = 0._wp
             do i = 1, size(ds)
-                ds(i) = sdfs_array(i)%evaluate(smallStepPos)
+                dsNew(i) = sdfs_array(i)%evaluate(smallStepPos)
             end do
-            packet%cnts = packet%cnts + size(ds)
-            new_layer = maxloc(ds,dim=1, mask=(ds<=0._wp))
+            packet%cnts = packet%cnts + size(dsNew)
+            new_layer = maxloc(dsNew,dim=1, mask=(dsNew<=0._wp))
+            old_layer = packet%layer
 
             if (new_layer == 0) then
                 ! There is no new layer and we are outside of the SDFs
@@ -231,46 +237,50 @@ subroutine tauint2(grid, packet, sdfs_array, dects, history)
 
             !carry out refelction/refraction
             if (n1 /= n2) then !check for fresnel reflection
-                print*, "n are different"
+                !print*, "n are different"
                 !Need to test and understand
                 !get correct sdf normal
-                if(ds(packet%layer) < 0._wp .and. ds(new_layer) < 0._wp)then
-                    oldlayer = minloc(abs([ds(packet%layer), ds(new_layer)]), dim=1)
-
-                elseif(dstmp(packet%layer) < 0._wp .and. dstmp(new_layer) < 0._wp)then
-                    oldlayer = maxloc([dstmp(packet%layer), dstmp(new_layer)], dim=1)
-
-                elseif(ds(packet%layer) > 0._wp .and. ds(new_layer) < 0._wp)then
-                    oldlayer = packet%layer
-
-                elseif(ds(packet%layer) > 0._wp .and. ds(new_layer) > 0._wp)then
-                    packet%tflag = .true.
-                    exit
-                else
+                !print*, packet%layer
+                !print*, pos
+                !print*, old_layer, new_layer
+                !print*, ds
+                !print*, dsNew
+
+                !print*, pos
+                !print*, smallStepPos
+                !print*, dir
+
+                layer = -1
+                if(dsNew(new_layer) < 0._wp .and. ds(new_layer) >= 0.0_wp) then
+                    !entered new layer
+                    layer = new_layer
+                else if (dsNew(old_layer) > 0._wp .and. ds(old_layer) <= 0.0_wp) then
+                    !left old_layer
+                    layer = old_layer
+                else 
                     error stop "This should not be reached!"
                 end if
-                if(oldlayer == 1)then
-                    oldlayer = packet%layer
-                else
-                    oldlayer = new_layer
-                end if
 
-                !print*, oldlayer
-                !print*, new_layer
-                !print*, pos
-                !print*, smallStepPos
-                N = calcNormal(pos, sdfs_array(oldlayer))
+                !print*, layer
+
+                N = calcNormal(pos, sdfs_array(layer))
                 !print*, N
 
                 rflag = .false.
                 call reflect_refract(dir, N, n1, n2, rflag, Ri)
-                packet%weight = packet%weight * Ri
-                ! Why on earth are we doing this
-                !tau = -log(ran2())
-                !taurun = 0._wp
+
                 if(.not.rflag)then
+                    !print*, "Transmitted"
                     !update layer and step across the boundary
+
+                    !move to boundary, then move a small step over the boundary in the new direction
+
                     packet%layer = new_layer
+                    oldpos = pos
+                    !comment out for phase screen
+                    call update_grids(grid, oldpos, dir, d_sdf, packet, sdfs_array(packet%layer)%getmua())
+                    t_sdf = d_sdf * sdfs_array(packet%layer)%getkappa()
+                    taurun = taurun + t_sdf
                     pos = smallStepPos
 
                     pointSep = sqrt((pos%x - startPos%x)**2 + (pos%y - startPos%y)**2 + (pos%z - startPos%z)**2)
@@ -280,17 +290,14 @@ subroutine tauint2(grid, packet, sdfs_array, dects, history)
                     end do
                     startPos = pos
                 else
+                    !print*, "Reflected"
                     !step back inside original sdf
-                    pos = oldpos
-                    startPos = oldpos
+                    oldpos = pos
+                    startPos = pos
 
                     !reflect so incrment bounce counter
                     packet%bounces = packet%bounces + 1
-                    ! Also why on earth are we doing this
-                    !if(packet%bounces > 1000)then
-                    !    packet%tflag=.true.
-                    !    exit
-                    !end if
+
                 end if
             else
                 ! n are equal, no change to the direction

src/kernelsMod.f90

@@ -505,7 +505,7 @@ subroutine pathlength_scatter2(input_file)
                 nscatt = nscatt + 1
                 packet%step = packet%step + 1
 
-                ! !Find next scattering location
+                ! Find next scattering location
                 call tauint2(state%grid, packet, array, dects, history)
             end do
 

src/parse/parse_source.f90

@@ -67,7 +67,7 @@ subroutine parse_source(table, packet, dict, spectrum, context, error)
 
             children => null()
 
-            if(state%source /= "unifrom" .and. state%source /= "point" .and. &
+            if(state%source /= "uniform" .and. state%source /= "point" .and. &
                 state%source /= "circular" .and. state%source /= "pencil")then
                 call get_value(child, "rotation", children, requested=.false., origin=origin)
                 if(associated(children))then