diff --git a/romanisim/image.py b/romanisim/image.py
index 1f5e9cda..ba32a3ae 100644
--- a/romanisim/image.py
+++ b/romanisim/image.py
@@ -29,10 +29,7 @@
 from romanisim import log
 
 import roman_datamodels
-try:
-    import roman_datamodels.maker_utils as maker_utils
-except ImportError:
-    import roman_datamodels.testing.utils as maker_utils
+import roman_datamodels.maker_utils as maker_utils
 
 
 # galsim fluxes are in photons / cm^2 / s
@@ -217,7 +214,8 @@ def trim_objlist(objlist, image):
 
 
 def add_objects_to_image(image, objlist, xpos, ypos, psf,
-                         flux_to_counts_factor, bandpass=None, filter_name=None,
+                         flux_to_counts_factor, convtimes=None,
+                         bandpass=None, filter_name=None,
                          rng=None, seed=None):
     """Add sources to an image.
 
@@ -227,17 +225,19 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
     Parameters
     ----------
     image : galsim.Image
-        Image to which to add sources with associated WCS.
+        Image to which to add sources with associated WCS. Updated in place.
     objlist : list[CatalogObject]
         Objects to add to image
     xpos, ypos : array_like
         x & y positions of sources (pixel) at which sources should be added
     psf : galsim.Profile
         PSF for image
-    flux_to_counts_factor : float
+    flux_to_counts_factor : float or list
         physical fluxes in objlist (whether in profile SEDs or flux arrays)
         should be multiplied by this factor to convert to total counts in the
         image
+    convtimes: array_like
+        Exposure times with unit scaling to convert to output rate units
     bandpass : galsim.Bandpass
         bandpass in which image is being rendered.  This is used only in cases
         where chromatic profiles & PSFs are being used.
@@ -279,6 +279,7 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
     for i, obj in enumerate(objlist):
         t0 = time.time()
         image_pos = galsim.PositionD(xpos[i], ypos[i])
+        pwcs = image.wcs.local(image_pos)
         profile = obj.profile
         if not chromatic:
             if obj.flux is None:
@@ -289,18 +290,22 @@ def add_objects_to_image(image, objlist, xpos, ypos, psf,
             psf0 = psf.at_position(xpos[i], ypos[i])
         else:
             psf0 = psf
-        final = galsim.Convolve(profile * flux_to_counts_factor, psf0)
+        factor = flux_to_counts_factor[i] if isinstance(flux_to_counts_factor, list) else flux_to_counts_factor
+        final = galsim.Convolve(profile * factor, psf0)
         if chromatic:
             stamp = final.drawImage(
-                bandpass, center=image_pos, wcs=image.wcs.local(image_pos),
+                bandpass, center=image_pos, wcs=pwcs,
                 method='phot', rng=rng)
         else:
             try:
                 stamp = final.drawImage(center=image_pos,
-                                        wcs=image.wcs.local(image_pos))
+                                        wcs=pwcs)
             except galsim.GalSimFFTSizeError:
                 log.warning(f'Skipping source {i} due to too '
                             f'large FFT needed for desired accuracy.')
+        if convtimes is not None:
+            stamp /= convtimes[i]
+
         bounds = stamp.bounds & image.bounds
         if bounds.area() > 0:
             image[bounds] += stamp[bounds]
diff --git a/romanisim/l3.py b/romanisim/l3.py
index 0b47022c..4ee19ed3 100644
--- a/romanisim/l3.py
+++ b/romanisim/l3.py
@@ -1,80 +1,104 @@
-"""Function for Level 3-like images.
+"""Roman WFI simulator functions for Level 3 mosaics.
 
+Based on galsim's implementation of Roman image simulation.  Uses galsim Roman modules
+for most of the real work.
 """
 
 import numpy as np
+import math
+import astropy.time
+from astropy import table
 import galsim
