#!/usr/bin/env python3

"""
      ====================================================
      Description: ACR1000 Assembly 43 Fuel pins
      Case:        Problem  (43 fuel pins)
      Written By:  Ahmed K Madani
      Conversion:  Adam Parler
      Date:        12/03/2023
      ====================================================
"""

import os
import openmc

try:
    os.mkdir(os.path.dirname(os.path.realpath(__file__)) + '/problem')
except OSError: pass

################## MATERIALS ##################

fuel_21_hot = openmc.Material(material_id=11)
fuel_21_hot.set_density("g/cc", 9.82)
fuel_21_hot.temperature = 1200 # K (42c=300K, 50c=293.6, 53c=587.2, 66=293.6, 70c=293.6, 71c=600K, 72c=900K, 73c=1200K, 74c=2500K)
fuel_21_hot.add_nuclide("U235", 2.1, percent_type="wo")
fuel_21_hot.add_nuclide("U238", 97.9, percent_type="wo")
fuel_21_hot.add_nuclide("O16", 13.4529, percent_type="wo")

fuel_21_mid = openmc.Material(material_id=12)
fuel_21_mid.set_density("g/cc", 9.82)
fuel_21_mid.temperature = 900
fuel_21_mid.add_nuclide("U235", 2.1, percent_type="wo")
fuel_21_mid.add_nuclide("U238", 97.9, percent_type="wo")
fuel_21_mid.add_nuclide("O16", 13.4529, percent_type="wo")

fuel_21_low = openmc.Material(material_id=13)
fuel_21_low.set_density("g/cc", 10.298)
fuel_21_low.temperature = 900
fuel_21_low.add_nuclide("U235", 2.1, percent_type="wo")
fuel_21_low.add_nuclide("U238", 97.9, percent_type="wo")
fuel_21_low.add_nuclide("O16", 13.4529, percent_type="wo")

fuel_21_dep = openmc.Material(material_id=14)
fuel_21_dep.set_density("g/cc", 10.298)
fuel_21_dep.temperature = 900
fuel_21_dep.add_nuclide("U235", 0.72, percent_type="wo")
fuel_21_dep.add_nuclide("U238", 99.28, percent_type="wo")
fuel_21_dep.add_nuclide("O16", 13.4529, percent_type="wo")
fuel_21_dep.add_nuclide("Dy160", 0.1735, percent_type="wo")
fuel_21_dep.add_nuclide("Dy161", 1.418, percent_type="wo")
fuel_21_dep.add_nuclide("Dy162", 1.915, percent_type="wo")
fuel_21_dep.add_nuclide("Dy163", 1.861, percent_type="wo")
fuel_21_dep.add_nuclide("Dy164", 2.1317, percent_type="wo")

clad = openmc.Material(material_id=2)
clad.set_density("g/cc", 7.48)
clad.temperature = 600
clad.add_nuclide("O16", 0.125, percent_type="wo")
clad.add_nuclide("Cr50", 0.004345, percent_type="wo")
clad.add_nuclide("Cr52", 0.083789, percent_type="wo")
clad.add_nuclide("Cr53", 0.009501, percent_type="wo")
clad.add_nuclide("Cr54", 0.002365, percent_type="wo")
clad.add_nuclide("Sn112", 0.014065, percent_type="wo")
clad.add_nuclide("Sn114", 0.00957, percent_type="wo")
clad.add_nuclide("Sn115", 0.00493, percent_type="wo")
clad.add_nuclide("Sn116", 0.21083, percent_type="wo")
clad.add_nuclide("Sn117", 0.11136, percent_type="wo")
clad.add_nuclide("Sn118", 0.35119, percent_type="wo")
clad.add_nuclide("Sn119", 0.124555, percent_type="wo")
clad.add_nuclide("Sn120", 0.47241, percent_type="wo")
clad.add_nuclide("Sn122", 0.067135, percent_type="wo")
clad.add_nuclide("Sn124", 0.083955, percent_type="wo")
clad.add_nuclide("Fe54", 0.0122745, percent_type="wo")
clad.add_nuclide("Fe56", 0.1926834, percent_type="wo")
clad.add_nuclide("Fe57", 0.0044499, percent_type="wo")
clad.add_nuclide("Fe58", 0.0005922, percent_type="wo")
clad.add_nuclide("Zr90", 50.48, percent_type="wo")
clad.add_nuclide("Zr91", 11, percent_type="wo")
clad.add_nuclide("Zr92", 16.83, percent_type="wo")
clad.add_nuclide("Zr94", 17.05, percent_type="wo")
clad.add_nuclide("Zr96", 2.75, percent_type="wo")

