diff --git a/app/main.f90 b/app/main.f90 index 01708cd..d9e4c98 100644 --- a/app/main.f90 +++ b/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 diff --git a/res/scat_test.toml b/res/scat_test.toml index 45a6a3f..a9484c8 100644 --- a/res/scat_test.toml +++ b/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" diff --git a/res/validation1.toml b/res/validation1.toml new file mode 100644 index 0000000..ee65453 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/res/validation2.toml b/res/validation2.toml new file mode 100644 index 0000000..dffa1bf --- /dev/null +++ b/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 \ No newline at end of file diff --git a/res/validation3.toml b/res/validation3.toml new file mode 100644 index 0000000..ae0a87d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/inttau2.f90 b/src/inttau2.f90 index acf9b79..b6da407 100644 --- a/src/inttau2.f90 +++ b/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 diff --git a/src/kernelsMod.f90 b/src/kernelsMod.f90 index e9686e8..c284327 100644 --- a/src/kernelsMod.f90 +++ b/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 diff --git a/src/parse/parse_source.f90 b/src/parse/parse_source.f90 index bc4f602..3535ef6 100644 --- a/src/parse/parse_source.f90 +++ b/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 diff --git a/src/surfaces.f90 b/src/surfaces.f90 index ca4adcd..cc3f4fc 100644 --- a/src/surfaces.f90 +++ b/src/surfaces.f90 @@ -98,7 +98,6 @@ function fresnel(I, N, n1, n2) result (tir) real(kind=wp) :: costt, sintt, sint2, cost2, tir, f1, f2 costt = abs(I .dot. N) - sintt = sqrt(1._wp - costt * costt) sint2 = n1/n2 * sintt if(sint2 > 1._wp)then @@ -112,9 +111,11 @@ function fresnel(I, N, n1, n2) result (tir) cost2 = sqrt(1._wp - sint2 * sint2) f1 = abs((n1*costt - n2*cost2) / (n1*costt + n2*cost2))**2 f2 = abs((n1*cost2 - n2*costt) / (n1*cost2 + n2*costt))**2 - + tir = 0.5_wp * (f1 + f2) - if(ieee_is_nan(tir) .or. tir > 1._wp .or. tir < 0._wp)print*,'TIR: ', tir, f1, f2, costt,sintt,cost2,sint2 + end if + if(ieee_is_nan(tir) .or. tir > 1._wp .or. tir < 0._wp) then + print*,'TIR: ', tir, f1, f2, costt,sintt,cost2,sint2 return end if end function fresnel diff --git a/src/writer.f90 b/src/writer.f90 index c42f1dd..42f0d78 100644 --- a/src/writer.f90 +++ b/src/writer.f90 @@ -242,13 +242,9 @@ subroutine write_hdr(u, sizes, type) integer :: i string = "" - do i = 1, size(sizes) - if(i == 1)then - string = str(sizes(i)) - else - string = trim(string) // " " // str(sizes(i)) - end if - end do + string = str(sizes(3)) + string = trim(string) // " " // str(sizes(2)) + string = trim(string) // " " // str(sizes(1)) write(u,"(A)")"NRRD0004" write(u,"(A)")"type: "//type diff --git a/tools/plotDetectors.py b/tools/plotDetectors.py index 81b2864..8e3f7df 100644 --- a/tools/plotDetectors.py +++ b/tools/plotDetectors.py @@ -4,7 +4,7 @@ folderName = "RSMCRT/data/detectors/" -filename = folderName + "detector_1.dat" +filename = folderName + "detector_2.dat" NoOpticalDepths = 5 OpticalDepths = [0.1, 1, 10, 30, 100] @@ -39,9 +39,9 @@ #print(count) totalCounts = sum(count) -#print(count) +print(count) #print("Total counts: " + str(totalCounts)) -print("Total Diffuse ...: " + str(totalCounts/1000000)) +print("Total Diffuse ... : " + str(totalCounts/1000000)) fig = plt.figure(1) @@ -51,4 +51,4 @@ ax1.set_xlabel("Radius (m)") ax1.set_ylabel("Counts (Arb. Units)") ax1.set_title("Radius Detector") -#plt.show() \ No newline at end of file +plt.show() \ No newline at end of file diff --git a/tools/plotGrid.py b/tools/plotGrid.py new file mode 100644 index 0000000..a3d6265 --- /dev/null +++ b/tools/plotGrid.py @@ -0,0 +1,126 @@ +from collections import OrderedDict +import os +import re + +import numpy as np + + +_types = {"float": np.float32, "double": np.float64, + "int32": np.int32, "int64": np.int64, + "uchar": np.ubyte} + + +def read_nrrd(file): + with open(file, "rb") as fh: + hdr = read_header(fh) + data = read_data(fh, hdr) + return data, hdr + + +def read_header(file): + + it = iter(file) + magic_line = next(it) + if hasattr(magic_line, "decode"): + need_decode = True + magic_line = magic_line.decode("ascii", "ignore") # type: ignore[assignment] + if not magic_line.startswith("NRRD"): # type: ignore[arg-type] + raise NotImplementedError + + header = OrderedDict() + + for line in it: + if need_decode: + line = line.decode("ascii", "ignore") # type: ignore[assignment] + + line = line.rstrip() + if line.startswith("#"): # type: ignore[arg-type] + continue + elif line == "": + break + + key, value = re.split(r"[:=?]", line, 1) # type: ignore[type-var] + key, value = key.strip(), value.strip() # type: ignore[attr-defined] + + value = _get_value_type(key, value) + + header[key] = value + + return header + + +def _get_value_type(key, value): + + if key in ["dimension", "space dimension"]: + return int(value) + elif key in ["endian", "encoding"]: + return value + elif key in ["type"]: + return _types[value] + elif key in ["sizes"]: + return [int(x) for x in value.split()] + else: + pass + # raise NotImplementedError + + +def read_data(file, header): + + size = np.array(header["sizes"]) + dtype = header["type"] + + total_data_points = size.prod(dtype=np.int64) + dtype_size = np.dtype(dtype).itemsize + file.seek(-1 * dtype_size * total_data_points, os.SEEK_END) + data = np.fromfile(file, dtype=dtype, sep="") + data = data.reshape((size)) + return data + + +if __name__ == '__main__': + import matplotlib.pyplot as plt + folderName = "RSMCRT/data/absorb/" + fileName = "absorb.nrrd" + file = folderName + fileName + grid, hdr = read_nrrd(file) + + + fig = plt.figure(1) + ax1 = fig.add_subplot() + + depths = np.linspace(-2.0, 2.0, int(hdr["sizes"][0])) + ymid = int(hdr["sizes"][1]/2) + zmid = int(hdr["sizes"][2]/2) + data = grid + + fluence = np.mean(np.mean(data, axis =2), axis = 1) + + ax1.plot(depths, fluence) + ax1.set_xlabel("Depth (cm)") + ax1.set_ylabel("Fluence (-)") + ax1.set_xlim([depths[-14],1.6]) + ax1.set_ylim([0, np.max(fluence)*1.1]) + + """ Used for validate 1 + c1 = 5.76 + k1 = 1.00 + c2 = 1.31 + k2 = 10.2 + delta = 0.047 + norm = 0.115 + #""" + + #""" Used for validate 2 + c1 = 6.27 + k1 = 1.00 + c2 = 1.18 + k2 = 14.4 + delta = 0.261 + norm = 0.0151 + #""" + + fittingFunction = norm * (c1* np.exp((depths-1.95)*k1/delta) - c2*np.exp((depths-1.95)*k2/delta)) + ax1.plot(depths, fittingFunction) + + plt.show() +