-from . import log
-from . import image, wcs, psf
-
-
-def add_objects_to_l3(l3_mos, source_cat, rng=None, seed=None):
+from galsim import roman
+
+from . import parameters
+from . import util
+import romanisim.wcs
+import romanisim.l1
+import romanisim.bandpass
+import romanisim.psf
+import romanisim.image
+import romanisim.persistence
+from romanisim import log
+import roman_datamodels.maker_utils as maker_utils
+import roman_datamodels.datamodels as rdm
+from roman_datamodels.stnode import WfiMosaic
+
+
+def add_objects_to_l3(l3_mos, source_cat, exptimes, xpos=None, ypos=None, coords=None, cps_conv=1.0, unit_factor=1.0,
+                      filter_name=None, coords_unit='rad', wcs=None, psf=None, rng=None, seed=None):
     """Add objects to a Level 3 mosaic
 
     Parameters
     ----------
-    l3_mos : MosaicModel
+    l3_mos : MosaicModel or galsim.Image
         Mosaic of images
     source_cat : list
         List of catalog objects to add to l3_mos
+    exptimes : list
+        Exposure times to scale back to rate units
+    xpos, ypos : array_like
+        x & y positions of sources (pixel) at which sources should be added
+    coords : array_like
+        ra & dec positions of sources (coords_unit) at which sources should be added
+    cps_conv : float
+        Factor to convert data to cps
+    unit_factor: float
+        Factor to convert counts data to MJy / sr
+    coords_unit : string
+        units of coords
+    wcs : galsim.GSFitsWCS
+        WCS corresponding to image
+    psf : galsim.Profile
+        PSF for image
+    rng : galsim.BaseDeviate
+        random number generator to use
+    seed : int
+        seed to use for random number generator
 
     Returns
     -------
     None
         l3_mos is updated in place
     """
-
-    if rng is None and seed is None:
-        seed = 143
-        log.warning(
-            'No RNG set, constructing a new default RNG from default seed.')
-    if rng is None:
-        rng = galsim.UniformDeviate(seed)
-
     # Obtain optical element
-    filter_name = l3_mos.meta.basic.optical_element
-
-    # Generate WCS
-    twcs = wcs.get_mosaic_wcs(l3_mos.meta)
-
-    # Create PSF
-    l3_psf = psf.make_psf(filter_name=filter_name, sca=2, chromatic=False, webbpsf=True)
-
-    # Generate x,y positions for sources
-    coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad]
-                       for o in source_cat])
-    sourcecountsall = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0)
-    xpos, ypos = sourcecountsall.wcs._xy(coords[:, 0], coords[:, 1])
-    xpos_idx = [round(x) for x in xpos]
-    ypos_idx = [round(y) for y in ypos]
-
-    # Create overall scaling factor map
-    Ct_all = (l3_mos.data.value / l3_mos.var_poisson)
-
-    # Cycle over sources and add them to the mosaic
-    for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)):
-        # Set scaling factor for injected sources
-        # Flux / sigma_p^2
-        Ct = (l3_mos.data[x][y].value / l3_mos.var_poisson[x][y].value)
-
-        # Create empty postage stamp galsim source image
-        sourcecounts = galsim.ImageF(l3_mos.data.shape[0], l3_mos.data.shape[1], wcs=twcs, xmin=0, ymin=0)
-
-        # Simulate source postage stamp
-        image.add_objects_to_image(sourcecounts, [source_cat[idx]], xpos=[xpos[idx]], ypos=[ypos[idx]],
-                                   psf=l3_psf, flux_to_counts_factor=Ct, bandpass=[filter_name],
-                                   filter_name=filter_name, rng=rng)
-
-        # Scale the source image back by its flux ratios
-        sourcecounts /= Ct
-
-        # Add sources to the original mosaic data array
-        l3_mos.data = (l3_mos.data.value + sourcecounts.array) * l3_mos.data.unit
-
-        # Note for the future - other noise sources (read and flat) need to be implemented
-
-    # Set new poisson variance
-    l3_mos.var_poisson = l3_mos.data.value / Ct_all
+    if filter_name is None:
+        filter_name = l3_mos.meta.basic.optical_element
+
+    # Generate WCS (if needed)
+    if wcs is None:
+        wcs = romanisim.wcs.get_mosaic_wcs(l3_mos.meta, shape=l3_mos.data.shape)
+
+    # Create PSF (if needed)
+    if psf is None:
+        psf = romanisim.psf.make_psf(filter_name=filter_name, sca=parameters.default_sca, chromatic=False, webbpsf=True)
+
+    # Create Image canvas to add objects to
+    if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)):
+        sourcecountsall = galsim.ImageF(l3_mos.data.value, wcs=wcs, xmin=0, ymin=0)
+    else:
+        sourcecountsall = l3_mos
+
+    # Create position arrays (if needed)
+    if any(pos is None for pos in [xpos, ypos]):
+        # Create coordinates (if needed)
+        if coords is None:
+            coords = np.array([[o.sky_pos.ra.rad, o.sky_pos.dec.rad]
+                               for o in source_cat])
+            coords_unit = 'rad'
+        # Generate x,y positions for sources
+        xpos, ypos = sourcecountsall.wcs.radecToxy(coords[:, 0], coords[:, 1], coords_unit)
+
+    # Add sources to the original mosaic data array
+    romanisim.image.add_objects_to_image(sourcecountsall, source_cat, xpos=xpos, ypos=ypos,
+                                         psf=psf, flux_to_counts_factor=[xpt * cps_conv for xpt in exptimes],
+                                         convtimes=[xpt / unit_factor for xpt in exptimes],
+                                         bandpass=[filter_name], filter_name=filter_name,
+                                         rng=rng, seed=seed)
+
+    # Save array with added sources
+    if isinstance(l3_mos, (rdm.MosaicModel, WfiMosaic)):
+        l3_mos.data = sourcecountsall.array * l3_mos.data.unit
 
     # l3_mos is updated in place, so no return
     return None
