From c5d3fe949c1354e89c228e67f75af26af1ae0502 Mon Sep 17 00:00:00 2001
From: Connor Moore <connor@hhmoore.ca>
Date: Mon, 6 Apr 2026 05:11:16 -0400
Subject: Added basic tracking and transport for main fortran file, makefile
 for building/running/plotting with gnuplot (wip)

---
 Makefile |  16 ++++++++
 main.f90 | 132 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 146 insertions(+), 2 deletions(-)
 create mode 100644 Makefile

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..96a64c3
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,16 @@
+FC = gfortran
+FFLAGS = -march=native -Ofast
+XNAME = monte_carlo.x
+NPROC := $(shell nproc)
+
+
+all: build run plot
+
+build: main.f90
+	$(FC) $(FFLAGS) main.f90 -o $(XNAME)
+
+run: build
+	OMP_NUM_THREADS=$(NPROC) ./$(XNAME)
+
+plot: run
+	echo test
diff --git a/main.f90 b/main.f90
index 05f5219..1be36c4 100644
--- a/main.f90
+++ b/main.f90
@@ -1,7 +1,135 @@
-program main
+program monte_carlo
 
     use, intrinsic :: iso_fortran_env
     implicit none
 
+    !> Initialize some variables for neutrons
+    integer(int64), parameter :: n = 1e6
+    integer(int64), parameter :: bins = 800
 
-end program main
+    !> And for geometry
+    real(real64), parameter :: L = 8.0_real64
+    real(real64), parameter :: dx = L / bins
+    real(real64), dimension(6), parameter :: boundaries = [0.0_real64, 2.0_real64, 3.0_real64, 5.0_real64, 6.0_real64, 8.0_real64]
+
+    !> These will contain the x-position and cosine of angle for each neutron
+    real(real64), dimension(n) :: px, mu
+    !> Track whether or not the neutron is alive
+    logical, dimension(n) :: is_alive
+
+    !> Some local trackers to be updated via subroutine call
+    real(real64) :: sigma_t, sigma_s, q, d, d_boundary, xi, next_boundary
+    integer(int64) :: i, region
+
+    !> Initialize RNG to be actually random
+    call random_seed()
+    call initialize_neutrons(n, px, mu)
+
+    !> The main monte carlo loop 
+    main: do i = 1,n
+
+        single: do while (is_alive(i))
+
+            !> Find what region the neutron is in and get XS
+            region = count(px(i).ge.boundaries)
+            call get_region_info(region, sigma_t, sigma_s, q)
+
+            !> Determine the distance before collision
+            if(sigma_t.eq.0.0_real64) then
+                !> If we are in a void, make it huge
+                d = 1e10_real64
+            else
+                !> Else, calculate with d = -ln(xi)/sigma_t
+                call random_number(xi)
+                d = -log(xi) / sigma_t
+            end if
+
+            !> Determine where the next boundary is and the distance
+            if(mu(i).gt.0.0_real64) then
+                d_boundary = (boundaries(region+1)-px(i))/mu(i)
+            else if(mu(i).lt.0.0_real64) then
+                d_boundary = (boundaries(region)-px(i))/mu(i)
+            else
+                !> In the rare chance it is perfectly parallel with the walls
+                d_boundary = 1e10_real64
+                is_alive(i) = .false.
+            end if
+
+            !> If it goes through a wall, restart the particle at the new boundary
+            if(d_boundary.le.d) then
+                px(i) = px(i) + d_boundary*mu(i)
+                if(px(i).le.0.0_real64) then
+                    px(i) = 0.0_real64
+                    mu(i) = -mu(i)
+                else if (px(i).ge.L) then
+                    is_alive = .false.
+                    exit single
+                end if
+                cycle single
+            end if
+
+            !> If not, determine what type of interaction happens
+            call random_number(xi)
+            if(xi.le.sigma_s/sigma_t) then
+                !> Scattering interaction
+                !> Sample new isotropic angle
+                call random_number(xi)
+                mu(i) = 2.0_real64*xi - 1.0_real64
+            else
+                !> Absorption interaction
+                is_alive(i) = .false.
+            end if
+
+        end do single
+
+
+    end do main
+
+
+contains
+        
+    pure elemental subroutine get_region_info(region, sigma_t, sigma_s, q)
+        integer(int64), intent(in) :: region
+        real(real64), intent(out) :: sigma_t, sigma_s, q
+
+        if(region.eq.1) then
+            sigma_t = 50.0_real64; sigma_s = 0.0_real64; q = 50.0_real64
+        else if(region.eq.2) then
+            sigma_t = 5.0_real64; sigma_s = 0.0_real64; q = 5.0_real64
+        else if(region.eq.3) then
+            sigma_t = 0.0_real64; sigma_s = 0.0_real64; q = 0.0_real64
+        else if(region.eq.4) then
+            sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 1.0_real64
+        else if(region.eq.5) then
+            sigma_t = 1.0_real64; sigma_s = 0.9_real64; q = 0.0_real64
+        end if
+    end subroutine get_region_info
+
+    impure subroutine initialize_neutrons(n, px, mu)
+        integer(int64), intent(in) :: n
+        real(real64), intent(out), dimension(n) :: px, mu
+        real(real64), dimension(n) :: r_region, r_position, r_mu
+
+        !> Make some uniform random numbers
+        call random_number(r_region)
+        call random_number(r_position)
+        call random_number(r_mu)
+
+        !> Consider the total integrated source-distance for regions
+        r_region = r_region * 106.0_real64
+
+        !> Assign a random position in each source region
+        where(r_region.lt.100)
+            px = r_position * 2.0_real64
+        else where(r_region.lt.105)
+            px = 2.0_real64 + r_position
+        else where
+            px = 5.0_real64 + r_position
+        end where
+
+        !> Assign random angles to each particle
+        mu = 2.0_real64 * r_mu - 1.0_real64
+
+    end subroutine initialize_neutrons
+
+end program monte_carlo
-- 
cgit v1.2.3