6 files changed, 308 insertions, 1 deletions
diff --git a/.gitignore b/.gitignore
index 5e3ad64..c037c97 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,4 @@
 *.xml
 *.h5
 *.pkl
-*.csv
 *-results/
diff --git a/apollonian-52.py b/apollonian-52.py
new file mode 100644
index 0000000..2b7fcd8
--- /dev/null
+++ b/apollonian-52.py
@@ -0,0 +1,197 @@
+import openmc
+import numpy as np
+import argparse as ap
+import csv
+
+#> Apollonian 52-element CANDU fuel bundle in OpenMC for AnnCon2026.
+#> Connor Moore, March 2026. <connor.moore@ontariotechu.net>
+
+### Argparser for parametric analysis ###
+
+parser = ap.ArgumentParser(prog="52-Element Apollonian Bundle Analysis",
+                           description="Variable pin radius for calculating ratios and criticality.")
+
+parser.add_argument("-fmr","--fuel-moderator-ratio",help="Ratio of the radius of fuel to whole circle. Default is 0.75",default="0.75")
+parser.add_argument("-bt","--bundle-type",help="Apollonian bundle type [normal/alt], default=normal",default="normal")
+args = parser.parse_args()
+
+### Material Definitions ###
+
+#> Materials marked PNNL-15870 are from the 2021 revision of "Compendium of Material
+#> Composition Data for Radiation Transport Modeling" published by the U.S. Department
+#> of Homeland Security and Pacific Northwest National Laboratory.
+#> <https://www.pnnl.gov/main/publications/external/technical_reports/PNNL-15870Rev2.pdf>
+
+#> The natural UO2 fuel
+mat_fuel = openmc.Material(name="Natural Uranium Fuel (UO2)")
+mat_fuel.add_element("U", 1.0, enrichment=0.71)
+mat_fuel.add_element("O", 2.0)
+mat_fuel.set_density("g/cc", 10.6)
+#> Density of 10.6 g/cc is from The Essential CANDU, Ch. 17 (Fuel), sec. 2.2, pp. 11
+
+#> Pressure tube and calandria tube use Zircaloy-2
+mat_zircaloy_2 = openmc.Material(name="Zircaloy-2 (PNNL-15870)")
+mat_zircaloy_2.add_nuclide("O16",   0.001194,"wo")
+mat_zircaloy_2.add_nuclide("O17",   0.000000,"wo")
+mat_zircaloy_2.add_nuclide("O18",   0.000003,"wo")
+mat_zircaloy_2.add_element("Cr",    0.000997,"wo")
+mat_zircaloy_2.add_element("Fe",    0.000997,"wo")
+mat_zircaloy_2.add_element("Ni",    0.000499,"wo")
+mat_zircaloy_2.add_nuclide("Zr90",  0.498109,"wo")
+mat_zircaloy_2.add_nuclide("Zr91",  0.109835,"wo")
+mat_zircaloy_2.add_nuclide("Zr92",  0.169730,"wo")
+mat_zircaloy_2.add_nuclide("Zr94",  0.175752,"wo")
+mat_zircaloy_2.add_nuclide("Zr96",  0.028918,"wo")
+mat_zircaloy_2.add_element("Sn",    0.013962,"wo")
+mat_zircaloy_2.set_density("g/cc",6.56)
+
+#> The fuel cladding uses Zircaloy-4
+mat_zircaloy_4 = openmc.Material(name="Zircaloy-4 (PNNL-15870")
+mat_zircaloy_4.add_nuclide("O16",   0.001193,"wo")
+mat_zircaloy_4.add_nuclide("O17",   0.000000,"wo")
+mat_zircaloy_4.add_nuclide("O18",   0.000003,"wo")
+mat_zircaloy_4.add_element("Cr",    0.000997,"wo")
+mat_zircaloy_4.add_element("Fe",    0.001993,"wo")
+mat_zircaloy_4.add_nuclide("Zr90",  0.497860,"wo")
+mat_zircaloy_4.add_nuclide("Zr91",  0.109780,"wo")
+mat_zircaloy_4.add_nuclide("Zr92",  0.169646,"wo")
+mat_zircaloy_4.add_nuclide("Zr94",  0.175665,"wo")
+mat_zircaloy_4.add_nuclide("Zr96",  0.028904,"wo")
+mat_zircaloy_4.add_element("Sn",    0.013955,"wo")
+mat_zircaloy_4.set_density("g/cc",6.56)
+
+#> Heavy water is used for moderation and cooling
+mat_d2o = openmc.Material(name="Heavy Water (PNNL-15870)")
+mat_d2o.add_nuclide("H2",   0.201133,"wo")
+mat_d2o.add_nuclide("O16",  0.796703,"wo")
+mat_d2o.add_nuclide("O17",  0.000323,"wo")
+mat_d2o.add_nuclide("O18",  0.001842,"wo")
+mat_d2o.add_s_alpha_beta("c_D_in_D2O")
+mat_d2o.set_density("g/cc",1.1044)
+
+materials = openmc.Materials([mat_fuel, mat_zircaloy_2, mat_zircaloy_4, mat_d2o])
+materials.export_to_xml()
+
+
+### Geometry Definition ###
+
+fuel_region_list = []
+clad_region_list = []
+fuel_rad_list = []
+fuel_area_list = []
+clad_area_list = []
+clad_circumference_list = []
+
+#> Define a function to create a fuel "pin" using a radius.
+
+def make_apollonian_pin(rg: float, rfm: float, x0: float, y0: float) -> None:
+    fuel_surf = openmc.ZCylinder(r=rg*rfm, x0=x0, y0=y0)
+    #> Hard-coded cladding thickness of 0.4 mm
+    clad_surf = openmc.ZCylinder(r=fuel_surf.r+0.04, x0=x0, y0=y0)
+    
+    fuel_region = -fuel_surf
+    clad_region = +fuel_surf & -clad_surf
+
+    fuel_region_list.append(fuel_region)
+    clad_region_list.append(clad_region)
+
+    fuel_area_list.append(np.pi*fuel_surf.r**2)
+    clad_area_list.append(np.pi*(clad_surf.r**2 - fuel_surf.r**2))
+
+    clad_circumference_list.append(2*np.pi*clad_surf.r)
+    return
+
+#> Import gasket points from file
+filename = "gasket-geometry/apollonian-gasket-52-normal.csv" if args.bundle_type=="normal" else "gasket-geometry/apollonian-gasket-52-alt.csv"
+rfm = eval(args.fuel_moderator_ratio)
+
+with open(filename, "r", newline="", encoding="utf-8") as file:
+    reader = csv.reader(file)
+
+    try:
+        header = next(reader)
+    except StopIteration:
+        header = []
+
+    for row in reader:
+        make_apollonian_pin(float(row[2]),rfm,float(row[0]),float(row[1]))
+
+#> Combine the regions to make a fuel cell and a cladding cell
+fuel_cell = openmc.Cell(name="UO2 Fuel Regions", region=openmc.Union(fuel_region_list), fill=mat_fuel)
+clad_cell = openmc.Cell(name="Zircaloy-4 Cladding Regions", region=openmc.Union(clad_region_list), fill=mat_zircaloy_4)
+
+#> Add the pressure tube
+pt_inner = openmc.ZCylinder(r=5.16890, x0=0.0, y0=0.0)
+pt_outer = openmc.ZCylinder(r=5.60320, x0=0.0, y0=0.0)
+
+pt_region = +pt_inner & -pt_outer
+pt_cell = openmc.Cell(name="Pressure Tube", region=pt_region, fill=mat_zircaloy_2)
+
+#> Pack the fuel with D2O
+pt_d2o_region = ~(fuel_cell.region | clad_cell.region) & -pt_inner
+pt_d2o_cell = openmc.Cell(name="D2O Coolant", region=pt_d2o_region, fill=mat_d2o)
+
+#> Add the calandria tube
+ct_inner = openmc.ZCylinder(r=6.44780, x0=0.0, y0=0.0)
+ct_outer = openmc.ZCylinder(r=6.58750, x0=0.0, y0=0.0)
+
+ct_region = +ct_inner & -ct_outer
+ct_cell = openmc.Cell(name="Calandria Tube", region=ct_region, fill=mat_zircaloy_2)
+
+#> The space between the calandria tube and pressure tube is considered void
+pt_ct_region = +pt_outer & -ct_inner
+pt_ct_cell = openmc.Cell(name="Annulus Gap", region=pt_ct_region, fill=None)
+
+#> The lattice pitch for the assembly is 28.575, so apply a reflecting boundary
+outer_boundary = openmc.model.RectangularPrism(width=28.575, height=28.575, origin=(0.0, 0.0), boundary_type="reflective")
+
+ct_d2o_region = +ct_outer & -outer_boundary
+ct_d2o_cell = openmc.Cell(name="D2O Moderator", region=ct_d2o_region, fill=mat_d2o)
+
+universe = openmc.Universe(cells=[fuel_cell, clad_cell, pt_cell, pt_d2o_cell, ct_cell, pt_ct_cell, ct_d2o_cell])
+
+geometry = openmc.Geometry(universe)
+geometry.export_to_xml()
+
+
+### Settings definition ###
+settings = openmc.Settings()
+settings.particles = 10000
+settings.batches = 200
+settings.inactive = 80
+
+#> Set up the source to sample inside the pressure tube region uniformly
+settings.source = openmc.IndependentSource()
+settings.source.space = openmc.stats.CylindricalIndependent(
+        r = openmc.stats.PowerLaw(a=0, b=pt_inner.r, n=1),
+        phi = openmc.stats.Uniform(0, 2*np.pi),
+        z = openmc.stats.Discrete([0.0], [1.0])
+        )
+settings.source.energy = openmc.stats.Watt()
+
+settings.export_to_xml()
+
+
+### Area and DTU Calculations ###
+
+V_fuel = sum(fuel_area_list)
+V_clad = sum(clad_area_list)
+C_clad = sum(clad_circumference_list)
+V_mod = np.pi*pt_inner.r**2 - (V_fuel + V_clad) + (28.575**2 - np.pi*ct_outer.r**2)
+
+d_fuel = mat_fuel.get_nuclide_atom_densities()
+N_fuel = sum(density for nuclide,density in d_fuel.items() if "U" in nuclide)
+
+N_mod = mat_d2o.get_nuclide_atom_densities()["H2"]
+DTU_ratio = (V_mod * N_mod) / (V_fuel * N_fuel)
+
+##> Print these for the final table
+print(f"Flow area, cm² = {V_mod}")
+print(f"Cladding ciricumference, cm = {C_clad}")
+print(f"Fuel mass per length, g/cm = {V_fuel*mat_fuel.density}")
+print(f"Cladding mass per length, g/cm = {V_clad*mat_zircaloy_4.density}")
+print(f"DTU ratio = {DTU_ratio}")
+
+
+#> Run the model!
+openmc.run()
diff --git a/candu-37.py b/candu-37.py
index e9642e7..b58dd68 100644
--- a/candu-37.py
+++ b/candu-37.py
@@ -108,9 +108,12 @@ def make_ring(n: int, r: float, alpha: float) -> None:
         theta = (i*(360/n) + alpha)*np.pi/180
         x = r*np.cos(theta)
         y = r*np.sin(theta)