diff --git a/romanisim/parameters.py b/romanisim/parameters.py
index bbf79227..6da7eda6 100644
--- a/romanisim/parameters.py
+++ b/romanisim/parameters.py
@@ -84,6 +84,17 @@
     "linearity": None,
     "readnoise": 5.0 * u.DN,
     "saturation": 55000 * u.DN,
+    # Taken from roman_wfi_photom_0046.asdf
+    "photom": {"pixelareasr": {"F062": 2.640600009241359e-13 * u.sr,
+                               "F087": 2.640600009241359e-13 * u.sr,
+                               "F106": 2.640600009241359e-13 * u.sr,
+                               "F129": 2.640600009241359e-13 * u.sr,
+                               "F146": 2.640600009241359e-13 * u.sr,
+                               "F158": 2.640600009241359e-13 * u.sr,
+                               "F184": 2.640600009241359e-13 * u.sr,
+                               "F213": 2.640600009241359e-13 * u.sr,
+                               }
+               }
 }
 
 nborder = 4  # number of border pixels used for reference pixels.
@@ -140,6 +151,9 @@
     "pixel_depth": 5,
 }
 
+# Centermost PSF to use for mosaic creation
+default_sca = 2
+
 # Persistence defaults
 # delete persistence records fainter than 0.01 electron / s
 # e.g., MA table 1 has ~144 s, so this is ~1 electron over the whole exposure.
diff --git a/romanisim/psf.py b/romanisim/psf.py
index 40294500..a669751f 100644
--- a/romanisim/psf.py
+++ b/romanisim/psf.py
@@ -201,8 +201,8 @@ def at_position(self, x, y):
         # x = [0, off] -> 1
         # x = [npix, infinity] -> 0
         # linearly between those, likewise for y.
-        out = (self.psf['ll'] * wleft * wlow +
-               self.psf['lr'] * (1 - wleft) * wlow +
-               self.psf['ul'] * wleft * (1 - wlow) +
-               self.psf['ur'] * (1 - wleft) * (1 - wlow))
+        out = (self.psf['ll'] * wleft * wlow
+               + self.psf['lr'] * (1 - wleft) * wlow
+               + self.psf['ul'] * wleft * (1 - wlow)
+               + self.psf['ur'] * (1 - wleft) * (1 - wlow))
         return out