water = openmc.Material(material_id=3)
water.set_density("g/cc", 0.71273)
water.temperature = 600
water.add_nuclide("H1", 2/3, percent_type="ao")
water.add_nuclide("O16", 1/3, percent_type="ao")
water.add_s_alpha_beta(name='c_H_in_H2O')

tube = openmc.Material(material_id=4)
tube.set_density("g/cc", 6.57)
tube.temperature = 600
tube.add_nuclide("Zr90", 50.16, percent_type="wo")
tube.add_nuclide("Zr91", 10.94, percent_type="wo")
tube.add_nuclide("Zr92", 16.72, percent_type="wo")
tube.add_nuclide("Zr94", 16.945, percent_type="wo")
tube.add_nuclide("Zr96", 2.73, percent_type="wo")
tube.add_nuclide("Nb93", 2.5, percent_type="wo")

graphite = openmc.Material(material_id=5)
graphite.set_density("g/cc", 0.0017)
graphite.add_nuclide("O16", 2/3, percent_type="ao")
graphite.add_nuclide("C0", 1/3, percent_type="ao")
graphite.add_s_alpha_beta(name='c_Graphite')

cal_zirc = openmc.Material(material_id=6)
cal_zirc.set_density("g/cc", 6.44)
cal_zirc.temperature = 600
cal_zirc.add_nuclide("O16", 0.125, percent_type="wo")
cal_zirc.add_nuclide("Cr50", 0.004345, percent_type="wo")
cal_zirc.add_nuclide("Cr52", 0.083789, percent_type="wo")
cal_zirc.add_nuclide("Cr53", 0.009501, percent_type="wo")
cal_zirc.add_nuclide("Cr54", 0.002365, percent_type="wo")
cal_zirc.add_nuclide("Sn112", 0.014065, percent_type="wo")
cal_zirc.add_nuclide("Sn114", 0.00957, percent_type="wo")
cal_zirc.add_nuclide("Sn115", 0.00493, percent_type="wo")
cal_zirc.add_nuclide("Sn116", 0.21083, percent_type="wo")
cal_zirc.add_nuclide("Sn117", 0.11136, percent_type="wo")
cal_zirc.add_nuclide("Sn118", 0.35119, percent_type="wo")
cal_zirc.add_nuclide("Sn119", 0.124555, percent_type="wo")
cal_zirc.add_nuclide("Sn120", 0.47241, percent_type="wo")
cal_zirc.add_nuclide("Sn122", 0.067135, percent_type="wo")
cal_zirc.add_nuclide("Sn124", 0.083955, percent_type="wo")
cal_zirc.add_nuclide("Fe54", 0.0122745, percent_type="wo")
cal_zirc.add_nuclide("Fe56", 0.1926834, percent_type="wo")
cal_zirc.add_nuclide("Fe57", 0.004431, percent_type="wo")
cal_zirc.add_nuclide("Fe58", 0.0005922, percent_type="wo")
cal_zirc.add_nuclide("Zr90", 50.16, percent_type="wo")
cal_zirc.add_nuclide("Zr91", 10.94, percent_type="wo")
cal_zirc.add_nuclide("Zr92", 16.72, percent_type="wo")
cal_zirc.add_nuclide("Zr94", 17.94, percent_type="wo")
cal_zirc.add_nuclide("Zr96", 2.73, percent_type="wo")

