diff --git a/WC_Layers/OpenMC_PhotonTransport.py b/WC_Layers/OpenMC_PhotonTransport.py
new file mode 100644
index 0000000..a1bdab4
--- /dev/null
+++ b/WC_Layers/OpenMC_PhotonTransport.py
@@ -0,0 +1,84 @@
+import openmc
+import numpy as np
+import yaml
+import argparse
+from Source_Mesh_Reader import extract_source_data
+from TwoLayers_Materials import alara_element_densities, make_materials
+
+def make_source(bounds, cells, mesh_file, source_mesh_index):
+      '''
+    Creates a list of OpenMC sources, complete with the relevant space and energy distributions
+    
+    inputs:
+        bounds : iterable of photon energy bounds (float)
+        cells: list of OpenMC Cell objects
+        mesh_file: .h5/.h5m mesh onto which photon source will be distributed
+        source_mesh_index: index specifying the photon source from which data is extracted
+        
+    output:
+        source_list: list of OpenMC independent sources
+        total_mesh: OpenMC Unstructured Mesh object
+    '''
+    source_list = []
+    total_mesh = openmc.UnstructuredMesh(mesh_file, library='moab')
+    for index, (lower_bound, upper_bound) in enumerate(zip(bounds[:-1],bounds[1:])):
+        mesh_dist = openmc.stats.MeshSpatial(total_mesh, strengths=esd[source_mesh_index][:,index], volume_normalized=False)
+        energy_dist = openmc.stats.Uniform(a=lower_bound, b=upper_bound)
+        source_list.append(openmc.IndependentSource(space=mesh_dist, energy=energy_dist, strength=np.sum(esd[source_mesh_index][:, index]), particle='photon', domains=cells))
+    return source_list, total_mesh
+
+def make_tallies(total_mesh, tallied_cells):
+      '''
+    Creates tallies and assigns energy, spatial, and particle filters.
+    
+    inputs: 
+        total_mesh: OpenMC unstructured mesh object
+        tallied_cells: OpenMC Cell/iterable of OpenMC Cell objects/iterable of Cell ID #
+        
+    outputs:
+        talls: OpenMC Tallies object
+    '''
+    particle_filter = openmc.ParticleFilter('photon')
+    total_filter = openmc.MeshFilter(total_mesh)
+    
+    photon_tally = openmc.Tally(tally_id=1, name="Photon tally")
+    photon_tally.scores = ['flux', 'elastic', 'absorption']
+    
+    # Implementing filter for neutron tally through shells with material
+    cell_filter = openmc.CellFilter(tallied_cells)
+    photon_tally.filters = [cell_filter, total_filter]
+
+    # Obtain coefficients to calculate effective dose
+    dose_energy, dose = openmc.data.dose_coefficients('photon', 'AP')
+    dose_filter = openmc.EnergyFunctionFilter(dose_energy, dose)
+
+    # Vitamin-J energy filter created to assign to the flux tally.
+    energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-42")
+
+    spectrum_tally = openmc.Tally(tally_id=2, name="Flux spectrum")
+    spectrum_tally.filters = [cell_filter, total_filter, energy_filter_flux, particle_filter, dose_filter]
+    spectrum_tally.scores = ['flux']
+    
+    talls = openmc.Tallies([photon_tally, spectrum_tally])
+    return talls
+  
+def make_settings(source_list, total_batches, inactive_batches, num_particles, run_mode):
+      '''
+    Creates an OpenMC Settings object
+    
+    inputs:
+        source_list: iterable of OpenMC SourceBase objects
+    outputs:
+        sets: OpenMC Settings object
+    '''
+    sets = openmc.Settings()
+    sets.batches = total_batches
+    sets.inactive = inactive_batches
+    sets.particles = num_particles
+    sets.source = source_list
+    sets.run_mode = run_mode
+    return sets
+
+def create_openmc_model(geom_object, mats_object, sets_object, talls_object):
+    model = openmc.model.Model(geometry = geom_object, materials = mats_object, settings = sets_object, tallies = talls_object)
+    return model
diff --git a/WC_Layers/PhotonTransport_Inputs.yaml b/WC_Layers/PhotonTransport_Inputs.yaml
new file mode 100644
index 0000000..40ac0fa
--- /dev/null
+++ b/WC_Layers/PhotonTransport_Inputs.yaml
@@ -0,0 +1,57 @@
+filename_dict :
+    alara_el_lib : elelib.std
+    mesh_file_openmc : OpenMC_Mesh.h5m
+
+file_indices: 
+    source_mesh_index : 0
+
+source_info :
+    phtn_e_bounds :
+        - 0
+        - 1.00e+4
+        - 2.00e+4
+        - 5.00e+4
+        - 1.00e+5 
+        - 2.00e+5 
+        - 3.00e+5 
+        - 4.00e+5 
+        - 6.00e+5 
+        - 8.00e+5 
+        - 1.00e+6 
+        - 1.22e+6 
+        - 1.44e+6 
+        - 1.66e+6 
+        - 2.00e+6 
+        - 2.50e+6 
+        - 3.00e+6 
+        - 4.00e+6 
+        - 5.00e+6 
+        - 6.50e+6 
+        - 8.00e+6 
+        - 1.00e+7 
+        - 1.20e+7 
+        - 1.40e+7 
+        - 2.00e+7
+    
+mat_info :
+    element_list : #add in order of radially innermost to outermost
+        - W
+        - C
+
+geom_info :
+     inner_radius : 995
+     thicknesses : #corresponds to each of the elements in element_list
+         - 5
+         - 5
+     outer_boundary_type : vacuum
+
+tally_info :
+     tallied_elements : #change according to desired tally region/material
+         - W
+         - C
+
+settings_info :
+     particles : 100000
+     batches : 10
+     inactive_batches : 1
+     run_mode : fixed source
diff --git a/WC_Layers/Photon_TallytoVtk.py b/WC_Layers/Photon_TallytoVtk.py
new file mode 100644
index 0000000..8a77639
--- /dev/null
+++ b/WC_Layers/Photon_TallytoVtk.py
@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct  3 23:18:20 2024
+
+@author: Anupama Rajendra
+"""
+import openmc
+import matplotlib.pyplot as plt
+
+def read_statepoint(sp_filename, flux_spectrum_tally_id, mesh_number, summed_axes_energy_filter, energy_filter_index):
+    with openmc.StatePoint(sp_filename) as sp:
+        flux_spectrum = sp.get_tally(id=flux_spectrum_tally_id) 
+        mesh=sp.meshes[mesh_number]
+        # Return tally data condensed into 1 dimension per filter
+        tally_data_reshaped = flux_spectrum.get_reshaped_data(value='mean')
+        flux_sum_energy = tally_data_reshaped.sum(axis=summed_axes_energy_filter)   
+        #Vitamin-J energy filter:
+        e_filter = flux_spectrum.filters[energy_filter_index]
+        #Lower bounds of the energy bins
+        e_filter_lower = e_filter.bins[:, 0]
+    return e_filter_lower, flux_sum_energy, tally_data_reshaped, mesh
+
+def plot_flux_data(e_filter_lower, flux_sum_energy, figure_filename):
+    # Plot flux spectrum
+    fix, ax = plt.subplots()
+    ax.loglog(e_filter_lower, flux_sum_energy, drawstyle='steps-post')
+    ax.set_xlabel('Energy [eV]')
+    ax.set_ylabel('Flux [photon-cm/source]')
+    ax.grid(True, which='both')
+    plt.savefig(figure_filename)
+    plt.show()
+
+def write_to_vtk(tally_data_reshaped, summed_axes_mesh_filter, vtk_filename, mesh):
+    mesh_data = tally_data_reshaped.sum(axis=summed_axes_mesh_filter)
+    mesh.write_data_to_vtk(filename=vtk_filename, datasets={"mean":mesh_data.flatten()})
diff --git a/WC_Layers/Source_Mesh_Reader.py b/WC_Layers/Source_Mesh_Reader.py
new file mode 100644
index 0000000..5b3d612
--- /dev/null
+++ b/WC_Layers/Source_Mesh_Reader.py
@@ -0,0 +1,51 @@
+import h5py
+import numpy as np
+import yaml
+import argparse
+
+def extract_source_data(source_mesh_list, num_elements, photon_groups):
+    '''
+    Identifies the location of the source density dataset within each mesh file.
+    
+    input: iterable of .h5m filenames (str), whose files contain photon source information
+    output: numpy array of source density data with rows = # of mesh elements and 
+    columns = # number of photon groups, with one array per source mesh
+    '''
+    
+    sd_list = np.ndarray((len(source_mesh_list), num_elements, photon_groups))
+    for source_index, source_name in enumerate(source_mesh_list):
+         file = h5py.File(source_name, 'r')
+         sd_list[source_index,:] = file['tstt']['elements']['Tet4']['tags']['source_density'][:]
+    return sd_list
+
+def save_source_density(sd_list, sd_filename):
+    '''
+    Saves source density data as a text file.
+    
+    sd_list: list of source density datasets from h5m files
+    sd_filename: user-provided filename that appears as a prefix for the text files (str)
+    '''
+    for sd_index, source_density in enumerate(sd_list):
+        with open(f'{sd_filename}_{sd_index + 1}.txt', 'w') as source:
+            for tet_element in source_density:
+                source.write(' '.join(map(str, tet_element)) + '\n')
+                
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--Mesh_Reader_YAML', default = 'Source_Mesh_Reader_Inputs.yaml', help="Path (str) to YAML file containing required inputs for Source_Mesh_Reader")
+    args = parser.parse_args()
+    smr_yaml = args.Mesh_Reader_YAML
+    
+    with open(smr_yaml, 'r') as yaml_file:
+        inputs = yaml.safe_load(yaml_file)
+        
+    source_mesh_list = inputs['source_meshes']
+    num_elements = inputs['num_elements']
+    photon_groups = inputs['photon_groups']
+    sd_filename = inputs['sd_filename']
+    
+    esd = extract_source_data(source_mesh_list, num_elements, photon_groups)
+    ssd = save_source_density(esd, sd_filename)
+
+if __name__ == "__main__":
+    main()
diff --git a/WC_Layers/Source_Mesh_Reader_Inputs.yaml b/WC_Layers/Source_Mesh_Reader_Inputs.yaml
new file mode 100644
index 0000000..b4780b1
--- /dev/null
+++ b/WC_Layers/Source_Mesh_Reader_Inputs.yaml
@@ -0,0 +1,6 @@
+num_elements : 6168
+photon_groups : 24
+source_meshes :
+    - source_mesh_1.h5m
+    - source_mesh_2.h5m
+sd_filename : source_density    
diff --git a/WC_Layers/TwoLayers_Geometry.py b/WC_Layers/TwoLayers_Geometry.py
new file mode 100644
index 0000000..295fe85
--- /dev/null
+++ b/WC_Layers/TwoLayers_Geometry.py
@@ -0,0 +1,23 @@
+import openmc
+
+def make_spherical_shells(materials, thicknesses, inner_radius, outer_boundary_type):    
+    '''
+    Creates a set of concentric spherical shells, each with its own material & inner/outer radius.
+    inputs:
+        materials : iterable of OpenMC Material objects
+        thicknesses: thickness of each OpenMC Material
+        inner_radius: the radius of the innermost spherical shell
+    '''
+    layers = zip(materials, thicknesses)
+    inner_sphere = openmc.Sphere(r = inner_radius)
+    cells = [openmc.Cell(fill = None, region = -inner_sphere)]
+    for (material, thickness) in layers:
+        outer_radius = inner_radius + thickness
+        outer_sphere = openmc.Sphere(r = outer_radius)
+        cells.append(openmc.Cell(fill = material, region = +inner_sphere & -outer_sphere))
+        outer_radius = inner_radius
+        outer_sphere = inner_sphere
+    outer_sphere.boundary_type = outer_boundary_type     
+    cells.append(openmc.Cell(fill = None, region = +outer_sphere)) 
+    geometry = openmc.Geometry(cells)    
+    return geometry
diff --git a/WC_Layers/TwoLayers_Materials.py b/WC_Layers/TwoLayers_Materials.py
new file mode 100644
index 0000000..cc151dc
--- /dev/null
+++ b/WC_Layers/TwoLayers_Materials.py
@@ -0,0 +1,35 @@
+import openmc
+
+def alara_element_densities(alara_fp):
+    '''
+    Creates a dictionary where keys = element names (str) and values = element density (float)
+    inputs:
+        alara_filepath : path to ALARA element library
+    '''
+    with open(alara_fp) as ALARA_Lib:
+        libLines = ALARA_Lib.readlines()
+    num_lines = len(libLines)
+    density_dict = {}
+    line_num = 0
+    while line_num < num_lines:
+        element_data = libLines[line_num].strip().split()
+        element_name = element_data[0].lower()
+        density_dict[element_name] = float(element_data[3])
+        line_num += int(element_data[4]) + 1
+    return density_dict
+    
+def make_materials(elements, density_dict):
+    '''
+    Creates an OpenMC Materials object using user-specified elements
+    inputs:
+        elements: iterable of element names (str)
+        density_dict: dictionary with keys = element names (str) and values = element density (float)
+    '''    
+    mats = openmc.Materials([])
+    for element_id, element in enumerate(elements):
+        mat = openmc.Material(material_id=element_id+1, name=element)
+        mat.add_element(element, 1.00)
+        mat.set_density('g/cm3', density_dict.get(element.lower()))
+        mats.append(mat)
+    return mats
+