diff --git a/romanisim/tests/test_image.py b/romanisim/tests/test_image.py
index 8c2f4de9..4a76b833 100644
--- a/romanisim/tests/test_image.py
+++ b/romanisim/tests/test_image.py
@@ -18,7 +18,7 @@
 import numpy as np
 import galsim
 from galsim import roman
-from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1, l3
+from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1
 from astropy.coordinates import SkyCoord
 from astropy import units as u
 from astropy.time import Time
@@ -30,7 +30,6 @@
 from metrics_logger.decorators import metrics_logger
 from romanisim import log
 from roman_datamodels.stnode import WfiScienceRaw, WfiImage
-import roman_datamodels.maker_utils as maker_utils
 import romanisim.bandpass
 
 
@@ -732,98 +731,6 @@ def test_inject_source_into_image():
         af.write_to(os.path.join(artifactdir, 'dms231.asdf'))
 
 
-@metrics_logger("DMS232")
-@pytest.mark.soctests
-def test_inject_sources_into_mosaic():
-    """Inject sources into a mosaic.
-    """
-
-    # Set constants and metadata
-    galsim.roman.n_pix = 200
-    rng_seed = 42
-    metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary)
-    metadata['basic']['optical_element'] = 'F158'
-
-    # Create WCS
-    twcs = wcs.get_mosaic_wcs(metadata)
-
-    # Create initial Level 3 mosaic
-
-    # Create Four-quadrant pattern of gaussian noise, centered around one
-    # Each quadrant's gaussian noise scales like total exposure time
-    # (total files contributed to each quadrant)
-
-    # Create gaussian noise generators
-    g1 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.01)
-    g2 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.02)
-    g3 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.05)
-    g4 = galsim.GaussianDeviate(rng_seed, mean=1.0, sigma=0.1)
-
-    # Create level 3 mosaic model
-    l3_mos = maker_utils.mk_level3_mosaic(shape=(galsim.roman.n_pix, galsim.roman.n_pix))
-
-    # Update metadata in the l3 model
-    for key in metadata.keys():
-        if key in l3_mos.meta:
-            l3_mos.meta[key].update(metadata[key])
-
-    # Populate the mosaic data array with gaussian noise from generators
-    g1.generate(l3_mos.data.value[0:100, 0:100])
-    g2.generate(l3_mos.data.value[0:100, 100:200])
-    g3.generate(l3_mos.data.value[100:200, 0:100])
-    g4.generate(l3_mos.data.value[100:200, 100:200])
-
-    # Define Poisson Noise of mosaic
-    l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2
-    l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2
-    l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2
-    l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2
-
-    # Create normalized psf source catalog (same source in each quadrant)
-    sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0],
-               "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], "F158": 4 * [1.0]}
-    sc_table = table.Table(sc_dict)
-
-    xpos, ypos = 50, 50
-    sc_table["ra"][0], sc_table["dec"][0] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value
-    xpos, ypos = 50, 150
-    sc_table['ra'][1], sc_table['dec'][1] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value
-    xpos, ypos = 150, 50
-    sc_table['ra'][2], sc_table['dec'][2] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value
-    xpos, ypos = 150, 150
-    sc_table['ra'][3], sc_table['dec'][3] = (twcs._radec(xpos, ypos) * u.rad).to(u.deg).value
-
-    source_cat = catalog.table_to_catalog(sc_table, ["F158"])
-
-    # Copy original Mosaic before adding sources
-    l3_mos_orig = l3_mos.copy()
-
-    # Add source_cat objects to mosaic
-    l3.add_objects_to_l3(l3_mos, source_cat, seed=rng_seed)
-
-    # Ensure that every data pixel value has increased or
-    # remained the same with the new sources injected
-    assert np.all(l3_mos.data.value >= l3_mos_orig.data.value)
-
-    # Ensure that every pixel's poisson variance has increased or
-    # remained the same with the new sources injected
-    # Numpy isclose is needed to determine equality, due to float precision issues
-    close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06)
-    assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask])
-
-    # Create log entry and artifacts
-    log.info('DMS232 successfully injected sources into a mosiac at points (50,50), (50,150), (150,50), (150,150).')
-
-    artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
-    if artifactdir is not None:
-        af = asdf.AsdfFile()
-        af.tree = {'l3_mos': l3_mos,
-                   'l3_mos_orig': l3_mos_orig,
-                   'source_cat_table': sc_table,
-                   }
-        af.write_to(os.path.join(artifactdir, 'dms232.asdf'))
-
-
 @metrics_logger("DMS228")
 @pytest.mark.soctests
 def test_image_input(tmpdir):
