diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..410edf52
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+/build/*
+/dist/*
+/roman_imsim.egg-info/*
\ No newline at end of file
diff --git a/config/was.yaml b/config/was.yaml
index 8c9a1640..589072d9 100644
--- a/config/was.yaml
+++ b/config/was.yaml
@@ -88,16 +88,42 @@ image:
     # read_noise: False
     # sky_subtract: False
 
-    ignore_noise: False
-    stray_light: True
-    thermal_background: True
-    reciprocity_failure: True
-    dark_current: True
-    nonlinearity: True
-    ipc: True
-    read_noise: True
+    add_effects:
+        Background:
+            model: simple_model
+            thermal_background: True
+            stray_light: True
+        QuantumEfficiency: 
+            model: lab_model
+        BrighterFatter: 
+            model: lab_model
+        Persistence:
+            model: lab_model
+        # ReciprocityFailure:
+        #     model: simple_model
+        DarkCurrent: 
+            model: lab_model
+        Saturation:
+            model: lab_model
+        Nonlinearity: 
+            model: lab_model
+        IPC:
+            model: lab_model
+        DeadPixel:
+            model: lab_model
+        VTPE:
+            model: lab_model
+        ReadNoise: 
+            model: lab_model
+        Gain:
+            model: lab_model
+        Bias:
+            model: lab_model
     sky_subtract: False
 
+    sca_filepath: /hpc/group/cosmology/phy-lsst/roman_sensors
+    save_diff: False
+
     # nobjects: 500
 
 stamp:
@@ -124,17 +150,14 @@ gal:
 
 input:
     obseq_data:
-        # file_name: /lus/grand/projects/RomanDESC/final_roman_runfiles/was/Roman_WAS_obseq_11_1_23.fits
         file_name: /hpc/group/cosmology/OpenUniverse2024/RomanWAS/Roman_WAS_obseq_11_1_23.fits
-        visit: 1
+        visit: 4906
         SCA: '@image.SCA'
     roman_psf:
         SCA: '@image.SCA'
         n_waves: 5
     sky_catalog:
-        # file_name: /lus/grand/projects/RomanDESC/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml 
-        # file_name: /hpc/group/cosmology/OpenUniverse2024/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml 
-        file_name: /hpc/home/yf194/Work/projects/roman_imsim/config/skyCatalog.yaml
+        file_name: /hpc/group/cosmology/OpenUniverse2024/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml
         edge_pix: 512
         mjd: { type: ObSeqData, field: mjd }
         exptime: { type: ObSeqData, field: exptime }
diff --git a/roman_imsim/__init__.py b/roman_imsim/__init__.py
index 72da6e7a..2cbd35d0 100644
--- a/roman_imsim/__init__.py
+++ b/roman_imsim/__init__.py
@@ -16,3 +16,6 @@
 from .skycat import *
 from .stamp import *
 from .wcs import *
+from .skycat import *
+from .photonOps import *
+from .bandpass import *
diff --git a/roman_imsim/detector_physics.py b/roman_imsim/detector_physics.py
index 7e58f8d4..a2028ae8 100644
--- a/roman_imsim/detector_physics.py
+++ b/roman_imsim/detector_physics.py
@@ -83,15 +83,6 @@ def __init__(self, params, visit, SCA):
         )[self.sca]
         self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)))
 
-        self.WCS = roman.getWCS(
-            world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]),
-            PA=obseq_data.ob["pa"],
-            date=self.date,
-            SCAs=self.sca,
-            PA_is_FPA=True,
-        )[self.sca]
-        self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)))
-
 
 class modify_image(object):
     """
@@ -121,8 +112,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
         if sca_filepath is not None:
             self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca])
             print("------- Using SCA files --------")
-            self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca])
-            print("------- Using SCA files --------")
         else:
             self.df = None
             print("------- Using simple detector model --------")
@@ -139,16 +128,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
             "truth", self.get_path_name(use_galsim=use_galsim)
         )
 
-        b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix)
-        old_filename = os.path.join(self.params["output"]["dir"], imfilename)
-        if not os.path.exists(
-            self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))
-        ):
-            os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)))
-        new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace(
-            "truth", self.get_path_name(use_galsim=use_galsim)
-        )
-
         b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix)
         im = fio.FITS(old_filename)[-1].read()
         im = galsim.Image(im, bounds=b, wcs=self.pointing.WCS)
@@ -160,7 +139,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
                 force_cvz = True
         self.setup_sky(im, self.pointing, rng, visit * sca, force_cvz=force_cvz)
 
-        img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim)
         img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim)
 
         write_fits(
@@ -172,7 +150,6 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_
             self.pointing.sca,
             sky_mean=sky_mean,
         )
-        write_fits(old_filename, new_filename, img, sky_noise, dq, self.pointing.sca, sky_mean=sky_mean)
 
     def get_path_name(self, use_galsim=False):
 
@@ -562,10 +539,6 @@ def bfe(self, im):
         # The img is clipped by the saturation level here to cap the brighter fatter effect
         # and avoid unphysical behavior
 
-        array_pad = self.saturate(im.copy()).array[4:-4, 4:-4]  # img of interest 4088x4088
-        array_pad = np.pad(
-            array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric"
-        )  # 4100x4100 array
         array_pad = self.saturate(im.copy()).array[4:-4, 4:-4]  # img of interest 4088x4088
         array_pad = np.pad(
             array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric"
@@ -697,8 +670,6 @@ def get_eff_sky_bg(self, pointing, radec):
         radec               : World coordinate position of image
         """
 
-        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
-        sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2
         sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
         sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2
 
@@ -739,9 +710,6 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False):
             self.dark_current_ = (
                 roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime
             )
-            self.dark_current_ = (
-                roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime
-            )
         if self.df is None:
             self.gain = roman.gain
         else:
@@ -755,8 +723,6 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False):
             radec = pointing.radec
         sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
         sky_level *= 1.0 + roman.stray_light_fraction
-        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
-        sky_level *= 1.0 + roman.stray_light_fraction
         # Make a image of the sky that takes into account the spatially variable pixel scale. Note
         # that makeSkyImage() takes a bit of time. If you do not care about the variable pixel
         # scale, you could simply compute an approximate sky level in e-/pix by multiplying
@@ -864,8 +830,6 @@ def dark_current(self, im):
             # opt for numpy random geneator instead for speed
             self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype)
             im.array[:, :] += self.im_dark
-            self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype)
-            im.array[:, :] += self.im_dark
 
         # NOTE: Sky level and dark current might appear like a constant background that can be
         # simply subtracted. However, these contribute to the shot noise and matter for the
@@ -960,11 +924,6 @@ def add_persistence(self, im, pointing):
         p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected])
         dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list])
         p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))]
-        dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0]
-        dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10]
-        p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected])
-        dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list])
-        p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))]
 
         if self.df is None:
             # iterate over previous exposures
@@ -986,7 +945,6 @@ def add_persistence(self, im, pointing):
                 x = x.clip(0)  # remove negative stimulus
 
                 im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime
-                im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime
 
         else:
 
@@ -1106,8 +1064,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
             num_grids = 4  # num_grids <= 8
             grid_size = 4096 // num_grids
 
-            array_pad = im.array[4:-4, 4:-4]  # it's an array instead of img
-            array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric")  # 4098x4098 array
             array_pad = im.array[4:-4, 4:-4]  # it's an array instead of img
             array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric")  # 4098x4098 array
 
@@ -1127,8 +1083,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
                         for i in range(3):
                             tmp = (t.dot(K[j, i, :, :])).dot(t.T)  # grid_sizexgrid_size
                             K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric")
-                            tmp = (t.dot(K[j, i, :, :])).dot(t.T)  # grid_sizexgrid_size
-                            K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric")
 
                     for dy in range(-1, 2):
                         for dx in range(-1, 2):
@@ -1177,11 +1131,6 @@ def add_read_noise(self, im):
                 self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype)
             )
             im.array[:, :] += self.im_read
-            read_noise = self.df["READ"][2, :, :].flatten()  # flattened 4096x4096 array
-            self.im_read = (
-                self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype)
-            )
-            im.array[:, :] += self.im_read
 
             # noise_array = self.rng_np.normal(loc=0., scale=read_noise)
             # 4088x4088 img
diff --git a/roman_imsim/effects/__init__.py b/roman_imsim/effects/__init__.py
new file mode 100644
index 00000000..7a1faf33
--- /dev/null
+++ b/roman_imsim/effects/__init__.py
@@ -0,0 +1,15 @@
+from .roman_effects import RomanEffects, setup_sky
+from .brighter_fatter import BrighterFatter
+from .nonlinearity import Nonlinearity
+from .background import Background
+from .quantum_efficiency import QuantumEfficiency
+from .persistence import Persistence
+from .reciprocity_failure import ReciprocityFailure
+from .dark_current import DarkCurrent
+from .saturation import Saturation
+from .ipc import IPC
+from .dead_pixel import DeadPixel
+from .vtpe import VTPE
+from .read_noise import ReadNoise
+from .gain import Gain
+from .bias import Bias
diff --git a/roman_imsim/effects/background.py b/roman_imsim/effects/background.py
new file mode 100644
index 00000000..b16b7374
--- /dev/null
+++ b/roman_imsim/effects/background.py
@@ -0,0 +1,54 @@
+import os
+import galsim
+from . import RomanEffects
+import galsim.roman as roman
+
+
+class Background(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+        self.thermal_background = (
+            self.params["thermal_background"] if "thermal_background" in self.params else False
+        )
+        self.stray_light = self.params["stray_light"] if "stray_light" in self.params else False
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        if self.save_diff:
+            orig = image.copy()
+
+        self.logger.warning("Simple model will be applied for background.")
+        pointing = self.pointing
+        # Build current specification sky level if sky level not given
+        if self.force_cvz:
+            radec = self.translate_cvz(pointing.radec)
+        else:
+            radec = pointing.radec
+        sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date)
+        self.logger.debug("Adding sky_level = %s", sky_level)
+
+        if self.stray_light:
+            self.logger.debug("Stray light fraction = %s", roman.stray_light_fraction)
+            sky_level *= 1.0 + roman.stray_light_fraction
+        # Create sky image
+        self.sky = galsim.Image(bounds=image.bounds, wcs=pointing.WCS)
+        pointing.WCS.makeSkyImage(self.sky, sky_level)
+        if self.thermal_background:
+            tb = roman.thermal_backgrounds[pointing.filter] * pointing.exptime
+            self.logger.debug("Adding thermal background: %s", tb)
+            self.sky += tb
+        self.sky.addNoise(self.noise)
+
+        # [TODO] Not entirely sure about this block, since the 'auto' option is meant to
+        # let the software choose which drawing method to use based on the total flux.
+        if self.base["image"]["draw_method"] not in ["phot", "auto"]:
+            image.addNoise(self.noise)
+
+        # Adding sky level to the image.
+        image += self.sky[self.sky.bounds & image.bounds]
+        if self.save_diff:
+            prev = image.copy()
+            diff = prev - orig
+            diff.write(os.path.join(self.diff_dir, "sky_a.fits"))
+        return image
diff --git a/roman_imsim/effects/bias.py b/roman_imsim/effects/bias.py
new file mode 100644
index 00000000..54a3a574
--- /dev/null
+++ b/roman_imsim/effects/bias.py
@@ -0,0 +1,27 @@
+import os
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class Bias(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("No bias will be applied.")
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No bias data provided; no bias will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for bias.")
+        bias = self.df["BIAS"][:, :]  # 4096x4096 img
+
+        image.array[:, :] += bias
+        return image
diff --git a/roman_imsim/effects/brighter_fatter.py b/roman_imsim/effects/brighter_fatter.py
new file mode 100644
index 00000000..bad18085
--- /dev/null
+++ b/roman_imsim/effects/brighter_fatter.py
@@ -0,0 +1,222 @@
+import os
+import numpy as np
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class BrighterFatter(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("No bfe effect will be applied.")
+        return image
+
+    def lab_model(self, image):
+        """
+        Apply brighter-fatter effect.
+        Brighter fatter effect is a non-linear effect that deflects photons due to the
+        the eletric field built by the accumulated charges. This effect exists in both
+        CCD and CMOS detectors and typically percent level change in charge.
+        The built-in electric field by the charges in pixels tends to repulse charges
+        to nearby pixels. Thus, the profile of more illuminous ojbect becomes broader.
+        This effect can also be understood effectly as change in pixel area and pixel
+        boundaries.
+        BFE is defined in terms of the Antilogus coefficient kernel of total pixel area change
+        in the detector effect charaterization file. Kernel of the total pixel area, however,
+        is not sufficient. Image simulation of the brighter fatter effect requires the shift
+        of the four pixel boundaries. Before we get better data, we solve for the boundary
+        shift components from the kernel of total pixel area by assumming several symmetric constraints.
+        Input
+        im                                      : Image
+        BFE[nbfe+Delta y, nbfe+Delta x, y, x]   : bfe coefficient kernel, nbfe=2
+        """
+        if self.sca_filepath is None:
+            self.logger.warning("No BFE kernel data file provided; no bfe effect will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for brighter-fatter effect.")
+
+        nbfe = 2  # kernel of bfe in shape (2 x nbfe+1)*(2 x nbfe+1)
+        bin_size = 128
+        n_max = 32
+        m_max = 32
+        num_grids = 4
+        n_sub = n_max // num_grids
+        m_sub = m_max // num_grids
+
+        # =======================================================================
+        # solve boundary shfit kernel aX components
+        # =======================================================================
+        a_area = self.df["BFE"][:, :, :, :]  # 5x5x32x32
+        a_components = np.zeros((4, 2 * nbfe + 1, 2 * nbfe + 1, n_max, m_max))  # 4x5x5x32x32
+
+        # solve aR aT aL aB for each a
+        for n in range(n_max):  # m_max and n_max = 32 (binned in 128x128)
+            for m in range(m_max):
+                a = a_area[:, :, n, m]  # a in (2 x nbfe+1)*(2 x nbfe+1)
+
+                # assume two parity symmetries
+                a = (a + np.fliplr(a) + np.flipud(a) + np.flip(a)) / 4.0
+
+                r = 0.5 * (3.25 / 4.25) ** (1.5) / 1.5  # source-boundary projection
+                B = (a[2, 2], a[3, 2], a[2, 3], a[3, 3], a[4, 2], a[2, 4], a[3, 4], a[4, 4])
+
+                A = np.array(
+                    [
+                        [-2, -2, 0, 0, 0, 0, 0],
+                        [0, 1, 0, -1, -2, 0, 0],
+                        [1, 0, -1, 0, -2, 0, 0],
+                        [0, 0, 0, 0, 2, -2, 0],
+                        [0, 0, 0, 1, 0, -2 * r, 0],
+                        [0, 0, 1, 0, 0, -2 * r, 0],
+                        [0, 0, 0, 0, 0, 1 + r, -1],
+                        [0, 0, 0, 0, 0, 0, 2],
+                    ]
+                )
+
+                s1, s2, s3, s4, s5, s6, s7 = np.linalg.lstsq(A, B, rcond=None)[0]
+
+                aR = np.array(
+                    [
+                        [0.0, -s7, -r * s6, r * s6, s7],
+                        [0.0, -s6, -s5, s5, s6],
+                        [0.0, -s3, -s1, s1, s3],
+                        [0.0, -s6, -s5, s5, s6],
+                        [0.0, -s7, -r * s6, r * s6, s7],
+                    ]
+                )
+
+                aT = np.array(
+                    [
+                        [0.0, 0.0, 0.0, 0.0, 0.0],
+                        [-s7, -s6, -s4, -s6, -s7],
+                        [-r * s6, -s5, -s2, -s5, -r * s6],
+                        [r * s6, s5, s2, s5, r * s6],
+                        [s7, s6, s4, s6, s7],
+                    ]
+                )
+
+                aL = aR[::-1, ::-1]
+                aB = aT[::-1, ::-1]
+
+                a_components[0, :, :, n, m] = aR[:, :]
+                a_components[1, :, :, n, m] = aT[:, :]
+                a_components[2, :, :, n, m] = aL[:, :]
+                a_components[3, :, :, n, m] = aB[:, :]
+
+        # =============================
+        # Apply bfe to image
+        # =============================
+
+        # pad and expand kernels
+        # The img is clipped by the saturation level here to cap the brighter fatter effect
+        # and avoid unphysical behavior
+
+        # array_pad = image.copy().array
+        # saturation_array = np.ones_like(array_pad) * self.saturation_level
+        # where_sat = np.where(array_pad > saturation_array)
+        # array_pad[ where_sat ] = saturation_array[ where_sat ]
+        # array_pad = array_pad[4:-4,4:-4]
+        saturate = self.cross_refer("Saturation")
+        array_pad = saturate.apply(image=image.copy()).array[4:-4, 4:-4]  # img of interest 4088x4088
+        array_pad = np.pad(
+            array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric"
+        )  # 4100x4100 array
+
+        # (4, 4096, 4096) in order of [aR, aT, aL, aB]
+        dQ_components = np.zeros((4, bin_size * n_max, bin_size * m_max))
+
+        # run in sub grids to reduce memory
+
+        # pad and expand kernels
+        t = np.zeros((bin_size * n_sub, n_sub))
+        for row in range(t.shape[0]):
+            t[row, row // (bin_size)] = 1
+
+        for gj in range(num_grids):
+            for gi in range(num_grids):
+
+                # (4,5,5,sub_grid,sub_grid)
+                a_components_pad = np.zeros(
+                    (4, 2 * nbfe + 1, 2 * nbfe + 1, bin_size * n_sub + 2 * nbfe, bin_size * m_sub + 2 * nbfe)
+                )
+
+                for comp in range(4):
+                    for j in range(2 * nbfe + 1):
+                        for i in range(2 * nbfe + 1):
+                            # sub_grid*sub_grid
+                            tmp = (
+                                t.dot(
+                                    a_components[
+                                        comp,
+                                        j,
+                                        i,
+                                        gj * n_sub : (gj + 1) * n_sub,
+                                        gi * m_sub : (gi + 1) * m_sub,
+                                    ]
+                                )
+                            ).dot(t.T)
+                            a_components_pad[comp, j, i, :, :] = np.pad(
+                                tmp, [(nbfe, nbfe), (nbfe, nbfe)], mode="symmetric"
+                            )
+
+                # convolve aX_ij with Q_ij
+                for comp in range(4):
+                    for dy in range(-nbfe, nbfe + 1):
+                        for dx in range(-nbfe, nbfe + 1):
+                            dQ_components[
+                                comp,
+                                gj * bin_size * n_sub : (gj + 1) * bin_size * n_sub,
+                                gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub,
+                            ] += (
+                                a_components_pad[
+                                    comp,
+                                    nbfe + dy,
+                                    nbfe + dx,
+                                    nbfe - dy : nbfe - dy + bin_size * n_sub,
+                                    nbfe - dx : nbfe - dx + bin_size * m_sub,
+                                ]
+                                * array_pad[
+                                    -dy
+                                    + nbfe
+                                    + gj * bin_size * n_sub : -dy
+                                    + nbfe
+                                    + (gj + 1) * bin_size * n_sub,
+                                    -dx
+                                    + nbfe
+                                    + gi * bin_size * m_sub : -dx
+                                    + nbfe
+                                    + (gi + 1) * bin_size * m_sub,
+                                ]
+                            )
+
+                    dj = int(np.sin(comp * np.pi / 2))
+                    di = int(np.cos(comp * np.pi / 2))
+
+                    dQ_components[
+                        comp,
+                        gj * bin_size * n_sub : (gj + 1) * bin_size * n_sub,
+                        gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub,
+                    ] *= 0.5 * (
+                        array_pad[
+                            nbfe + gj * bin_size * n_sub : nbfe + (gj + 1) * bin_size * n_sub,
+                            nbfe + gi * bin_size * m_sub : nbfe + (gi + 1) * bin_size * m_sub,
+                        ]
+                        + array_pad[
+                            dj + nbfe + gj * bin_size * n_sub : dj + nbfe + (gj + 1) * bin_size * n_sub,
+                            di + nbfe + gi * bin_size * m_sub : di + nbfe + (gi + 1) * bin_size * m_sub,
+                        ]
+                    )
+
+        image.array[:, :] -= dQ_components.sum(axis=0)
+        image.array[:, 1:] += dQ_components[0][:, :-1]
+        image.array[1:, :] += dQ_components[1][:-1, :]
+        image.array[:, :-1] += dQ_components[2][:, 1:]
+        image.array[:-1, :] += dQ_components[3][1:, :]
+
+        return image
diff --git a/roman_imsim/effects/dark_current.py b/roman_imsim/effects/dark_current.py
new file mode 100644
index 00000000..252ffb43
--- /dev/null
+++ b/roman_imsim/effects/dark_current.py
@@ -0,0 +1,39 @@
+import os
+import fitsio as fio
+import galsim
+import galsim.roman as roman
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class DarkCurrent(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for dark current.")
+        exptime = self.pointing.exptime
+        self.dark_current_ = roman.dark_current * exptime
+        self.im_dark = image.copy()
+        dark_current_ = self.dark_current_
+        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(self.rng, dark_current_))
+        image.addNoise(dark_noise)
+        self.im_dark = image - self.im_dark
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No dark current file provided; no dark current will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        exptime = self.pointing.exptime
+        self.logger.warning("Lab measured model will be applied for dark current.")
+        self.dark_current_ = roman.dark_current * exptime + self.df["DARK"][:, :].flatten() * exptime
+        dark_current_ = self.dark_current_.clip(0)
+        # opt for numpy random geneator instead for speed
+        self.im_dark = self.rng_np.poisson(dark_current_).reshape(image.array.shape).astype(image.dtype)
+        image.array[:, :] += self.im_dark
+        return image
diff --git a/roman_imsim/effects/dead_pixel.py b/roman_imsim/effects/dead_pixel.py
new file mode 100644
index 00000000..320bddf7
--- /dev/null
+++ b/roman_imsim/effects/dead_pixel.py
@@ -0,0 +1,27 @@
+import os
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class DeadPixel(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("No dead pixel will be applied.")
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No bad pixel data file provided; no dead pixel will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for dead pixel.")
+        dead_mask = self.df["BADPIX"][:, :] & 1  # 4096x4096 array
+        image.array[dead_mask > 0] = 0
+
+        return image
diff --git a/roman_imsim/effects/gain.py b/roman_imsim/effects/gain.py
new file mode 100644
index 00000000..c90cb55f
--- /dev/null
+++ b/roman_imsim/effects/gain.py
@@ -0,0 +1,35 @@
+import os
+import numpy as np
+import fitsio as fio
+import galsim.roman as roman
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class Gain(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for gain.")
+        image /= roman.gain
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No gain data provided; galsim.roman.gain will be applied.")
+            image /= roman.gain
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for gain.")
+        gain_map = self.df["GAIN"][:, :]  # 32x32 img
+
+        t = np.zeros((4096, 32))
+        for row in range(t.shape[0]):
+            t[row, row // 128] = 1
+        gain_expand = (t.dot(gain_map)).dot(t.T)  # 4096x4096 gain img
+        image.array[:, :] /= gain_expand
+        return image
diff --git a/roman_imsim/effects/ipc.py b/roman_imsim/effects/ipc.py
new file mode 100644
index 00000000..ebd13c70
--- /dev/null
+++ b/roman_imsim/effects/ipc.py
@@ -0,0 +1,66 @@
+import os
+import numpy as np
+import fitsio as fio
+import galsim.roman as roman
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class IPC(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for IPC effect.")
+        kernel = roman.ipc_kernel
+        image.applyIPC(kernel, edge_treatment="extend", fill_value=None)
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No IPC kernel data file provided; no IPC effect will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for IPC effect.")
+        # pad the array by one pixel at the four edges
+        num_grids = 4  # num_grids <= 8
+        grid_size = 4096 // num_grids
+
+        array_pad = image.array[4:-4, 4:-4]  # it's an array instead of img
+        array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric")  # 4098x4098 array
+
+        K = self.df["IPC"][:, :, :, :]  # 3,3,512, 512
+
+        t = np.zeros((grid_size, 512))
+        for row in range(t.shape[0]):
+            t[row, row // (grid_size // 512)] = 1
+
+        array_out = np.zeros((4096, 4096))
+        # split job in sub_grids to reduce memory
+        for gj in range(num_grids):
+            for gi in range(num_grids):
+                K_pad = np.zeros((3, 3, grid_size + 2, grid_size + 2))
+
+                for j in range(3):
+                    for i in range(3):
+                        tmp = (t.dot(K[j, i, :, :])).dot(t.T)  # grid_sizexgrid_size
+                        K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric")
+
+                for dy in range(-1, 2):
+                    for dx in range(-1, 2):
+
+                        array_out[
+                            gj * grid_size : (gj + 1) * grid_size, gi * grid_size : (gi + 1) * grid_size
+                        ] += (
+                            K_pad[1 + dy, 1 + dx, 1 - dy : 1 - dy + grid_size, 1 - dx : 1 - dx + grid_size]
+                            * array_pad[
+                                1 - dy + gj * grid_size : 1 - dy + (gj + 1) * grid_size,
+                                1 - dx + gi * grid_size : 1 - dx + (gi + 1) * grid_size,
+                            ]
+                        )
+
+        image.array[:, :] = array_out
+        return image
diff --git a/roman_imsim/effects/nonlinearity.py b/roman_imsim/effects/nonlinearity.py
new file mode 100644
index 00000000..f2f43862
--- /dev/null
+++ b/roman_imsim/effects/nonlinearity.py
@@ -0,0 +1,47 @@
+import os
+import fitsio as fio
+import galsim.roman as roman
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class Nonlinearity(RomanEffects):
+    """
+    Applying a quadratic non-linearity.
+
+    Note that users who wish to apply some other nonlinearity function (perhaps for other NIR
+    detectors, or for CCDs) can use the more general nonlinearity routine, which uses the
+    following syntax:
+    final_image.applyNonlinearity(NLfunc=NLfunc)
+    with NLfunc being a callable function that specifies how the output image pixel values
+    should relate to the input ones.
+    """
+
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.info("Galsim.roman.NLfunc will be applied for simulating non-linearity effect.")
+        image.applyNonlinearity(NLfunc=roman.NLfunc)
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning(
+                "No non-linearity data file provided; no non-linearity effect will be applied."
+            )
+            return image
+
+        self.logger.warning("Lab measured model will be applied for non-linearity effect.")
+
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        image.array[:, :] -= (
+            self.df["CNL"][0, :, :][0] * image.array**2
+            + self.df["CNL"][1, :, :][0] * image.array**3
+            + self.df["CNL"][2, :, :][0] * image.array**4
+        )
+
+        return image
diff --git a/roman_imsim/effects/persistence.py b/roman_imsim/effects/persistence.py
new file mode 100644
index 00000000..d814be15
--- /dev/null
+++ b/roman_imsim/effects/persistence.py
@@ -0,0 +1,129 @@
+import os
+import numpy as np
+import fitsio as fio
+import galsim
+from . import RomanEffects
+from .utils import sca_number_to_file, get_pointing
+
+
+class Persistence(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+        p_list = np.array([get_pointing(self.base, i, self.sca) for i in range(self.visit - 10, self.visit)])
+        dt_list = np.array([(self.pointing.date - p.date).total_seconds() for p in p_list])
+        self.p_pers = p_list[np.where((dt_list > 0) & (dt_list < self.pointing.exptime * 10))]
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for persistence effect.")
+        for p in self.p_pers:
+            dt = (
+                self.pointing.date - p.date
+            ).total_seconds() - self.pointing.exptime / 2  # avg time since end of exposures
+            # self.base['output']['file_name']['items'] = [p.filter, p.visit, p.sca]
+            # imfilename = ParseValue(self.base['output'], 'file_name', self.base, str)[0]
+            imfilename = self.base["output"]["file_name"]["format"] % (p.filter, p.visit, p.sca)
+            fn = os.path.join(self.base["output"]["dir"], imfilename)
+
+            # [TODO]
+            if not os.path.exists(fn):
+                continue
+
+            # apply all the effects that occured before persistence on the previouse exposures
+            # since max of the sky background is of order 100, it is thus negligible for persistence
+            bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096)
+            x = galsim.Image(bound_pad)
+            x.array[4:-4, 4:-4] = galsim.Image(fio.FITS(fn)[0].read()).array[:, :]
+
+            recip_failure = self.cross_refer("ReciprocityFailure")
+            x = recip_failure.apply(image=x)
+
+            x.array.clip(0)  # remove negative stimulus
+
+            image.array[:, :] += (
+                galsim.roman.roman_detectors.fermi_linear(x.array, dt) * self.pointing.exptime
+            )
+
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No persistence data file provided; no persistence effect will be applied.")
+            return image
+
+        self.logger.warning("Lab measured model will be applied for persistence effect.")
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        # setup parameters for persistence
+        Q01 = self.df["PERSIST"].read_header()["Q01"]
+        Q02 = self.df["PERSIST"].read_header()["Q02"]
+        Q03 = self.df["PERSIST"].read_header()["Q03"]
+        Q04 = self.df["PERSIST"].read_header()["Q04"]
+        Q05 = self.df["PERSIST"].read_header()["Q05"]
+        Q06 = self.df["PERSIST"].read_header()["Q06"]
+        alpha = self.df["PERSIST"].read_header()["ALPHA"]
+
+        # iterate over previous exposures
+        for p in self.p_pers:
+            dt = (
+                self.pointing.date - p.date
+            ).total_seconds() - self.pointing.exptime / 2  # avg time since end of exposures
+            # linear time dependence (approximate until we get t1 and Delat t of the data)
+            fac_dt = (self.pointing.exptime / 2.0) / dt
+            # self.base['output']['file_name']['items'] = [p.filter, p.visit, p.sca]
+            # imfilename = ParseValue(self.base['output'], 'file_name', self.base, str)[0]
+
+            imfilename = self.base["output"]["file_name"]["format"] % (p.filter, p.visit, p.sca)
+            fn = os.path.join(self.base["output"]["dir"], imfilename)
+
+            # [TODO]
+            if not os.path.exists(fn):
+                continue
+
+            # apply all the effects that occured before persistence on the previouse exposures
+            # since max of the sky background is of order 100, it is thus negligible for persistence
+            # same for brighter fatter effect
+            bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096)
+            x = galsim.Image(bound_pad)
+            x.array[4:-4, 4:-4] = galsim.Image(fio.FITS(fn)[0].read()).array[:, :]
+            # x = self.qe(x).array[:,:]
+
+            qe = self.cross_refer("QuantumEfficiency")
+            x = qe.apply(image=x).array[:, :]
+
+            x = x.clip(0.1)  # remove negative and zero stimulus
+
+            # Do linear interpolation
+            a = np.zeros(x.shape)
+            a += ((x < Q01)) * x / Q01
+            a += ((x >= Q01) & (x < Q02)) * (Q02 - x) / (Q02 - Q01)
+            image.array[:, :] += a * self.df["PERSIST"][0, :, :][0] * fac_dt
+
+            a = np.zeros(x.shape)
+            a += ((x >= Q01) & (x < Q02)) * (x - Q01) / (Q02 - Q01)
+            a += ((x >= Q02) & (x < Q03)) * (Q03 - x) / (Q03 - Q02)
+            image.array[:, :] += a * self.df["PERSIST"][1, :, :][0] * fac_dt
+
+            a = np.zeros(x.shape)
+            a += ((x >= Q02) & (x < Q03)) * (x - Q02) / (Q03 - Q02)
+            a += ((x >= Q03) & (x < Q04)) * (Q04 - x) / (Q04 - Q03)
+            image.array[:, :] += a * self.df["PERSIST"][2, :, :][0] * fac_dt
+
+            a = np.zeros(x.shape)
+            a += ((x >= Q03) & (x < Q04)) * (x - Q03) / (Q04 - Q03)
+            a += ((x >= Q04) & (x < Q05)) * (Q05 - x) / (Q05 - Q04)
+            image.array[:, :] += a * self.df["PERSIST"][3, :, :][0] * fac_dt
+
+            a = np.zeros(x.shape)
+            a += ((x >= Q04) & (x < Q05)) * (x - Q04) / (Q05 - Q04)
+            a += ((x >= Q05) & (x < Q06)) * (Q06 - x) / (Q06 - Q05)
+            image.array[:, :] += a * self.df["PERSIST"][4, :, :][0] * fac_dt
+
+            a = np.zeros(x.shape)
+            a += ((x >= Q05) & (x < Q06)) * (x - Q05) / (Q06 - Q05)
+            a += ((x >= Q06)) * (x / Q06) ** alpha  # avoid fractional power of negative values
+            image.array[:, :] += a * self.df["PERSIST"][5, :, :][0] * fac_dt
+
+        return image
diff --git a/roman_imsim/effects/quantum_efficiency.py b/roman_imsim/effects/quantum_efficiency.py
new file mode 100644
index 00000000..7c5951ac
--- /dev/null
+++ b/roman_imsim/effects/quantum_efficiency.py
@@ -0,0 +1,21 @@
+import os
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class QuantumEfficiency(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No QE data file provided; a default value of QE = 1 will be used.")
+            return image
+
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+        self.logger.warning("Lab measured model will be applied for quantum efficiency.")
+        image.array[:, :] *= self.df["RELQE1"][:, :]  # 4096x4096 array
+        return image
diff --git a/roman_imsim/effects/read_noise.py b/roman_imsim/effects/read_noise.py
new file mode 100644
index 00000000..baac8224
--- /dev/null
+++ b/roman_imsim/effects/read_noise.py
@@ -0,0 +1,34 @@
+import os
+import fitsio as fio
+import galsim
+import galsim.roman as roman
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class ReadNoise(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for read noise.")
+        self.read_noise = galsim.GaussianNoise(self.rng, sigma=roman.read_noise)
+        self.im_read = image.copy()
+        image.addNoise(self.read_noise)
+        self.im_read = image - self.im_read
+        # self.sky.addNoise(self.read_noise)
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No read noise data file provided; no read noise will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for read noise.")
+        rdn = self.df["READ"][2, :, :].flatten()  # flattened 4096x4096 array
+        self.im_read = self.rng_np.normal(loc=0.0, scale=rdn).reshape(image.array.shape).astype(image.dtype)
+        image.array[:, :] += self.im_read
+        return image
diff --git a/roman_imsim/effects/reciprocity_failure.py b/roman_imsim/effects/reciprocity_failure.py
new file mode 100644
index 00000000..ef577045
--- /dev/null
+++ b/roman_imsim/effects/reciprocity_failure.py
@@ -0,0 +1,19 @@
+from . import RomanEffects
+import galsim.roman as roman
+
+
+class ReciprocityFailure(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.alpha = self.params["alpha"] if "alpha" in self.params else roman.reciprocity_alpha
+        self.base_flux = self.params["base_flux"] if "base_flux" in self.params else 1.0
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        # Add reciprocity effect
+        self.logger.warning("Simple model will be applied for reciprocity failure effect.")
+        exptime = self.pointing.exptime
+        image.addReciprocityFailure(exp_time=exptime, alpha=self.alpha, base_flux=self.base_flux)
+        return image
diff --git a/roman_imsim/effects/roman_effects.py b/roman_imsim/effects/roman_effects.py
new file mode 100644
index 00000000..94c51037
--- /dev/null
+++ b/roman_imsim/effects/roman_effects.py
@@ -0,0 +1,131 @@
+import galsim as galsim
+import numpy as np
+from .utils import get_pointing
+import roman_imsim.effects as Effects
+
+
+class RomanEffects(object):
+    """
+    Class to simulate non-idealities and noise of roman detector images.
+    """
+
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        self.params = params
+        self.base = base
+        self.visit = int(self.base["input"]["obseq_data"]["visit"])
+        self.sca = base["image"]["SCA"]
+        self.filter = base["image"]["filter"]
+        self.sca_filepath = base["image"]["sca_filepath"]
+        self.rng_iter = rng_iter if rng_iter else self.visit * self.sca
+
+        self.rng = rng
+        self.noise = galsim.PoissonNoise(self.rng)
+        self.rng_np = np.random.default_rng(self.rng_iter)
+        self.pointing = get_pointing(self.base, self.visit, self.sca)
+        self.exptime = self.pointing.exptime
+        self.logger = logger
+
+        self.force_cvz = False
+        if "force_cvz" in self.base["image"]["wcs"]:
+            if self.base["image"]["wcs"]["force_cvz"]:
+                self.force_cvz = True
+
+        self.save_diff = False
+        if "save_diff" in self.base["image"]:
+            self.save_diff = bool(self.base["image"]["save_diff"])
+            if "diff_dir" in self.base["output"]:
+                self.diff_dir = self.base["output"]["diff_dir"]
+            else:
+                self.diff_dir = self.base["output"]["dir"]
+
+    def is_model_valid(self):
+        self.model = getattr(self, self.params["model"], None)
+        if self.model is None:
+            self.logger.warning(
+                "%s hasn't been implemented yet, the simple model will be applied for %s"
+                % (str(self.params["model"]), str(self.__class__.__name__))
+            )
+            self.model = self.simple_model
+
+    def simple_model(self, image):
+        self.logger.info("Applying the default model...")
+        return image
+
+    def cross_refer(self, effect_name):
+        if effect_name not in self.base["image"]["add_effects"]:
+            try:
+                effect = getattr(Effects, effect_name)(
+                    {"model": "simple_model"}, self.base, self.logger, self.rng
+                )
+            except Exception as e:
+                self.logger.warning(e)
+                # self.logger.warning("Effect %s is not implemented!"%(effect_name))
+                return None
+        else:
+            try:
+                params = self.base["image"]["add_effects"][effect_name]
+                effect = getattr(Effects, effect_name)(params, self.base, self.logger, self.rng)
+            except Exception as e:
+                self.logger.warning(e)
+                # self.logger.warning("Effect %s is not implemented!"%(effect_name))
+                return None
+        return effect
+
+    def apply(self, image):
+        image = self.model(image)
+        return image
+
+    def set_diff(self, im=None):
+        if self.save_diff:
+            self.pre = im.copy()
+            self.pre.write("bg.fits", dir=self.diff_dir)
+        return
+
+    def diff(self, msg, im=None, verbose=True):
+        if self.save_diff:
+            diff = im - self.pre
+            diff.write("%s_diff.fits" % msg, dir=self.diff_dir)
+            self.pre = im.copy()
+            im.write("%s_cumul.fits" % msg, dir=self.diff_dir)
+        return
+
+
+class setup_sky(object):
+    def __init__(self, base, logger, rng, rng_iter=None):
+        self.base = base
+        self.logger = logger
+        self.rng = rng
+
+        self.visit = int(self.base["input"]["obseq_data"]["visit"])
+        self.sca = base["image"]["SCA"]
+        self.pointing = get_pointing(self.base, self.visit, self.sca)
+
+    def get_sky_image(self):
+        bounds = galsim.BoundsI(xmin=1, ymin=1, xmax=4088, ymax=4088)
+        self.sky_img = galsim.Image(bounds=bounds, wcs=self.pointing.WCS)
+        self.pointing.WCS.makeSkyImage(self.sky_img, 0.0)
+
+        bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096)
+        im_pad = galsim.Image(bound_pad)
+        im_pad.array[4:-4, 4:-4] = self.sky_img.array[:, :]
+
+        effects_list = list(self.base["image"]["add_effects"])
+        if "Background" not in effects_list:
+            return self.sky_img
+
+        bkg_idx = effects_list.index("Background")
+        for i in range(bkg_idx, len(effects_list)):
+            effect_name = effects_list[i]
+            args = (self.base["image"]["add_effects"][effect_name], self.base, self.logger, self.rng)
+            effect = getattr(Effects, effect_name)(*args)
+            im_pad = effect.apply(image=im_pad)
+
+        im_pad.quantize()
+        # output 4088x4088 img in uint16
+        self.sky_img = im_pad.array[4:-4, 4:-4]
+        self.sky_img = galsim.Image(self.sky_img, dtype=np.uint16)
+
+        return self.sky_img
+
+    def save_sky_img(self, outdir=".", sky_img_name="sky_img.fits"):
+        self.sky_img.write(sky_img_name, dir=outdir)
diff --git a/roman_imsim/effects/saturation.py b/roman_imsim/effects/saturation.py
new file mode 100644
index 00000000..b2d1ef4d
--- /dev/null
+++ b/roman_imsim/effects/saturation.py
@@ -0,0 +1,36 @@
+import os
+import numpy as np
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class Saturation(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+        self.saturation_level = (
+            self.params["saturation_level"] if "saturation_level" in self.params else 100000
+        )
+
+    def simple_model(self, image):
+        self.logger.warning("Simple model will be applied for saturation.")
+        saturation_array = np.ones_like(image.array) * self.saturation_level
+        where_sat = np.where(image.array > saturation_array)
+        image.array[where_sat] = saturation_array[where_sat]
+        return image
+
+    def lab_model(self, image):
+        if self.sca_filepath is None:
+            self.logger.warning("No saturation data file provided; no saturation effect will be applied.")
+            return image
+
+        self.logger.warning("Lab measured model will be applied for saturation effect.")
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        saturation_array = self.df["SATURATE"][:, :]  # 4096x4096 array
+        where_sat = np.where(image.array > saturation_array)
+        image.array[where_sat] = saturation_array[where_sat]
+        return image
diff --git a/roman_imsim/effects/utils.py b/roman_imsim/effects/utils.py
new file mode 100644
index 00000000..a350f76a
--- /dev/null
+++ b/roman_imsim/effects/utils.py
@@ -0,0 +1,59 @@
+import galsim as galsim
+import galsim.roman as roman
+import numpy as np
+from roman_imsim.obseq import ObSeqDataLoader
+
+sca_number_to_file = {
+    1: "SCA_22066_211227_v001.fits",
+    2: "SCA_21815_211221_v001.fits",
+    3: "SCA_21946_211225_v001.fits",
+    4: "SCA_22073_211229_v001.fits",
+    5: "SCA_21816_211222_v001.fits",
+    6: "SCA_20663_211102_v001.fits",
+    7: "SCA_22069_211228_v001.fits",
+    8: "SCA_21641_211216_v001.fits",
+    9: "SCA_21813_211219_v001.fits",
+    10: "SCA_22078_211230_v001.fits",
+    11: "SCA_21947_211226_v001.fits",
+    12: "SCA_22077_211230_v001.fits",
+    13: "SCA_22067_211227_v001.fits",
+    14: "SCA_21814_211220_v001.fits",
+    15: "SCA_21645_211228_v001.fits",
+    16: "SCA_21643_211218_v001.fits",
+    17: "SCA_21319_211211_v001.fits",
+    18: "SCA_20833_211116_v001.fits",
+}
+
+
+class get_pointing(object):
+    """
+    Class to store stuff about the telescope
+    """
+
+    def __init__(self, params, visit, SCA):
+
+        self.params = params
+        file_name = params["input"]["obseq_data"]["file_name"]
+        obseq_data = ObSeqDataLoader(file_name, visit, SCA, logger=None)
+        self.filter = obseq_data.ob["filter"]
+        self.sca = obseq_data.ob["sca"]
+        self.visit = obseq_data.ob["visit"]
+        self.date = obseq_data.ob["date"]
+        self.exptime = obseq_data.ob["exptime"]
+        self.bpass = roman.getBandpasses()[self.filter]
+        self.WCS = roman.getWCS(
+            world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]),
+            PA=obseq_data.ob["pa"],
+            date=self.date,
+            SCAs=self.sca,
+            PA_is_FPA=True,
+        )[self.sca]
+        self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)))
+
+
+def translate_cvz(orig_radec, field_ra=9.5, field_dec=-44, cvz_ra=61.24, cvz_dec=-48.42):
+    ra = orig_radec.ra / galsim.degrees - field_ra
+    dec = orig_radec.dec / galsim.degrees - field_dec
+    ra += cvz_ra / np.cos(cvz_dec * np.pi / 180)
+    dec += cvz_dec
+    return galsim.CelestialCoord(ra * galsim.degrees, dec * galsim.degrees)
diff --git a/roman_imsim/effects/vtpe.py b/roman_imsim/effects/vtpe.py
new file mode 100644
index 00000000..92e1e1c8
--- /dev/null
+++ b/roman_imsim/effects/vtpe.py
@@ -0,0 +1,55 @@
+import os
+import numpy as np
+import fitsio as fio
+from . import RomanEffects
+from .utils import sca_number_to_file
+
+
+class VTPE(RomanEffects):
+    def __init__(self, params, base, logger, rng, rng_iter=None):
+        super().__init__(params, base, logger, rng, rng_iter)
+
+        self.is_model_valid()
+
+    def simple_model(self, image):
+        self.logger.warning("No vertical trailing pixel effect will be applied.")
+        return image
+
+    def lab_model(self, image):
+        """
+        Apply vertical trailing pixel effect.
+        The vertical trailing pixel effect (VTPE) is a non-linear effect that is
+        related to readout patterns.
+        Q'[j,i] = Q[j,i] + f(  Q[j,i] - Q[j-1, i]  ),
+        where f( dQ ) = dQ ( a + b * ln(1 + |dQ|/dQ0) )
+        Input
+        im           : image
+        VTPE[0,512,512]  : coefficient a binned in 8x8
+        VTPE[1,512,512]  : coefficient a
+        VTPE[2,512,512]  : coefficient dQ0
+        """
+
+        if self.sca_filepath is None:
+            self.logger.warning("No VTPE data file provided; no VTPE will be applied.")
+            return image
+        self.df = fio.FITS(os.path.join(self.sca_filepath, sca_number_to_file[self.sca]))
+
+        self.logger.warning("Lab measured model will be applied for VTPE.")
+
+        # expand 512x512 arrays to 4096x4096
+        t = np.zeros((4096, 512))
+        for row in range(t.shape[0]):
+            t[row, row // 8] = 1
+        a_vtpe = t.dot(self.df["VTPE"][0, :, :][0]).dot(t.T)
+        # NaN check
+        if np.isnan(a_vtpe).any():
+            self.logger.warning("vtpe skipped due to NaN in file")
+            return image
+        b_vtpe = t.dot(self.df["VTPE"][1, :, :][0]).dot(t.T)
+        dQ0 = t.dot(self.df["VTPE"][2, :, :][0]).dot(t.T)
+
+        dQ = image.array - np.roll(image.array, 1, axis=0)
+        dQ[0, :] *= 0
+
+        image.array[:, :] += dQ * (a_vtpe + b_vtpe * np.log(1.0 + np.abs(dQ) / dQ0))
+        return image
diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py
index 7764770d..9fbe5ef9 100644
--- a/roman_imsim/sca.py
+++ b/roman_imsim/sca.py
@@ -8,6 +8,9 @@
 from galsim.image import Image
 
 
+import roman_imsim.effects as RomanEffects
+
+
 class RomanSCAImageBuilder(ScatteredImageBuilder):
 
     def setup(self, config, base, image_num, obj_num, ignore, logger):
@@ -52,16 +55,15 @@ def setup(self, config, base, image_num, obj_num, ignore, logger):
         req = {"SCA": int, "filter": str, "mjd": float, "exptime": float}
         opt = {
             "draw_method": str,
-            "stray_light": bool,
-            "thermal_background": bool,
-            "reciprocity_failure": bool,
-            "dark_current": bool,
-            "nonlinearity": bool,
-            "ipc": bool,
-            "read_noise": bool,
+            "add_effects": dict,
+            "sca_filepath": str,
             "sky_subtract": bool,
             "ignore_noise": bool,
+            "save_diff": bool,
         }
+
+        logger.warning("opt dict = %s" % (str(opt)))
+
         params = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore + extra_ignore)[0]
 
         self.sca = params["SCA"]
@@ -71,14 +73,6 @@ def setup(self, config, base, image_num, obj_num, ignore, logger):
         self.exptime = params["exptime"]
 
         self.ignore_noise = params.get("ignore_noise", False)
-        # self.exptime = params.get('exptime', roman.exptime)  # Default is roman standard exposure time.
-        self.stray_light = params.get("stray_light", False)
-        self.thermal_background = params.get("thermal_background", False)
-        self.reciprocity_failure = params.get("reciprocity_failure", False)
-        self.dark_current = params.get("dark_current", False)
-        self.nonlinearity = params.get("nonlinearity", False)
-        self.ipc = params.get("ipc", False)
-        self.read_noise = params.get("read_noise", False)
         self.sky_subtract = params.get("sky_subtract", False)
 
         # If draw_method isn't in image field, it may be in stamp.  Check.
@@ -95,10 +89,18 @@ def setup(self, config, base, image_num, obj_num, ignore, logger):
         # # GalSim expects a wcs in the image field.
         # config['wcs'] = wcs
 
+        self.rng = galsim.config.GetRNG(config, base)
+        self.visit = int(base["input"]["obseq_data"]["visit"])
+
+        self.sca_filepath = params.get("sca_filepath", None)
+
         # If user hasn't overridden the bandpass to use, get the standard one.
         if "bandpass" not in config:
             base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger)
 
+        self.base = base
+        self.logger = logger
+
         return roman.n_pix, roman.n_pix
 
     # def getBandpass(self, filter_name):
@@ -218,93 +220,30 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, logger)
             return
 
         base["current_noise_image"] = base["current_image"]
-        wcs = base["wcs"]
-        bp = base["bandpass"]
-        rng = galsim.config.GetRNG(config, base)
+        # rng = galsim.config.GetRNG(config, base)
         logger.info("image %d: Start RomanSCA detector effects", base.get("image_num", 0))
 
-        # Things that will eventually be subtracted (if sky_subtract) will have their expectation
-        # value added to sky_image.  So technically, this includes things that aren't just sky.
-        # E.g. includes dark_current and thermal backgrounds.
-        sky_image = image.copy()
-        sky_level = roman.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center))
-        logger.debug("Adding sky_level = %s", sky_level)
-        if self.stray_light:
-            logger.debug("Stray light fraction = %s", roman.stray_light_fraction)
-            sky_level *= 1.0 + roman.stray_light_fraction
-        wcs.makeSkyImage(sky_image, sky_level)
-
-        # The other background is the expected thermal backgrounds in this band.
-        # These are provided in e-/pix/s, so we have to multiply by the exposure time.
-        if self.thermal_background:
-            tb = roman.thermal_backgrounds[self.filter] * self.exptime
-            logger.debug("Adding thermal background: %s", tb)
-            sky_image += roman.thermal_backgrounds[self.filter] * self.exptime
-
-        # The image up to here is an expectation value.
-        # Realize it as an integer number of photons.
-        poisson_noise = galsim.noise.PoissonNoise(rng)
-        if self.draw_method == "phot":
-            logger.debug("Adding poisson noise to sky photons")
-            sky_image1 = sky_image.copy()
-            sky_image1.addNoise(poisson_noise)
-            image.quantize()  # In case any profiles used InterpolatedImage, in which case
-            # the image won't necessarily be integers.
-            image += sky_image1
-        else:
-            logger.debug("Adding poisson noise")
-            image += sky_image
-            image.addNoise(poisson_noise)
-
-        # Apply the detector effects here.  Not all of these are "noise" per se, but they
-        # happen interspersed with various noise effects, so apply them all in this step.
-
-        # Note: according to Gregory Mosby & Bernard J. Rauscher, the following effects all
-        # happen "simultaneously" in the photo diodes: dark current, persistence,
-        # reciprocity failure (aka CRNL), burn in, and nonlinearity (aka CNL).
-        # Right now, we just do them in some order, but this could potentially be improved.
-        # The order we chose is historical, matching previous recommendations, but Mosby and
-        # Rauscher don't seem to think those recommendations are well-motivated.
-
-        # TODO: Add burn-in and persistence here.
-
-        if self.reciprocity_failure:
-            logger.debug("Applying reciprocity failure")
-            roman.addReciprocityFailure(image)
-
-        if self.dark_current:
-            dc = roman.dark_current * self.exptime
-            logger.debug("Adding dark current: %s", dc)
-            sky_image += dc
-            dark_noise = galsim.noise.DeviateNoise(galsim.random.PoissonDeviate(rng, dc))
-            image.addNoise(dark_noise)
-
-        if self.nonlinearity:
-            logger.debug("Applying classical nonlinearity")
-            roman.applyNonlinearity(image)
-
-        # Mosby and Rauscher say there are two read noises.  One happens before IPC, the other
-        # one after.
-        # TODO: Add read_noise1
-        if self.ipc:
-            logger.debug("Applying IPC")
-            roman.applyIPC(image)
-
-        if self.read_noise:
-            logger.debug("Adding read noise %s", roman.read_noise)
-            image.addNoise(galsim.GaussianNoise(rng, sigma=roman.read_noise))
-
-        logger.debug("Applying gain %s", roman.gain)
-        image /= roman.gain
-
-        # Make integer ADU now.
-        image.quantize()
+        # create padded image
+        bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096)
+        im_pad = galsim.Image(bound_pad)
+        im_pad.array[4:-4, 4:-4] = image.array[:, :]
+
+        effects_list = self.base["image"]["add_effects"].keys()
+        for effect_name in effects_list:
+            args = (self.base["image"]["add_effects"][effect_name], self.base, self.logger, self.rng)
+            effect = getattr(RomanEffects, effect_name)(*args)
+            im_pad = effect.apply(image=im_pad)
+
+        im_pad.quantize()
+        # output 4088x4088 img in uint16
+        image.array[:, :] = im_pad.array[4:-4, 4:-4]
 
         if self.sky_subtract:
             logger.debug("Subtracting sky image")
-            sky_image /= roman.gain
-            sky_image.quantize()
+            sky = RomanEffects.setup_sky(self.base, self.logger, self.rng)
+            sky_image = sky.get_sky_image()
             image -= sky_image
+            sky.save_sky_img(outdir=self.base["output"]["dir"])
 
 
 # Register this as a valid type