+        #print((float(x),float(y)))
 
         make_pin(rf=r_fuel, rc=r_clad, x0=x, y0=y)
 
+    return
+
 
 #> Start with the innermost pin
 make_pin(rf=r_fuel, rc=r_clad, x0=0.0, y0=0.0)
@@ -184,6 +187,7 @@ settings.export_to_xml()
 
 V_fuel = len(fuel_region_list)*np.pi*r_fuel**2
 V_clad = len(clad_region_list)*np.pi*(r_clad**2 - r_fuel**2)
+C_clad = len(clad_region_list)*2*np.pi*r_clad
 V_mod = np.pi*pt_inner.r**2 - (V_fuel + V_clad) + (28.575**2 - np.pi*ct_outer.r**2)
 
 d_fuel = mat_fuel.get_nuclide_atom_densities()
@@ -194,6 +198,7 @@ DTU_ratio = (V_mod * N_mod) / (V_fuel * N_fuel)
 
 #> Print these for the final table
 print(f"Flow area, cm² = {V_mod}")
+print(f"Cladding ciricumference, cm = {C_clad}")
 print(f"Fuel mass per length, g/cm = {V_fuel*mat_fuel.density}")
 print(f"Cladding mass per length, g/cm = {V_clad*mat_zircaloy_4.density}")
 print(f"DTU ratio = {DTU_ratio}")