diff --git a/romanisim/tests/test_l1.py b/romanisim/tests/test_l1.py
index d4a7aea4..5f944872 100644
--- a/romanisim/tests/test_l1.py
+++ b/romanisim/tests/test_l1.py
@@ -29,7 +29,7 @@
 
 read_pattern_list = [
     [[1 + x for x in range(10)], [11], [12 + x for x in range(10)], [30],
-         [40 + x for x in range(5)], [50 + x for x in range(100)]],
+     [40 + x for x in range(5)], [50 + x for x in range(100)]],
     [[1]],
     [[1], [10]],
 ]
diff --git a/romanisim/tests/test_l3.py b/romanisim/tests/test_l3.py
new file mode 100644
index 00000000..123e4d63
--- /dev/null
+++ b/romanisim/tests/test_l3.py
@@ -0,0 +1,150 @@
+"""Unit tests for mosaic module.
+
+"""
+
+import os
+import copy
+import math
+import numpy as np
+import galsim
+from galsim import roman
+from romanisim import image, parameters, catalog, psf, util, wcs, persistence, ramp, l1, l3
+from astropy.coordinates import SkyCoord
+from astropy import units as u
+from astropy.time import Time
+from astropy import table
+import asdf
+import webbpsf
+from astropy.modeling.functional_models import Sersic2D
+import pytest
+from metrics_logger.decorators import metrics_logger
+from romanisim import log
+from roman_datamodels.stnode import WfiScienceRaw, WfiImage
+import roman_datamodels.maker_utils as maker_utils
+import romanisim.bandpass
+
+
+@metrics_logger("DMS232")
+@pytest.mark.soctests
+def test_inject_sources_into_mosaic():
+    """Inject sources into a mosaic.
+    """
+
+    # Set constants and metadata
+    galsim.roman.n_pix = 200
+    rng_seed = 42
+    metadata = copy.deepcopy(parameters.default_mosaic_parameters_dictionary)
+    filter_name = 'F158'
+    metadata['basic']['optical_element'] = filter_name
+
+    # Create WCS
+    twcs = wcs.get_mosaic_wcs(metadata, shape=(galsim.roman.n_pix, galsim.roman.n_pix))
+
+    # Create initial Level 3 mosaic
+
+    # Create Four-quadrant pattern of gaussian noise, centered around one
+    # Each quadrant's gaussian noise scales like total exposure time
+    # (total files contributed to each quadrant)
+
+    # Create gaussian noise generators
+    g1 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.01e-7)
+    g2 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.02e-7)
+    g3 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.05e-7)
+    g4 = galsim.GaussianDeviate(rng_seed, mean=1.0e-7, sigma=0.1e-7)
+
+    # Create level 3 mosaic model
+    l3_mos = maker_utils.mk_level3_mosaic(shape=(galsim.roman.n_pix, galsim.roman.n_pix))
+    l3_mos['meta']['wcs'] = twcs
+
+    # Update metadata in the l3 model
+    for key in metadata.keys():
+        if key in l3_mos.meta:
+            l3_mos.meta[key].update(metadata[key])
+
+    # Obtain unit conversion factors
+    # Need to convert from counts / pixel to MJy / sr
+    # Flux to counts
+    cps_conv = romanisim.bandpass.get_abflux(filter_name)
+    # Unit factor
+    unit_factor = ((3631 * u.Jy) / (romanisim.bandpass.get_abflux(filter_name) * 10e6
+                                    * parameters.reference_data['photom']["pixelareasr"][filter_name])).to(u.MJy / u.sr)
+
+    # Populate the mosaic data array with gaussian noise from generators
+    g1.generate(l3_mos.data.value[0:100, 0:100])
+    g2.generate(l3_mos.data.value[0:100, 100:200])
+    g3.generate(l3_mos.data.value[100:200, 0:100])
+    g4.generate(l3_mos.data.value[100:200, 100:200])
+
+    # Define Poisson Noise of mosaic
+    l3_mos.var_poisson.value[0:100, 0:100] = 0.01**2
+    l3_mos.var_poisson.value[0:100, 100:200] = 0.02**2
+    l3_mos.var_poisson.value[100:200, 0:100] = 0.05**2
+    l3_mos.var_poisson.value[100:200, 100:200] = 0.1**2
+
+    # Create normalized psf source catalog (same source in each quadrant)
+    mag_flux = 1e-10
+    sc_dict = {"ra": 4 * [0.0], "dec": 4 * [0.0], "type": 4 * ["PSF"], "n": 4 * [-1.0],
+               "half_light_radius": 4 * [0.0], "pa": 4 * [0.0], "ba": 4 * [1.0], filter_name: 4 * [mag_flux]}
+    sc_table = table.Table(sc_dict)
+
+    # Set locations
+    xpos_idx = [50, 50, 150, 150]
+    ypos_idx = [50, 150, 50, 150]
+
+    # Populate flux scaling ratio and catalog
+    Ct = []
+    for idx, (x, y) in enumerate(zip(xpos_idx, ypos_idx)):
+        # Set scaling factor for injected sources
+        # Flux / sigma_p^2
+        if l3_mos.var_poisson[y, x].value != 0:
+            Ct.append(math.fabs(l3_mos.data[y, x].value / l3_mos.var_poisson[y, x].value))
+        else:
+            Ct.append(1.0)
+
+        sc_table["ra"][idx], sc_table["dec"][idx] = (twcs._radec(x, y) * u.rad).to(u.deg).value
+
+    source_cat = catalog.table_to_catalog(sc_table, [filter_name])
+
+    # Copy original Mosaic before adding sources as sources are added in place
+    l3_mos_orig = l3_mos.copy()
+    l3_mos_orig.data = l3_mos.data.copy()
+    l3_mos_orig.var_poisson = l3_mos.var_poisson.copy()
+
+    # Add source_cat objects to mosaic
+    l3.add_objects_to_l3(l3_mos, source_cat, Ct, cps_conv=cps_conv, unit_factor=unit_factor.value, seed=rng_seed)
+
+    # Create overall scaling factor map
+    Ct_all = np.divide(l3_mos_orig.data.value, l3_mos_orig.var_poisson.value,
+                       out=np.ones(l3_mos_orig.data.shape),
+                       where=l3_mos_orig.var_poisson.value != 0)
+
+    # Set new poisson variance
+    l3_mos.var_poisson = (l3_mos.data.value / Ct_all) * l3_mos.var_poisson.unit
+
+    # Ensure that every data pixel value has increased or
+    # remained the same with the new sources injected
+    assert np.all(l3_mos.data.value >= l3_mos_orig.data.value)
+
+    # Ensure that every pixel's poisson variance has increased or
+    # remained the same with the new sources injected
+    # Numpy isclose is needed to determine equality, due to float precision issues
+    close_mask = np.isclose(l3_mos.var_poisson.value, l3_mos_orig.var_poisson.value, rtol=1e-06)
+    assert False in close_mask
+    assert np.all(l3_mos.var_poisson.value[~close_mask] > l3_mos_orig.var_poisson.value[~close_mask])
+
+    # Ensure total added flux matches expected added flux
+    total_rec_flux = np.sum(l3_mos.data - l3_mos_orig.data) / unit_factor
+    total_theo_flux = 4 * mag_flux * cps_conv
+    assert np.isclose(total_rec_flux, total_theo_flux, rtol=4e-02)
+
+    # Create log entry and artifacts
+    log.info('DMS232 successfully injected sources into a mosaic at points (50,50), (50,150), (150,50), (150,150).')
+
+    artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
+    if artifactdir is not None:
+        af = asdf.AsdfFile()
+        af.tree = {'l3_mos': l3_mos,
+                   'l3_mos_orig': l3_mos_orig,
+                   'source_cat_table': sc_table,
+                   }
+        af.write_to(os.path.join(artifactdir, 'dms232.asdf'))
diff --git a/romanisim/util.py b/romanisim/util.py
index 4776ec71..a4da51ca 100644
--- a/romanisim/util.py
+++ b/romanisim/util.py
@@ -335,6 +335,90 @@ def sample_king_distances(rc, rt, npts, rng=None):
     return radii
 
 