heavy_water = openmc.Material(material_id=7)
heavy_water.set_density("g/cc", 1.07801)
heavy_water.temperature = 1200
heavy_water.add_nuclide("O16", 1, percent_type="ao")
heavy_water.add_nuclide("H2", 1.996, percent_type="ao")
heavy_water.add_nuclide("H1", 0.004,percent_type="ao")
heavy_water.add_s_alpha_beta(name='c_H_in_H2O')
heavy_water.add_s_alpha_beta(name='c_D_in_D2O')

materials_file = openmc.Materials([fuel_21_hot, fuel_21_mid, fuel_21_low, fuel_21_dep, clad, water, tube, graphite, cal_zirc, heavy_water])
materials_file.export_to_xml(path="problem/materials.xml")

################## GEOMETRY ##################

surf1  = openmc.ZCylinder(x0=0.0,    y0=  4.384, r = 0.536)
surf2  = openmc.ZCylinder(x0=1.293,  y0=  4.189, r = 0.536)
surf3  = openmc.ZCylinder(x0=2.47,   y0=  3.622, r = 0.536)
surf4  = openmc.ZCylinder(x0=3.429,  y0=  2.732,  r = 0.536)
surf5  = openmc.ZCylinder(x0=4.082,  y0=  1.600,  r = 0.536)
surf6  = openmc.ZCylinder(x0=4.372,  y0=  0.325,  r = 0.536)
surf7  = openmc.ZCylinder(x0=4.273,  y0= -0.980,  r = 0.536)
surf8  = openmc.ZCylinder(x0=3.795,  y0= -2.200,  r = 0.536)
surf9  = openmc.ZCylinder(x0=2.979,  y0= -3.220,  r = 0.536)
surf10 = openmc.ZCylinder(x0=1.898,  y0= -3.950,  r = 0.536)
surf11 = openmc.ZCylinder(x0=0.648,  y0= -4.340,  r = 0.536)
surf12 = openmc.ZCylinder(x0=-1.293, y0=  4.189,  r = 0.536)
surf13 = openmc.ZCylinder(x0=-2.47,  y0=  3.622,  r = 0.536)
surf14 = openmc.ZCylinder(x0=-3.429, y0=  2.732,  r = 0.536)
surf15 = openmc.ZCylinder(x0=-4.082, y0=  1.600,  r = 0.536)
surf16 = openmc.ZCylinder(x0=-4.372, y0=  0.325,  r = 0.536)
surf17 = openmc.ZCylinder(x0=-4.273, y0= -0.980,  r = 0.536)
surf18 = openmc.ZCylinder(x0=-3.795, y0= -2.200,  r = 0.536)
surf19 = openmc.ZCylinder(x0=-2.979, y0= -3.220,  r = 0.536)
surf20 = openmc.ZCylinder(x0=-1.898, y0= -3.950,  r = 0.536)
surf21 = openmc.ZCylinder(x0=-0.648, y0= -4.340,  r = 0.536)
surf22 = openmc.ZCylinder(x0=0.684,  y0=  2.998, r = 0.536)
surf23 = openmc.ZCylinder(x0=1.918,  y0=  2.404, r = 0.536)
surf24 = openmc.ZCylinder(x0=2.771,  y0=  1.333, r = 0.536)
surf25 = openmc.ZCylinder(x0=3.075,  y0=  0.000, r = 0.536)
surf26 = openmc.ZCylinder(x0=2.77,   y0= -1.340, r = 0.536)
surf27 = openmc.ZCylinder(x0=1.915,  y0= -2.410, r = 0.536)
surf28 = openmc.ZCylinder(x0=0.681,  y0= -3.000, r = 0.536)
surf29 = openmc.ZCylinder(x0=-0.69,  y0= -3.000, r = 0.536)
surf30 = openmc.ZCylinder(x0=-1.92,  y0= -2.400, r = 0.536)
surf31 = openmc.ZCylinder(x0=-2.77,  y0= -1.330, r = 0.536)
surf32 = openmc.ZCylinder(x0=-3.07,  y0= -0.006, r = 0.536)
surf33 = openmc.ZCylinder(x0=-2.77,  y0=  1.340, r = 0.536)
surf34 = openmc.ZCylinder(x0=-1.91,  y0=  2.408, r = 0.536)
surf35 = openmc.ZCylinder(x0=-0.68,  y0=  3.000, r = 0.536)
surf36 = openmc.ZCylinder(x0=0.0,    y0=  1.730, r = 0.631)
surf37 = openmc.ZCylinder(x0=1.353,  y0=  1.078, r = 0.631)
surf38 = openmc.ZCylinder(x0=1.686,  y0= -0.390,  r = 0.631)
surf39 = openmc.ZCylinder(x0=0.749,  y0= -1.560,  r = 0.631)
surf40 = openmc.ZCylinder(x0=-0.75,  y0= -1.560,  r = 0.631)
surf41 = openmc.ZCylinder(x0=-1.69,  y0= -0.380,  r = 0.631)
surf42 = openmc.ZCylinder(x0=-1.35,  y0=  1.082,  r = 0.631)
surf43 = openmc.ZCylinder(x0=0.00,   y0=  0.00,  r = 0.631)
surf44 = openmc.ZCylinder(x0=0.0,    y0=  4.384,  r = 0.575)
surf45 = openmc.ZCylinder(x0=1.293,  y0=  4.189,  r = 0.575)
surf46 = openmc.ZCylinder(x0=2.47,   y0=  3.622,  r = 0.575)
surf47 = openmc.ZCylinder(x0=3.429,  y0=  2.732,  r = 0.575)
surf48 = openmc.ZCylinder(x0=4.082,  y0=  1.600,  r = 0.575)
surf49 = openmc.ZCylinder(x0=4.372,  y0=  0.325,  r = 0.575)
surf50 = openmc.ZCylinder(x0=4.273,  y0= -0.980,  r = 0.575)
surf51 = openmc.ZCylinder(x0=3.795,  y0= -2.200,  r = 0.575)
surf52 = openmc.ZCylinder(x0=2.979,  y0= -3.220,  r = 0.575)
surf53 = openmc.ZCylinder(x0=1.898,  y0= -3.950,  r = 0.575)
surf54 = openmc.ZCylinder(x0=0.648,  y0= -4.340,  r = 0.575)
surf55 = openmc.ZCylinder(x0=-1.293, y0=  4.189,  r = 0.575)
surf56 = openmc.ZCylinder(x0=-2.47,  y0=  3.622,  r = 0.575)
surf57 = openmc.ZCylinder(x0=-3.429, y0=  2.732,  r = 0.575)
surf58 = openmc.ZCylinder(x0=-4.082, y0=  1.600,  r = 0.575)
surf59 = openmc.ZCylinder(x0=-4.372, y0=  0.325,  r = 0.575)
surf60 = openmc.ZCylinder(x0=-4.273, y0= -0.980,  r = 0.575)
surf61 = openmc.ZCylinder(x0=-3.795, y0= -2.200,  r = 0.575)
surf62 = openmc.ZCylinder(x0=-2.979, y0= -3.220,  r = 0.575)
surf63 = openmc.ZCylinder(x0=-1.898, y0= -3.950,  r = 0.575)
surf64 = openmc.ZCylinder(x0=-0.648, y0= -4.340,  r = 0.575)
surf65 = openmc.ZCylinder(x0=0.684,  y0=  2.998, r = 0.575)
surf66 = openmc.ZCylinder(x0=1.918,  y0=  2.404, r = 0.575)
surf67 = openmc.ZCylinder(x0=2.771,  y0=  1.333, r = 0.575)
surf68 = openmc.ZCylinder(x0=3.075,  y0=  0.000, r = 0.575)
surf69 = openmc.ZCylinder(x0=2.77,   y0= -1.340, r = 0.575)
surf70 = openmc.ZCylinder(x0=1.915,  y0= -2.410, r = 0.575)
surf71 = openmc.ZCylinder(x0=0.681,  y0= -3.000, r = 0.575)
surf72 = openmc.ZCylinder(x0=-0.69,  y0= -3.000, r = 0.575)
surf73 = openmc.ZCylinder(x0=-1.92,  y0= -2.400, r = 0.575)
surf74 = openmc.ZCylinder(x0=-2.77,  y0= -1.330, r = 0.575)
surf75 = openmc.ZCylinder(x0=-3.07,  y0= -0.006, r = 0.575)
surf76 = openmc.ZCylinder(x0=-2.77,  y0=  1.340, r = 0.575)
surf77 = openmc.ZCylinder(x0=-1.91,  y0=  2.408, r = 0.575)
surf78 = openmc.ZCylinder(x0=-0.68,  y0=  3.000, r = 0.575)
surf79 = openmc.ZCylinder(x0=0.0,    y0=  1.730, r = 0.675)
surf80 = openmc.ZCylinder(x0=1.353,  y0=  1.078, r = 0.675)
surf81 = openmc.ZCylinder(x0=1.686,  y0= -0.390,  r = 0.675)
surf82 = openmc.ZCylinder(x0=0.749,  y0= -1.560,  r = 0.675)
surf83 = openmc.ZCylinder(x0=-0.75,  y0= -1.560,  r = 0.675)
surf84 = openmc.ZCylinder(x0=-1.69,  y0= -0.380,  r = 0.675)
surf85 = openmc.ZCylinder(x0=-1.35,  y0=  1.082,  r = 0.675)
surf86 = openmc.ZCylinder(x0=0.00,   y0=  0.0,  r = 0.675)
surf87 = openmc.ZCylinder(x0=0.00,   y0=  0.0,  r = 5.17)
surf88 = openmc.ZCylinder(x0=0.00,   y0=  0.0,  r = 5.8189)
surf89 = openmc.ZCylinder(x0=0.00,   y0=  0.0,  r = 7.55)
surf90 = openmc.ZCylinder(x0=0.00,   y0=  0.0,  r = 7.8)

