Fixed statistics, pdf_tex plot for presenting

1 files changed, 32 insertions, 23 deletions

diff --git a/main.f90 b/main.f90
index d3588b5..cd61bc5 100644
--- a/main.f90
+++ b/main.f90
@@ -1,11 +1,10 @@
 program monte_carlo
 
-    use omp_lib
     use, intrinsic :: iso_fortran_env
     implicit none
 
     !> Initialize some variables for neutrons
-    integer(int64), parameter :: n = 1e8
+    integer(int64), parameter :: n = 1e7
     integer(int64), parameter :: bins = 800
     real(real64), parameter :: eps = 1.0e-6_real64
     real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64)
@@ -20,16 +19,18 @@ program monte_carlo
     !> Track whether or not the neutron is alive
     logical, allocatable, dimension(:) :: is_alive
     !> And for tallying the flux
-    real(real64), allocatable, dimension(:) :: flux_tally
+    real(real64), allocatable, dimension(:) :: neutron_bin_length, flux_tally, flux_sq_tally, flux_std_dev
 
     !> Some local trackers to be updated via subroutine call
     real(real64) :: sigma_t, sigma_s, d_flight, d_boundary, d_step, xi, avg_r1_flux
     integer(int64) :: i, region
 
     !> Allocate arrays
-    allocate(x(n), u(n), is_alive(n), flux_tally(bins))
+    allocate(x(n), u(n), is_alive(n), neutron_bin_length(bins), flux_tally(bins), flux_sq_tally(bins), flux_std_dev(bins))
     is_alive = .true.
     flux_tally = 0.0_real64
+    flux_sq_tally = 0.0_real64
+    flux_std_dev = 0.0_real64
 
     !> Initialize RNG to be actually random
     call random_seed()
@@ -40,6 +41,8 @@ program monte_carlo
 
     !> The main monte carlo loop 
     main: do i = 1,n
+        !> Initialize empty neutron length
+        neutron_bin_length = 0.0_real64
 
         !> Simulate the neutron until it is lost
         single: do while(is_alive(i))
@@ -49,7 +52,7 @@ program monte_carlo
             call get_region_info(region, sigma_t, sigma_s)
 
             if(region .ge. size(boundaries) .or. region .le. 0) then
-                !> We are outside of the problem. Not good. Skip.
+                !> Put neutrons back inside if they get lost
                 u(i) = eps
             end if
 
@@ -79,7 +82,7 @@ program monte_carlo
             d_step = min(d_flight,d_boundary)
 
             !> Tally the track-length contributions for the bin
-            call tally_segment(x(i), u(i), d_step, bins, bin_dx, flux_tally)
+            call tally_segment(x(i), u(i), d_step, bins, bin_dx, neutron_bin_length)
 
             !> Update the neutron position
             x(i) = x(i) + d_step * u(i)
@@ -90,7 +93,6 @@ program monte_carlo
                 if(u(i) .gt. 0.0_real64 .and. abs(x(i)-L) .lt. eps) then
                     !> Bin neutron at vacuum boundary at x=8.0
                     is_alive(i) = .false.
-                    cycle main
                 else if(u(i) .lt. 0.0_real64 .and. x(i) .lt. eps) then
                     !> Reflect on specular boundary at x=0.0
                     u(i) = -u(i)
@@ -113,23 +115,30 @@ program monte_carlo
                 else
                     !> The particle is absorbed and lost
                     is_alive(i) = .false.
-                    cycle main
                 end if
             end if
 
         end do single
 
+        !> Update tallies
+        flux_tally = flux_tally + neutron_bin_length
+        flux_sq_tally = flux_sq_tally + neutron_bin_length**2
+
     end do main
 
+    !> sqrt(x²/n - (x/n)²)/(dx²)
+    flux_std_dev = sqrt((flux_sq_tally/real(n, real64) - (flux_tally/real(n, real64))**2)) / bin_dx**2
+
     !> Average region 1 flux
     avg_r1_flux = sum(flux_tally(1:floor(boundaries(2)/bin_dx)))/real(floor(boundaries(2)/bin_dx), real64)
 
     !> Normalize to region 1
     flux_tally = flux_tally / avg_r1_flux
+    flux_std_dev = flux_std_dev / avg_r1_flux
 
     open(unit=10, file='flux.out', status='replace')
     write: do i = 1,bins
-        write(10,'(2E15.7)') (i-0.5_real64)*bin_dx, flux_tally(i)
+        write(10,'(3E15.7)') (i-0.5_real64)*bin_dx, flux_tally(i), flux_std_dev(i)
     end do write
 
     !> Deallocate everyting
@@ -181,12 +190,12 @@ contains
 
     end subroutine initialize_neutrons
 
-    impure subroutine tally_segment(x, u, d_step, bins, bin_dx, flux_tally)
+    pure subroutine tally_segment(x, u, d_step, bins, bin_dx, flux_tally)
         integer(int64), intent(in) :: bins
         real(real64), intent(in) :: x, u, d_step, bin_dx
         real(real64), intent(in out), dimension(bins) :: flux_tally
         real(real64) :: x_start, x_end, dx
-        integer(int64) :: i_start, i_end, i
+        integer(int64) :: start_bin, end_bin, bin
 
         x_start = x
         x_end = x + d_step*u
@@ -194,28 +203,28 @@ contains
             call swap(x_start, x_end)
         end if
 
-        i_start = ceiling(x_start/bin_dx)
-        i_end   = ceiling(x_end/bin_dx)
+        start_bin = ceiling(x_start/bin_dx)
+        end_bin   = ceiling(x_end/bin_dx)
 
         !> Clamp to valid bin range
-        i_start = max(1, min(bins, i_start))
-        i_end   = max(1, min(bins, i_end))
+        start_bin = max(1, min(bins, start_bin))
+        end_bin   = max(1, min(bins, end_bin))
 
-        if (i_start .eq. i_end) then
-            flux_tally(i_start) = flux_tally(i_start) + d_step
+        if (start_bin .eq. end_bin) then
+            flux_tally(start_bin) = flux_tally(start_bin) + d_step
         else
             !> First partial bin
-            dx = (i_start*bin_dx - x_start)/abs(u)
-            flux_tally(i_start) = flux_tally(i_start) + dx
+            dx = (start_bin*bin_dx - x_start)/abs(u)
+            flux_tally(start_bin) = flux_tally(start_bin) + dx
 
             !> Full bins in between
-            do i = i_start+1, i_end-1
-                flux_tally(i) = flux_tally(i) + bin_dx/abs(u)
+            do bin = start_bin+1, end_bin-1
+                flux_tally(bin) = flux_tally(bin) + bin_dx/abs(u)
             end do
 
             !> Last partial bin
-            dx = (x_end - (i_end-1)*bin_dx)/abs(u)
-            flux_tally(i_end) = flux_tally(i_end) + dx
+            dx = (x_end - (end_bin-1)*bin_dx)/abs(u)
+            flux_tally(end_bin) = flux_tally(end_bin) + dx
         end if
 
     end subroutine tally_segment