+def decode_context_times(context, exptimes):
+    """
+    Get 0-based indices of input images that contributed to (resampled)
+    output pixel with coordinates ``x`` and ``y``.
+
+    Parameters
+    ----------
+    context: numpy.ndarray
+        A 3D `~numpy.ndarray` of integral data type.
+
+    exptimes: list of floats
+        Exposure times for each component image.
+
+
+    Returns
+    -------
+    total_exptimes: numpy.ndarray
+        A 2D array of total exposure time for each pixel.
+    """
+
+    if context.ndim != 3:
+        raise ValueError("'context' must be a 3D array.")
+
+    """
+    Context decoding example:
+    An example context array for an output image of array shape ``(5, 6)``
+    obtained by resampling 80 input images.
+
+    .. code-block:: none
+
+        con = np.array(
+            [[[0, 0, 0, 0, 0, 0],
+              [0, 0, 0, 36196864, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 537920000, 0, 0, 0]],
+             [[0, 0, 0, 0, 0, 0,],
+              [0, 0, 0, 67125536, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 163856, 0, 0, 0]],
+             [[0, 0, 0, 0, 0, 0],
+              [0, 0, 0, 8203, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 0, 0, 0, 0],
+              [0, 0, 32865, 0, 0, 0]]],
+            dtype=np.int32
+        )
+        decode_context(con, [3, 2], [1, 4])
+        [array([ 9, 12, 14, 19, 21, 25, 37, 40, 46, 58, 64, 65, 67, 77]),
+        array([ 9, 20, 29, 36, 47, 49, 64, 69, 70, 79])]
+    """
+
+    nbits = 8 * context.dtype.itemsize
+
+    total_exptimes = np.zeros(context.shape[1:])
+
+    for x in range(total_exptimes.shape[0]):
+        for y in range(total_exptimes.shape[1]):
+            files = [v & (1 << k) for v in context[:, x, y] for k in range(nbits)]
+            tot_time = 0
+            files = [file for file in files if (file != 0)]
+
+            for im_idx in files:
+                tot_time += exptimes[im_idx - 1]
+
+            total_exptimes[x][y] = tot_time
+
+    def sum_times(x):
+        tot_time = 0
+        files = [x & (1 << k) for k in range(nbits)]
+        files = [file for file in files if (file != 0)]
+        for im_idx in files:
+            tot_time += exptimes[im_idx - 1]
+        return tot_time
+
+    vectorized_sum_times = np.vectorize(sum_times)
+
+    total_exptimes = vectorized_sum_times(context[:,])
+    total_exptimes = total_exptimes.reshape(total_exptimes.shape[1:])
+
+    return total_exptimes
+
+
 def default_image_meta(time=None, ma_table=1, filter_name='F087',
                        detector='WFI01', coord=None):
     """Return some simple default metadata for input to image.simulate
diff --git a/romanisim/wcs.py b/romanisim/wcs.py
index d1f2c6de..c0f7c289 100644
--- a/romanisim/wcs.py
+++ b/romanisim/wcs.py
@@ -17,6 +17,7 @@
 """
 
 import warnings