surf91 = openmc.XPlane(x0=-11, boundary_type="reflective")
surf92 = openmc.XPlane(x0=11, boundary_type="reflective")
surf93 = openmc.YPlane(y0=-11, boundary_type="reflective")
surf94 = openmc.YPlane(y0=11, boundary_type="reflective")
surf95 = openmc.ZPlane(z0=-24.765, boundary_type="reflective")
surf96 = openmc.ZPlane(z0=24.765, boundary_type="reflective")

univ1 = openmc.Universe(universe_id=1)
cell1  = openmc.Cell(cell_id=1,  fill=fuel_21_hot, region=-surf1 )
cell2  = openmc.Cell(cell_id=2,  fill=fuel_21_hot, region=-surf2 )
cell3  = openmc.Cell(cell_id=3,  fill=fuel_21_hot, region=-surf3 )
cell4  = openmc.Cell(cell_id=4,  fill=fuel_21_hot, region=-surf4 )
cell5  = openmc.Cell(cell_id=5,  fill=fuel_21_hot, region=-surf5 )
cell6  = openmc.Cell(cell_id=6,  fill=fuel_21_hot, region=-surf6 )
cell7  = openmc.Cell(cell_id=7,  fill=fuel_21_hot, region=-surf7 )
cell8  = openmc.Cell(cell_id=8,  fill=fuel_21_hot, region=-surf8 )
cell9  = openmc.Cell(cell_id=9,  fill=fuel_21_hot, region=-surf9 )
cell10 = openmc.Cell(cell_id=10, fill=fuel_21_hot, region=-surf10)
cell11 = openmc.Cell(cell_id=11, fill=fuel_21_hot, region=-surf11)
cell12 = openmc.Cell(cell_id=12, fill=fuel_21_hot, region=-surf12)
cell13 = openmc.Cell(cell_id=13, fill=fuel_21_hot, region=-surf13)
cell14 = openmc.Cell(cell_id=14, fill=fuel_21_hot, region=-surf14)
cell15 = openmc.Cell(cell_id=15, fill=fuel_21_hot, region=-surf15)
cell16 = openmc.Cell(cell_id=16, fill=fuel_21_hot, region=-surf16)
cell17 = openmc.Cell(cell_id=17, fill=fuel_21_hot, region=-surf17)
cell18 = openmc.Cell(cell_id=18, fill=fuel_21_hot, region=-surf18)
cell19 = openmc.Cell(cell_id=19, fill=fuel_21_hot, region=-surf19)
cell20 = openmc.Cell(cell_id=20, fill=fuel_21_hot, region=-surf20)
cell21 = openmc.Cell(cell_id=21, fill=fuel_21_hot, region=-surf21)
cell22 = openmc.Cell(cell_id=22, fill=fuel_21_mid, region=-surf22)
cell23 = openmc.Cell(cell_id=23, fill=fuel_21_mid, region=-surf23)
cell24 = openmc.Cell(cell_id=24, fill=fuel_21_mid, region=-surf24)
cell25 = openmc.Cell(cell_id=25, fill=fuel_21_mid, region=-surf25)
cell26 = openmc.Cell(cell_id=26, fill=fuel_21_mid, region=-surf26)
cell27 = openmc.Cell(cell_id=27, fill=fuel_21_mid, region=-surf27)
cell28 = openmc.Cell(cell_id=28, fill=fuel_21_mid, region=-surf28)
cell29 = openmc.Cell(cell_id=29, fill=fuel_21_mid, region=-surf29)
cell30 = openmc.Cell(cell_id=30, fill=fuel_21_mid, region=-surf30)
cell31 = openmc.Cell(cell_id=31, fill=fuel_21_mid, region=-surf31)
cell32 = openmc.Cell(cell_id=32, fill=fuel_21_mid, region=-surf32)
cell33 = openmc.Cell(cell_id=33, fill=fuel_21_mid, region=-surf33)
cell34 = openmc.Cell(cell_id=34, fill=fuel_21_mid, region=-surf34)
cell35 = openmc.Cell(cell_id=35, fill=fuel_21_mid, region=-surf35)
cell36 = openmc.Cell(cell_id=36, fill=fuel_21_low, region=-surf36)
cell37 = openmc.Cell(cell_id=37, fill=fuel_21_low, region=-surf37)
cell38 = openmc.Cell(cell_id=38, fill=fuel_21_low, region=-surf38)
cell39 = openmc.Cell(cell_id=39, fill=fuel_21_low, region=-surf39)
cell40 = openmc.Cell(cell_id=40, fill=fuel_21_low, region=-surf40)
cell41 = openmc.Cell(cell_id=41, fill=fuel_21_low, region=-surf41)
cell42 = openmc.Cell(cell_id=42, fill=fuel_21_low, region=-surf42)
cell43 = openmc.Cell(cell_id=43, fill=fuel_21_dep, region=-surf43)
cell44 = openmc.Cell(cell_id=44, fill=clad, region=+surf1 & -surf44)
cell45 = openmc.Cell(cell_id=45, fill=clad, region=+surf2 & -surf45)
cell46 = openmc.Cell(cell_id=46, fill=clad, region=+surf3 & -surf46)
cell47 = openmc.Cell(cell_id=47, fill=clad, region=+surf4 & -surf47)
cell48 = openmc.Cell(cell_id=48, fill=clad, region=+surf5 & -surf48)
cell49 = openmc.Cell(cell_id=49, fill=clad, region=+surf6 & -surf49)
cell50 = openmc.Cell(cell_id=50, fill=clad, region=+surf7 & -surf50)
cell51 = openmc.Cell(cell_id=51, fill=clad, region=+surf8 & -surf51)
cell52 = openmc.Cell(cell_id=52, fill=clad, region=+surf9 & -surf52)
cell53 = openmc.Cell(cell_id=53, fill=clad, region=+surf10 & -surf53)
cell54 = openmc.Cell(cell_id=54, fill=clad, region=+surf11 & -surf54)
cell55 = openmc.Cell(cell_id=55, fill=clad, region=+surf12 & -surf55)
cell56 = openmc.Cell(cell_id=56, fill=clad, region=+surf13 & -surf56)
cell57 = openmc.Cell(cell_id=57, fill=clad, region=+surf14 & -surf57)
cell58 = openmc.Cell(cell_id=58, fill=clad, region=+surf15 & -surf58)
cell59 = openmc.Cell(cell_id=59, fill=clad, region=+surf16 & -surf59)
cell60 = openmc.Cell(cell_id=60, fill=clad, region=+surf17 & -surf60)
cell61 = openmc.Cell(cell_id=61, fill=clad, region=+surf18 & -surf61)
cell62 = openmc.Cell(cell_id=62, fill=clad, region=+surf19 & -surf62)
cell63 = openmc.Cell(cell_id=63, fill=clad, region=+surf20 & -surf63)
cell64 = openmc.Cell(cell_id=64, fill=clad, region=+surf21 & -surf64)
cell65 = openmc.Cell(cell_id=65, fill=clad, region=+surf22 & -surf65)
cell66 = openmc.Cell(cell_id=66, fill=clad, region=+surf23 & -surf66)
cell67 = openmc.Cell(cell_id=67, fill=clad, region=+surf24 & -surf67)
cell68 = openmc.Cell(cell_id=68, fill=clad, region=+surf25 & -surf68)
cell69 = openmc.Cell(cell_id=69, fill=clad, region=+surf26 & -surf69)
cell70 = openmc.Cell(cell_id=70, fill=clad, region=+surf27 & -surf70)
cell71 = openmc.Cell(cell_id=71, fill=clad, region=+surf28 & -surf71)
cell72 = openmc.Cell(cell_id=72, fill=clad, region=+surf29 & -surf72)
cell73 = openmc.Cell(cell_id=73, fill=clad, region=+surf30 & -surf73)
cell74 = openmc.Cell(cell_id=74, fill=clad, region=+surf31 & -surf74)
cell75 = openmc.Cell(cell_id=75, fill=clad, region=+surf32 & -surf75)
cell76 = openmc.Cell(cell_id=76, fill=clad, region=+surf33 & -surf76)
cell77 = openmc.Cell(cell_id=77, fill=clad, region=+surf34 & -surf77)
cell78 = openmc.Cell(cell_id=78, fill=clad, region=+surf35 & -surf78)
cell79 = openmc.Cell(cell_id=79, fill=clad, region=+surf36 & -surf79)
cell80 = openmc.Cell(cell_id=80, fill=clad, region=+surf37 & -surf80)
cell81 = openmc.Cell(cell_id=81, fill=clad, region=+surf38 & -surf81)
cell82 = openmc.Cell(cell_id=82, fill=clad, region=+surf39 & -surf82)
cell83 = openmc.Cell(cell_id=83, fill=clad, region=+surf40 & -surf83)
cell84 = openmc.Cell(cell_id=84, fill=clad, region=+surf41 & -surf84)
cell85 = openmc.Cell(cell_id=85, fill=clad, region=+surf42 & -surf85)
cell86 = openmc.Cell(cell_id=86, fill=clad, region=+surf43 & -surf86)