diff --git a/driver-37.sh b/driver-37.sh
index 31a81f5..c365623 100755
--- a/driver-37.sh
+++ b/driver-37.sh
@@ -22,4 +22,5 @@ for r in $(seq 0.1 0.05 0.8); do
     dtu=$(cat $filename | grep "DTU ratio" | cut -d "=" -f 2)
 
     echo $keff +/- $pm, $flow_area cm² flow, $fuel_mass g/cm fuel, $clad_mass g/cm clad, $dtu DTU ratio | cowsay
+    echo $keff $pm $flow_area $fuel_mass $clad_mass $dtu >> 37-results/37_keff_table
 done
diff --git a/gasket-geometry/apollonian-gasket-52-alt.csv b/gasket-geometry/apollonian-gasket-52-alt.csv
new file mode 100644
index 0000000..2038ce1
--- /dev/null
+++ b/gasket-geometry/apollonian-gasket-52-alt.csv
@@ -0,0 +1,52 @@
+Center_X,Center_Y,Radius
+-2.752424,3.496048,0.495000
+-4.403878,0.635645,0.495000
+-1.651454,-4.131693,0.495000
+1.651454,-4.131693,0.495000
+2.752424,3.496048,0.495000
+4.403878,0.635645,0.495000
+-1.559446,3.579838,0.495000
+-0.0,0.878799,0.495000
+1.559446,3.579838,0.495000
+-3.879954,-0.4394,0.495000
+-2.320508,-3.140439,0.495000
+-0.761062,-0.4394,0.495000
+0.761062,-0.4394,0.495000
+2.320508,-3.140439,0.495000
+3.879954,-0.4394,0.495000
+-3.360139,2.540208,0.495000
+-3.879954,1.639862,0.495000
+-2.840324,1.639862,0.495000
+0.0,-3.279723,0.495000
+-0.519815,-4.180069,0.495000
+0.519815,-4.180069,0.495000
+3.360139,2.540208,0.495000
+2.840324,1.639862,0.495000
+3.879954,1.639862,0.495000
+0.0,4.500185,0.495000
+-0.499815,3.63448,0.495000
+0.499815,3.63448,0.495000
+-1.076952,2.634851,0.495000
+-1.576767,1.769146,0.495000
+-0.577137,1.769146,0.495000
+1.076952,2.634851,0.495000
+0.577137,1.769146,0.495000
+1.576767,1.769146,0.495000
+-2.320508,0.480947,0.495000
+-2.820323,-0.384758,0.495000
+-1.820693,-0.384758,0.495000
+-3.39746,-1.384387,0.495000
+-3.897275,-2.250092,0.495000
+-2.897645,-2.250092,0.495000
+-1.243556,-1.384387,0.495000
+-1.743371,-2.250092,0.495000
+-0.743741,-2.250092,0.495000
+2.320508,0.480947,0.495000
+1.820693,-0.384758,0.495000
+2.820323,-0.384758,0.495000
+1.243556,-1.384387,0.495000
+0.743741,-2.250092,0.495000
+1.743371,-2.250092,0.495000
+3.39746,-1.384387,0.495000
+2.897645,-2.250092,0.495000
+3.897275,-2.250092,0.495000
diff --git a/gasket-geometry/apollonian-gasket-52-normal.csv b/gasket-geometry/apollonian-gasket-52-normal.csv
new file mode 100644
index 0000000..78ddebe
--- /dev/null
+++ b/gasket-geometry/apollonian-gasket-52-normal.csv
@@ -0,0 +1,53 @@
+Center_X,Center_Y,Radius
+0.0,0.0,0.358984
+-2.752424,3.496048,0.550485
+-4.403878,0.635645,0.550485
+-1.651454,-4.131693,0.550485
+1.651454,-4.131693,0.550485
+2.752424,3.496048,0.550485
+4.403878,0.635645,0.550485
+-1.559446,3.579838,0.519815
+-0.0,0.878799,0.519815
+1.559446,3.579838,0.519815
+-3.879954,-0.4394,0.519815
+-2.320508,-3.140439,0.519815
+-0.761062,-0.4394,0.519815
+0.761062,-0.4394,0.519815
+2.320508,-3.140439,0.519815
+3.879954,-0.4394,0.519815
+-3.360139,2.540208,0.519815
+-3.879954,1.639862,0.519815
+-2.840324,1.639862,0.519815
+0.0,-3.279723,0.519815
+-0.519815,-4.180069,0.519815
+0.519815,-4.180069,0.519815
+3.360139,2.540208,0.519815
+2.840324,1.639862,0.519815
+3.879954,1.639862,0.519815
+0.0,4.500185,0.499815
+-0.499815,3.63448,0.499815
+0.499815,3.63448,0.499815
+-1.076952,2.634851,0.499815
+-1.576767,1.769146,0.499815
+-0.577137,1.769146,0.499815
+1.076952,2.634851,0.499815
+0.577137,1.769146,0.499815
+1.576767,1.769146,0.499815
+-2.320508,0.480947,0.499815
+-2.820323,-0.384758,0.499815
+-1.820693,-0.384758,0.499815
+-3.39746,-1.384387,0.499815
+-3.897275,-2.250092,0.499815
+-2.897645,-2.250092,0.499815
+-1.243556,-1.384387,0.499815
+-1.743371,-2.250092,0.499815
+-0.743741,-2.250092,0.499815
+2.320508,0.480947,0.499815
+1.820693,-0.384758,0.499815
+2.820323,-0.384758,0.499815
+1.243556,-1.384387,0.499815
+0.743741,-2.250092,0.499815
+1.743371,-2.250092,0.499815
+3.39746,-1.384387,0.499815
+2.897645,-2.250092,0.499815
+3.897275,-2.250092,0.499815