+import math
 import numpy as np
 import astropy.coordinates
 from astropy import units as u
@@ -219,15 +220,15 @@ def make_wcs(targ_pos,
     # V2V3 are in arcseconds, while SphericalToCartesian expects degrees,
     # so again start by scaling by 3600
     tel2sky = ((models.Scale(1 / 3600) & models.Scale(1 / 3600))
-                | gwcs.geometry.SphericalToCartesian(wrap_lon_at=wrap_v2_at)
-                | rot
-                | gwcs.geometry.CartesianToSpherical(wrap_lon_at=wrap_lon_at))
+               | gwcs.geometry.SphericalToCartesian(wrap_lon_at=wrap_v2_at)
+               | rot
+               | gwcs.geometry.CartesianToSpherical(wrap_lon_at=wrap_lon_at))
     tel2sky.name = 'v23tosky'
 
     detector = cf.Frame2D(name='detector', axes_order=(0, 1),
                           unit=(u.pix, u.pix))
     v2v3 = cf.Frame2D(name="v2v3", axes_order=(0, 1),
-        axes_names=("v2", "v3"), unit=(u.arcsec, u.arcsec))
+                      axes_names=("v2", "v3"), unit=(u.arcsec, u.arcsec))
     world = cf.CelestialFrame(reference_frame=astropy.coordinates.ICRS(),
                               name='world')
 