cell87 = openmc.Cell(cell_id=87, fill=water, region=-surf87 & +surf44 & +surf43 & +surf45 & +surf46 & +surf47 & +surf48 & +surf49 & +surf50 & +surf51 & +surf52 & +surf53 & +surf54 & +surf55 & +surf56 & +surf57 & +surf58 & +surf59 & +surf60 & +surf61 & +surf62 & +surf63 & +surf64 & +surf65 & +surf66 & +surf67 & +surf68 & +surf69 & +surf70 & +surf71 & +surf72 & +surf73 & +surf74 & +surf75 & +surf76 & +surf77 & +surf78 & +surf79 & +surf80 & +surf81 & +surf82 & +surf83 & +surf84 & +surf85 & +surf86)

cell88 = openmc.Cell(cell_id=88, fill=tube, region=+surf87 & -surf88)
cell89 = openmc.Cell(cell_id=89, fill=graphite, region=-surf89 & +surf88)
cell90 = openmc.Cell(cell_id=90, fill=cal_zirc, region=+surf89)

univ1.add_cells([cell1, cell2, cell3, cell4, cell5, cell6, cell7, cell8, cell9, cell10, cell11, cell12, cell13, cell14, cell15, cell16, cell17, cell18, cell19, cell20, cell21, cell22, cell23, cell24, cell25, cell26, cell27, cell28, cell29, cell30, cell31, cell32, cell33, cell34, cell35, cell36, cell37, cell38, cell39, cell40, cell41, cell42, cell43, cell44, cell45, cell46, cell47, cell48, cell49, cell50, cell51, cell52, cell53, cell54, cell55, cell56, cell57, cell58, cell59, cell60, cell61, cell62, cell63, cell64, cell65, cell66, cell67, cell68, cell69, cell70, cell71, cell72, cell73, cell74, cell75, cell76, cell77, cell78, cell79, cell80, cell81, cell82, cell83, cell84, cell85, cell86, cell87, cell88, cell89, cell90])