@@ -416,7 +417,7 @@ def convert_wcs_to_gwcs(wcs):
         return wcs_from_fits_header(wcs.header.header)
 
 
-def get_mosaic_wcs(mosaic):
+def get_mosaic_wcs(mosaic, shape=None):
     """Get a WCS object for a given sca or set of CRDS parameters.
 
     Parameters
@@ -431,13 +432,16 @@ def get_mosaic_wcs(mosaic):
 
     # If sent a dictionary, create a temporary model for data interface
     if (type(mosaic) is not roman_datamodels.datamodels.MosaicModel):
-        mosaic_node = maker_utils.mk_level3_mosaic()
+        if shape is None:
+            mosaic_node = maker_utils.mk_level3_mosaic()
+        else:
+            mosaic_node = maker_utils.mk_level3_mosaic(shape=shape)
         for key in mosaic.keys():
             if isinstance(mosaic[key], dict):
                 mosaic_node['meta'][key].update(mosaic[key])
             else:
                 mosaic_node['meta'][key] = mosaic[key]
-        mosaic_node = roman_datamodels.datamodels.MosaicModel(mosaic_node)
+        # mosaic_node = roman_datamodels.datamodels.MosaicModel(mosaic_node)
     else:
         mosaic_node = mosaic
 
@@ -445,11 +449,14 @@ def get_mosaic_wcs(mosaic):
         mosaic_node.meta.wcsinfo.ra_ref * u.deg,
         mosaic_node.meta.wcsinfo.dec_ref * u.deg)
 
+    if shape is None:
+        shape = (mosaic_node.data.shape[0],
+                 mosaic_node.data.shape[1])
+
     # Create a tangent plane WCS for the mosaic
     # The affine parameters below should be reviewed and updated
     affine = galsim.AffineTransform(
-        0.1, 0, 0, 0.1, origin=galsim.PositionI(mosaic_node.data.shape[0] / 2.0,
-                                                mosaic_node.data.shape[1] / 2.0,),
+        0.1, 0, 0, 0.1, origin=galsim.PositionI(math.ceil(shape[0] / 2.0), math.ceil(shape[1] / 2.0)),
         world_origin=galsim.PositionD(0, 0))
     wcs = galsim.TanWCS(affine,
                         util.celestialcoord(world_pos))
diff --git a/setup.py b/setup.py
index d2bb86d4..b9c44806 100644
--- a/setup.py
+++ b/setup.py
@@ -8,7 +8,7 @@
 Options.docstrings = True
 Options.annotate = False
 
-# importing these extension modules is tested in `.github/workflows/build.yml`; 
+# importing these extension modules is tested in `.github/workflows/build.yml`;
 # when adding new modules here, make sure to add them to the `test_command` entry there
 extensions = [Extension('romanisim.ramp_fit_casertano',
                         ['romanisim/ramp_fit_casertano.pyx'],