univ2 = openmc.Universe()
cell91 = openmc.Cell(fill=univ1, region=-surf90)
cell92 = openmc.Cell(cell_id=92, fill=heavy_water, region=+surf90)
univ2.add_cells([cell91, cell92])

cell93 = openmc.Cell(cell_id=93, fill=univ2, region=+surf91 & -surf92 & +surf93 & -surf94 & +surf95 & -surf96)
univ3 = openmc.Universe(cells=[cell93])

geometry_file = openmc.Geometry(univ3)
geometry_file.export_to_xml(path="problem/geometry.xml")

################## TALLIES ##################

filter1 = openmc.CellFilter(bins=[93, 87])
filter2 = openmc.CellFilter(bins=[93, 87])
filter3 = openmc.CellFilter(bins=[93, 87])
filter4 = openmc.CellFilter(bins=[93, 87])
filter5 = openmc.CellFilter(bins=[93, 87])
filter6 = openmc.CellFilter(bins=[93, 87])
filter7 = openmc.CellFilter(bins=[93, 87])
filter8 = openmc.DelayedGroupFilter(bins=[1, 2, 3, 4, 5, 6])

tally1 = openmc.Tally()
tally1.filters = [filter1]
tally1.scores = ["flux"]

tally2 = openmc.Tally()
tally2.filters = [filter1]
tally2.scores = ["nu-fission"]

tally3 = openmc.Tally()
tally3.filters = [filter1]
tally3.scores = ["delayed-nu-fission"]

tally4 = openmc.Tally()
tally4.filters = [filter1]
tally4.scores = ["fission"]

tally5 = openmc.Tally()
tally5.filters = [filter1]
tally5.scores = ["absorption"]

tally6 = openmc.Tally()
tally6.filters = [filter1]
tally6.scores = ["scatter"]

tally7 = openmc.Tally()
tally7.filters = [filter1]
tally7.scores = ["total"]

tally8 = openmc.Tally()
tally8.filters = [filter8]
tally8.scores = ["delayed-nu-fission"]

tallies_file = openmc.Tallies([tally1, tally2, tally3, tally4, tally5, tally6, tally7, tally8])
tallies_file.export_to_xml(path="problem/tallies.xml")

################## PLOTS ##################

plot = openmc.Plot()
plot.pixels = (1000, 1000)
plot.width = (25.0, 25.0)
# plot._colors = (123, 123, 231)

plot_file = openmc.Plots([plot])
plot_file.export_to_xml(path="problem/plots.xml")

################## SETTINGS ##################

settings_file = openmc.Settings()
settings_file.temperature = {'method': 'interpolation',
                             'tolerance': 550}
settings_file.source = openmc.IndependentSource(space=openmc.stats.Box(
                        (-11, -11, -24.765), ( 11,  11,  24.765)))
settings_file.batches = 600
settings_file.inactive = 50
settings_file.particles = 10000

settings_file.export_to_xml(path="problem/settings.xml")

openmc.plot_geometry(cwd='problem')
openmc.run(cwd='problem')