Repository: spacetelescope/jwst

Issue:
MIRI Subarrays 390 Hz EMI
_Issue [JP-3248](https://jira.stsci.edu/browse/JP-3248) was created on JIRA by [Michael Regan](https://jira.stsci.edu/secure/ViewProfile.jspa?name=mregan):_

The majority of the MIRI subarrays have an 390 Hz EMI noise pattern in the raw data. This 390 Hz maps exactly to 256 pixels and has a constant amplitude of +- 4 DN. The effect of this is to imprint this into the rate images. For short integrations in LRSSLITLESS the correlated noise from this is quite apparent in the rate images. For longer integrations the net effect is to increase the read noise by ~20%.

The long term plan is a change to the sizes and locations of the subarrays to get the frame times to be in phase with the 390 Hz like the full frame images. For the previous and near term observations this can be fixed by adding a new step to the pipeline. A prototype has been implemented and shown to remove the 390 Hz to a fraction of a DN.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Patch coverage: **`13.79%`** and project coverage change: **`-0.05%`** :warning:
> Comparison is base [(`9b22d5b`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/9b22d5b05094fcb99fd2604b1e0cbfd28ccea78f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 76.39% compared to head [(`1d1d26e`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 76.35%.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #7857      +/-   ##
==========================================
- Coverage   76.39%   76.35%   -0.05%     
==========================================
  Files         459      461       +2     
  Lines       37074    37103      +29     
==========================================
+ Hits        28324    28329       +5     
- Misses       8750     8774      +24     
```

| [Flag](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | | *Carryforward flag |
|---|---|---|---|
| [nightly](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | `77.37% <ø> (ø)` | | Carriedforward from [9b22d5b](https://app.codecov.io/gh/spacetelescope/jwst/commit/9b22d5b05094fcb99fd2604b1e0cbfd28ccea78f?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

*This pull request uses carry forward flags. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to find out more.

| [Files Changed](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [jwst/emicorr/emicorr\_step.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9lbWljb3JyL2VtaWNvcnJfc3RlcC5weQ==) | `11.11% <11.11%> (ø)` | |
| [jwst/emicorr/\_\_init\_\_.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9lbWljb3JyL19faW5pdF9fLnB5) | `50.00% <50.00%> (ø)` | |

... and [1 file with indirect coverage changes](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/7857?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

penaguerrero: Need to request a new pmap for CRDS so that the code knows about any potential entry for the new step.

penaguerrero: The docs now build but I had to remove the % and type the word instead... I don't like that solution tho.
hbushouse: I committed more reformatting of that rst file and the percent sign now works OK. I think it was really the bad formatting of the numbered list that was screwing it up. See https://jwst-pipeline--7857.org.readthedocs.build/en/7857/jwst/emicorr/description.html
hbushouse: Started a regtest run at https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1097
hbushouse: Everything here looks good, but seeing a number of errors in the regtest run. Will need to wait till they finish to investigate. Hopefully they're unrelated.
hbushouse: Rats. All the errors in the regtest run are due to the use of CRDS-OPS, which doesn't yet have definitions for the new "emicorr" reference file type, hence all the calls to retrieve the ref file throw an error. Need to switch to using CRDS_TEST. Will start a new regtest using that.
penaguerrero: the docs look good except for the formatting of the reference file description, I think i need to add a white line so that the ASDF tree looks like we want
penaguerrero: ok, the docs built and the tree looks ok but not pretty: https://jwst-pipeline--7857.org.readthedocs.build/en/7857/jwst/emicorr/reference_files.html

hbushouse: > ok, the docs built and the tree looks ok but not pretty: https://jwst-pipeline--7857.org.readthedocs.build/en/7857/jwst/emicorr/reference_files.html

You could try enclosing the tree text in a literal block. To do that, end the preceding paragraph with "::" and then have your ascii text tree indented, with a blank line before and after, e.g.
```
This is my intro text::

  This block of text, which is indented,
  should get rendered as a literal, with no formatting.

  Keep indenting each line that you want included in the
  literal block.

Then revert to normal text.
```
penaguerrero: yay! that worked! https://jwst-pipeline--7857.org.readthedocs.build/en/7857/jwst/emicorr/reference_files.html
thanks!
hbushouse: Still getting errors in regtests, even when using emicorr ref file from CRDS_TEST. Here's what I've found:
```
test_miri_image_detector1, in emicorr.py:

                        if colstop == 258:
>                           times_this_int[k, j, nx4-1] = times_this_int[k, j, nx4-1] + ref_pix_sample.astype('ulonglong')
E                           AttributeError: 'int' object has no attribute 'astype'
```

test_miri_dark is complaining that it can't find an emicorr ref file in CRDS. This might be due to the test data file having a DATE-OBS that's older than 2022-04-01, which is the USEAFTER on the emicorr ref file. Which, in my opinion is stupid. Why not have the USEAFTER go all the way back to cover all possible existing MIRI data. Surely the EMI problem did not just start in April 2022.
hbushouse: @penaguerrero Please try running jwst/regtest/test_miri_image.py locally to debug why it's crashing.
penaguerrero: ok! Reg tests ran with some failures but not associated to this PR: https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1099/#showFailuresLink

hbushouse: > ok! Reg tests ran with some failures but not associated to this PR: https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1099/#showFailuresLink

This run is again showing all of the "E           crds.core.exceptions.CrdsLookupError: Error determining best reference for 'emicorr'  =   Unknown reference type 'emicorr'" errors, even for MIRI datasets, presumably because it was run using CRDS Ops, but the new emicorr ref file only exists in CRDS_TEST. Need to run again with "CRDS_PATH=/grp/crds/jwst/test" set in the "ENV_VARS" box that comes up after you click on "Build with Parameters" to start the run.
hbushouse: Woot!  Regtest run  shows differences (but no crashes) in the MIRI tests that use subarrays due to new keywords:
```
     Extra keyword 'R_MIREMI' in a: 'crds://jwst_miri_emicorr_0001.asdf'
     Extra keyword 'S_MIREMI' in a: 'COMPLETE'
```
and slight changes to pixel values, so that looks like the step is running!

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index f734aa2c3..3e604d1e9 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -24,6 +24,11 @@ documentation
   with ``tweakreg`` step and ``transfer_wcs_correction`` function. Fixed
   reported typos in ``tweakreg`` documentation. [#8084]
 
+emicorr
+----------
+
+- Added new step for removing EMI from all MIRI data. [#7857]
+
 extract_1d
 ----------
 
@@ -99,6 +104,8 @@ pipeline
   ``msa_flagging`` immediately follows ``assign_wcs``, so that both steps have
   been applied before calling ``nsclean``. [#8000]
 
+- Added emicorr step to calwebb_detector1 and calwebb_dark. [#7857]
+
 pixel_replace
 -------------
 
diff --git a/docs/jwst/emicorr/arguments.rst b/docs/jwst/emicorr/arguments.rst
new file mode 100644
index 000000000..9dedea0ec
--- /dev/null
+++ b/docs/jwst/emicorr/arguments.rst
@@ -0,0 +1,25 @@
+Step Arguments
+==============
+The ``emicorr`` step has the following step-specific arguments.
+
+``--nints_to_phase`` (integer, default=None)
+    Number of integrations to phase.
+
+``--nbins`` (integer, default=None)
+    Number of bins in one phased wave.
+
+``--scale_reference`` (boolean, default=True)
+    If True, the reference wavelength will be scaled to the
+    data's phase amplitude.
+
+``--user_supplied_reffile`` (boolean, default=None)
+    This is to specify an ASDF-format user-created reference file.
+
+``--save_intermediate_results``  (string, default=False)
+    This is a boolean flag to specify whether to write a step output
+    file with the EMI correction. Additionaly, if the flag
+    ``user_supplied_reffile`` is None and no CRDS reference file was
+    found, all the on-the-fly frequencies phase amplitudes will be
+    saved to an ASDF output with the same format as an EMI reference
+    file.
+
diff --git a/docs/jwst/emicorr/description.rst b/docs/jwst/emicorr/description.rst
new file mode 100644
index 000000000..cfb181334
--- /dev/null
+++ b/docs/jwst/emicorr/description.rst
@@ -0,0 +1,70 @@
+Description
+===========
+
+:Class: `jwst.emicorr.EmiCorrStep`
+:Alias: emicorr
+
+Overview
+--------
+The ``emicorr`` step corrects for known noise patterns in the raw MIRI data.
+The majority of the MIRI subarrays have an 390 Hz or other electromagnetic
+interference (EMI) noise pattern in the raw data. The known frequencies to
+correct for are in the EMI reference file, under the key ``frequencies``.
+The effect of these EMI frequencies is to imprint each into the rate
+images. For short integrations in ``LRSSLITLESS`` the correlated noise from
+this is quite apparent in the rate images. For longer integrations the net
+effect is to increase the read noise by about 20\%.
+
+The process to do the correction is the following (repeated
+recursively for each discrete EMI frequency desired):
+
+#. Read image data.
+
+#. Make very crude slope image and fixed pattern "super" bias for each
+   integration, ignoring everything (nonlin, saturation, badpix, etc).
+
+#. Subtract scaled slope image and bias from each frame of each integration.
+
+#. Calculate phase of every pixel in the image at the desired EMI frequency
+   (e.g. 390 Hz) relative to the first pixel in the image.
+
+#. Make a binned, phased amplitude (pa) wave from the cleaned data (plot
+   cleaned vs phase, then bin by phase).
+
+#. Measure the phase shift between the binned pa wave and the input
+   reference wave. [#f1]_
+
+#. Use look-up table to get the aligned reference wave value for each pixel
+   (the noise array corresponding to the input image). [#f1]_
+
+#. Subtract the noise array from the input image and return the cleaned result.
+
+.. [#f1] Alternately, use the binned pa wave instead of the reference wave to "self-correct".
+
+The long term plan is a change to the sizes and locations of the subarrays
+to get the frame times to be in phase with the known noise frequencies like
+the full frame images. For the previous and near term observations this can
+be fixed through application of the ``emicorr`` step.
+
+An EMICORR reference file can be used to correct for all known noise
+patterns. The reference file is expected to be in ASDF format, containing
+the exact frequency numbers, the corresponding 1D array for the phase
+amplitudes, and a ``subarray_cases`` dictionary that contains
+the frequencies to correct for according to subarray, read pattern, and
+detector. If there is no reference file found in CRDS, the step has a set
+of default frequencies and subarray cases for which the correction is
+applied.
+
+Input
+-----
+The input file is the ``_uncal`` file after the
+:ref:`dq_init_step <dq_init_step>` step has been
+applied, in the in the :ref:`calwebb_detector1 <calwebb_detector1>`
+pipeline.
+
+Output
+------
+The output will be a partially-processed ``RampModel`` with the
+corrected data in the SCI extension, meaning, the effect of the
+EMI frequencies (either the default values or the ones in the
+reference file) removed. All other extensions will remain the same.
diff --git a/docs/jwst/emicorr/index.rst b/docs/jwst/emicorr/index.rst
new file mode 100644
index 000000000..a7c5ff4f4
--- /dev/null
+++ b/docs/jwst/emicorr/index.rst
@@ -0,0 +1,14 @@
+.. _emicorr_step:
+
+===================
+MIRI EMI Correction
+===================
+
+.. toctree::
+   :maxdepth: 2
+
+   description.rst
+   arguments.rst
+   reference_files.rst
+
+.. automodapi:: jwst.emicorr
diff --git a/docs/jwst/emicorr/reference_files.rst b/docs/jwst/emicorr/reference_files.rst
new file mode 100644
index 000000000..0c8802a54
--- /dev/null
+++ b/docs/jwst/emicorr/reference_files.rst
@@ -0,0 +1,6 @@
+Reference Files
+===============
+
+The ``emicorr`` step can use an EMI reference file.
+
+.. include:: ../references_general/emi_reffile.inc
diff --git a/docs/jwst/package_index.rst b/docs/jwst/package_index.rst
index ebace4b1f..9a7521124 100644
--- a/docs/jwst/package_index.rst
+++ b/docs/jwst/package_index.rst
@@ -19,6 +19,7 @@ Package Index
    cube_build/index.rst
    dark_current/index.rst
    dq_init/index.rst
+   emicorr/index.rst
    exp_to_source/index.rst
    extract_1d/index.rst
    extract_2d/index.rst
diff --git a/docs/jwst/pipeline/calwebb_dark.rst b/docs/jwst/pipeline/calwebb_dark.rst
index 6be2c1ad9..65b51c867 100644
--- a/docs/jwst/pipeline/calwebb_dark.rst
+++ b/docs/jwst/pipeline/calwebb_dark.rst
@@ -19,6 +19,8 @@ pipeline, the order of steps is a bit different for MIRI exposures.
 +---------------------------------------+-----------------------------------------+
 | :ref:`dq_init <dq_init_step>`         | :ref:`dq_init <dq_init_step>`           |
 +---------------------------------------+-----------------------------------------+
+|                                       | :ref:`emicorr <emicorr_step>`           |
++---------------------------------------+-----------------------------------------+
 | :ref:`saturation <saturation_step>`   | :ref:`saturation <saturation_step>`     |
 +---------------------------------------+-----------------------------------------+
 | :ref:`ipc <ipc_step>` [1]_            | :ref:`ipc <ipc_step>`                   |
diff --git a/docs/jwst/pipeline/calwebb_detector1.rst b/docs/jwst/pipeline/calwebb_detector1.rst
index 0060592df..e66d135af 100644
--- a/docs/jwst/pipeline/calwebb_detector1.rst
+++ b/docs/jwst/pipeline/calwebb_detector1.rst
@@ -54,6 +54,8 @@ on either the first group or frame zero pixel values.
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 | :ref:`dq_init <dq_init_step>`                                        | |check| | |check| | :ref:`dq_init <dq_init_step>`           | |check| | |check| |
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
+|                                                                      |         |         | :ref:`emicorr <emicorr_step>`           | |check| | |check| |
++----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 | :ref:`saturation <saturation_step>`                                  | |check| | |check| | :ref:`saturation <saturation_step>`     | |check| | |check| |
 +----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+
 | :ref:`ipc <ipc_step>` [1]_                                           |         |         | :ref:`ipc <ipc_step>`                   |         |         |
diff --git a/docs/jwst/references_general/emi_reffile.inc b/docs/jwst/references_general/emi_reffile.inc
new file mode 100644
index 000000000..5aa9be601
--- /dev/null
+++ b/docs/jwst/references_general/emi_reffile.inc
@@ -0,0 +1,70 @@
+.. _emi_reffile:
+
+EMI Reference File
+------------------
+
+:REFTYPE: EMI
+:Data model: `~jwst.datamodels.EmiModel`
+
+The EMI reference file contains data necessary for removing
+contaminating MIRI EMI frequencies.
+
+.. include:: ../references_general/emi_selection.inc
+
+.. include:: ../includes/standard_keywords.inc
+
+EMI Reference File Format
++++++++++++++++++++++++++
+
+MIRI EMI reference files are in ASDF format.  The EMI
+reference file contains the frequencies for which the image will
+be corrected.
+Example file contents::
+
+  root (AsdfObject)
+    ├─frequencies (dict)
+      │ ├─Hz10 (dict)
+      │ │ ├─frequency (float): 10.039216
+      │ │ └─phase_amplitudes (NDArrayType): shape=(500,), dtype=float32
+      │ ├─Hz10_slow_MIRIFULONG (dict)
+      │ │ ├─frequency (float): 10.039216
+      │ │ └─phase_amplitudes (NDArrayType): shape=(250,), dtype=float32
+    └─subarray_cases (dict)
+      ├─BRIGHTSKY (dict)
+      │ ├─frameclocks (int): 23968
+      │ ├─freqs (dict)
+      │ │ ├─FASTR1 (list)
+      │ │ │ └─[0] (str): Hz10
+      │ │ └─SLOWR1 (dict)
+      │ │   ├─MIRIFULONG (list)
+      │ │   │ └─[0] (str): Hz10_slow_MIRIFULONG
+      │ │   ├─MIRIFUSHORT (list)
+      │ │   │ └─[0] (str): Hz10_slow_MIRIFUSHORT
+      │ │   └─MIRIMAGE (list)
+      │ │     └─[0] (str): Hz10_slow_MIRIMAGE
+      │ └─rowclocks (int): 82
+      ├─MASK1140 (dict)
+      │ ├─frameclocks (int): 23968
+      │ ├─freqs (dict)
+      │ │ ├─FASTR1 (list)
+      │ │ │ ├─[0] (str): Hz390
+      │ │ │ └─[1] (str): Hz10
+      │ │ └─SLOWR1 (dict)
+      │ │   ├─MIRIFULONG (list)
+      │ │   │ ├─[0] (str): Hz390
+      │ │   │ └─[1] (str): Hz10_slow_MIRIFULONG
+      │ │   ├─MIRIFUSHORT (list)
+      │ │   │ ├─[0] (str): Hz390
+      │ │   │ └─[1] (str): Hz10_slow_MIRIFUSHORT
+      │ │   └─MIRIMAGE (list)
+      │ │     ├─[0] (str): Hz390
+      │ └─rowclocks (int): 82
+
+
+Frequency Selection
+^^^^^^^^^^^^^^^^^^^
+
+The frequency to be corrected for will be selected according to the
+dictionary contained in the key ``subarray_cases`` in the reference
+file. This contains the subarray names and the names of the
+corresponding frequencies to be used in the correction.
diff --git a/docs/jwst/references_general/emi_selection.inc b/docs/jwst/references_general/emi_selection.inc
new file mode 100644
index 000000000..9d029cff2
--- /dev/null
+++ b/docs/jwst/references_general/emi_selection.inc
@@ -0,0 +1,13 @@
+Reference Selection Keywords for EMICORR
+++++++++++++++++++++++++++++++++++++++++
+
+CRDS selects appropriate EMICORR references based on the
+following keywords.  EMICORR is not applicable for instruments
+not in the table.  All keywords used for file selection are
+*required*.
+
+========== ============================
+Instrument Keywords
+========== ============================
+MIRI       INSTRUME, DATE-OBS
+========== ============================
diff --git a/jwst/emicorr/__init__.py b/jwst/emicorr/__init__.py
new file mode 100644
index 000000000..4ad0ccef5
--- /dev/null
+++ b/jwst/emicorr/__init__.py
@@ -0,0 +1,3 @@
+from .emicorr_step import EmiCorrStep
+
+__all__ = ['EmiCorrStep']
diff --git a/jwst/emicorr/emicorr.py b/jwst/emicorr/emicorr.py
new file mode 100644
index 000000000..f29564156
--- /dev/null
+++ b/jwst/emicorr/emicorr.py
@@ -0,0 +1,795 @@
+#
+#  Module for applying MIRI EMI correction
+#
+
+import numpy as np
+import logging
+
+from stdatamodels.jwst import datamodels
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+
+default_emi_freqs = {
+    "Hz390": 390.625,
+    "Hz390_sub128": 390.625,
+    "Hz10": 10.039216,
+    "Hz10_slow_MIRIMAGE": 10.039216,
+    "Hz10_slow_MIRIFULONG": 10.039216,
+    "Hz10_slow_MIRIFUSHORT": 10.039216
+}
+
+default_subarray_cases = {
+
+    # 390Hz out-of-phase - these all need 390hz correction
+
+    "SLITLESSPRISM": {
+        "rowclocks": 28,
+        "frameclocks": 15904,
+        "freqs":  {"FASTR1": ["Hz390", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    "MASKLYOT": {
+        "rowclocks": 90,
+        "frameclocks": 32400,
+        "freqs":  {"FASTR1": ["Hz390", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    "SUB64": {
+        "rowclocks": 28,
+        "frameclocks": 8512,
+        "freqs":  {"FASTR1": ["Hz390", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    "SUB128": {
+        "rowclocks": 44,
+        "frameclocks": 11904,
+        "freqs":  {"FASTR1": ["Hz390_sub128", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    "MASK1140": {
+        "rowclocks": 82,
+        "frameclocks": 23968,
+        "freqs":  {"FASTR1": ["Hz390", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    "MASK1550": {
+        "rowclocks": 82,
+        "frameclocks": 23968,
+        "freqs":  {"FASTR1": ["Hz390", "Hz10"],
+                   "SLOWR1":  {"MIRIMAGE" : ["Hz390", "Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz390", "Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz390", "Hz10_slow_MIRIFUSHORT"]}}},
+
+    # 390Hz already in-phase for these, but may need corr for other
+    # frequencies (e.g. 10Hz heater noise)
+
+    "FULL_FASTR1": {
+        "rowclocks": 271,
+        "frameclocks": 277504,
+        "freqs": {"FASTR1" : ["Hz10"]}},
+
+    "FULL_SLOWR1": {
+        "rowclocks": 2333,
+        "frameclocks": 2388992,
+        "freqs": {"SLOWR1":  {"MIRIMAGE" : ["Hz10_slow_MIRIMAGE"],
+                              "MIRIFULONG" : ["Hz10_slow_MIRIFULONG"],
+                              "MIRIFUSHORT" : ["Hz10_slow_MIRIFUSHORT"]}}},
+
+    "BRIGHTSKY": {
+        "rowclocks": 162,
+        "frameclocks": 86528,
+        "freqs": {"FASTR1" : ["Hz10"],
+                  "SLOWR1": {"MIRIMAGE" : ["Hz10_slow_MIRIMAGE"],
+                            "MIRIFULONG" : ["Hz10_slow_MIRIFULONG"],
+                            "MIRIFUSHORT" : ["Hz10_slow_MIRIFUSHORT"]}}},
+
+    "SUB256": {
+        "rowclocks": 96,
+        "frameclocks": 29952,
+        "freqs": {"FASTR1" : ["Hz10"],
+                  "SLOWR1": {"MIRIMAGE" : ["Hz10_slow_MIRIMAGE"],
+                            "MIRIFULONG" : ["Hz10_slow_MIRIFULONG"],
+                            "MIRIFUSHORT" : ["Hz10_slow_MIRIFUSHORT"]}}},
+}
+
+
+def do_correction(input_model, emicorr_model, save_onthefly_reffile, **pars):
+    """
+    EMI-correct a JWST data model using an emicorr model
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodels.JwstDataModel`
+        Input science data model to be emi-corrected
+
+    emicorr_model : `~jwst.datamodels.EmiModel`
+        Data model containing emi correction
+
+    save_onthefly_reffile : str or None
+        Full path and root name of on-the-fly reference file
+
+    pars : dict
+        Optional user-specified parameters to modify how emicorr
+        will operate.  Valid parameters include:
+            save_intermediate_results - saves the output into a file and the
+                                        reference file (if created on-the-fly)
+            user_supplied_reffile - reference file supplied by the user
+
+    Returns
+    -------
+    output_model : JWST data model
+        emi-corrected science data model
+
+    """
+    save_intermediate_results = pars['save_intermediate_results']
+    user_supplied_reffile = pars['user_supplied_reffile']
+    nints_to_phase = pars['nints_to_phase']
+    nbins = pars['nbins']
+    scale_reference = pars['scale_reference']
+
+    output_model = apply_emicorr(input_model, emicorr_model, save_onthefly_reffile,
+                        save_intermediate_results=save_intermediate_results,
+                        user_supplied_reffile=user_supplied_reffile,
+                        nints_to_phase=nints_to_phase,
+                        nbins=nbins,
+                        scale_reference=scale_reference
+                        )
+
+    return output_model
+
+
+def apply_emicorr(input_model, emicorr_model, save_onthefly_reffile,
+        save_intermediate_results=False, user_supplied_reffile=None,
+        nints_to_phase=None, nbins=None, scale_reference=True):
+    """
+    -> NOTE: This is translated from IDL code fix_miri_emi.pro
+
+    EMI-corrects data by the following procedure:
+         [repeat recursively for each discrete EMI frequency desired]
+    1) Read image data.
+    2) Make very crude slope image and fixed pattern "super"bias for each
+       integration, ignoring everything (nonlin, saturation, badpix, etc).
+    3) Subtract scaled slope image and bias from each frame of each integration.
+    4) Calculate phase of every pixel in the image at the desired emi frequency
+       (e.g. 390 Hz) relative to the first pixel in the image.
+    5) Make a binned, phased amplitude (pa) wave from the cleaned data (plot
+       cleaned vs phase, then bin by phase).
+    6)* Measure the phase shift between the binned pa wave and the input
+       reference wave
+    7)* Use look-up table to get the aligned reference wave value for each pixel
+       (the noise array corresponding to the input image).
+
+     6a-7a) Alternately, use the binned pa wave instead of the reference wave to
+       "self-correct"
+
+    8) Subtract the noise array from the input image and return the cleaned result.
+
+         [repeat for next frequency using cleaned output as input]
+      [removing highest amplitude emi first will probably give best results]
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodels.JwstDataModel`
+        input science data model to be emi-corrected
+
+    emicorr_model : `~jwst.datamodels.EmiModel`
+         ImageModel of emi
+
+    save_onthefly_reffile : str or None
+        Full path and root name of on-the-fly reference file
+
+    save_intermediate_results : bool
+        Saves the output into a file and the reference file (if created on-the-fly)
+
+    user_supplied_reffile : str
+        Reference file supplied by the user
+
+    nints_to_phase : int
+        Number of integrations to phase
+
+    nbins : int
+        Number of bins in one phased wave
+
+    scale_reference : bool
+        If True, the reference wavelength will be scaled to the data's phase amplitude
+
+    Returns
+    -------
+    output_model : JWST data model
+        input science data model which has been emi-corrected
+    """
+    # get the subarray case and other info
+    detector = input_model.meta.instrument.detector
+    subarray = input_model.meta.subarray.name
+    readpatt = input_model.meta.exposure.readpatt
+    xsize = input_model.meta.subarray.xsize   # SUBSIZE1 keyword
+    xstart = input_model.meta.subarray.xstart   # SUBSTRT1 keyword
+
+    # get the subarray case from either the ref file or set default values
+    freqs_numbers = []
+    if emicorr_model is not None:
+        log.info('Using reference file to get subarray case.')
+        subname, rowclocks, frameclocks, freqs2correct = get_subarcase(emicorr_model, subarray, readpatt, detector)
+        reference_wave_list = []
+        for fnme in freqs2correct:
+            freq, ref_wave = get_frequency_info(emicorr_model, fnme)
+            freqs_numbers.append(freq)
+            reference_wave_list.append(ref_wave)
+    else:
+        log.info('Using default subarray case corrections.')
+        subname, rowclocks, frameclocks, freqs2correct = get_subarcase(default_subarray_cases, subarray, readpatt, detector)
+        for fnme in freqs2correct:
+            freq = get_frequency_info(default_emi_freqs, fnme)
+            freqs_numbers.append(freq)
+
+    if rowclocks is None:
+        # no subarray match found, print to log and skip correction
+        return subname
+
+    # get the number of samples, 10us sample times per pixel (1 for fastmode, 9 for slowmode)
+    nsamples = input_model.meta.exposure.nsamples
+
+    # Initialize the output model as a copy of the input
+    output_model = input_model.copy()
+
+    # create the dictionary to store the frequencies and corresponding phase amplitudes
+    if save_intermediate_results and save_onthefly_reffile is not None:
+        freq_pa_dict = {'frequencies': {}, 'subarray_cases': {}}
+
+    # Loop over the frequencies to correct
+    log.info('With configuration: Subarray={}, Read_pattern={}, Detector={}'.format(subarray, readpatt, detector))
+    log.info('Will correct data for the following {} frequencies: '.format(len(freqs2correct)))
+    log.info('   {}'.format(freqs2correct))
+    for fi, frequency_name in enumerate(freqs2correct):
+        frequency = freqs_numbers[fi]
+        log.info('Correcting for frequency: {} Hz  ({} out of {})'.format(frequency, fi+1, len(freqs2correct)))
+
+        # Read image data and set up some variables
+        orig_data = output_model.data
+        data = orig_data.copy()
+        nints, ngroups, ny, nx = np.shape(data)
+        if nints_to_phase is None:
+            nints_to_phase = nints
+
+        # Correspondance of array order in IDL
+        # sz[0] = 4 in idl
+        # sz[1] = nx
+        # sz[2] = ny
+        # sz[3] = ngroups
+        # sz[4] = nints
+        nx4 = int(nx/4)
+
+        dd_all = np.ones((nints, ngroups, ny, nx4))
+        log.info('Subtracting self-superbias from each group of each integration')
+        for ninti in range(nints_to_phase):
+            log.debug('  Working on integration: {}'.format(ninti+1))
+
+            # Remove source signal and fixed bias from each integration ramp
+            # (linear is good enough for phase finding)
+
+            # do linear fit for source + sky
+            s0, mm0 = sloper(data[ninti, 1:ngroups-1, :, :])
+
+            # subtract source+sky from each frame of this ramp
+            for ngroupi in range(ngroups):
+                data[ninti, ngroupi, ...] = orig_data[ninti, ngroupi, ...] - (s0 * ngroupi)
+
+            # make a self-superbias
+            m0 = minmed(data[ninti, 1:ngroups-1, :, :])
+
+            # subtract self-superbias from each frame of this ramp
+            for ngroupi in range(ngroups):
+                data[ninti, ngroupi, ...] = data[ninti, ngroupi, ...] - m0
+
+                # de-interleave each frame into the 4 separate output channels and
+                # average (or median) them together for S/N
+                d0 = data[ninti, ngroupi, :, 0:nx:4]
+                d1 = data[ninti, ngroupi, :, 1:nx:4]
+                d2 = data[ninti, ngroupi, :, 2:nx:4]
+                d3 = data[ninti, ngroupi, :, 3:nx:4]
+                dd = (d0 + d1 + d2 + d3)/4.
+
+                # fix a bad ref col
+                dd[:, 1] = (dd[:, 0] + dd[:, 3])/2
+                dd[:, 2] = (dd[:, 0] + dd[:, 3])/2
+                # This is the quad-averaged, cleaned, input image data for the exposure
+                dd_all[ninti, ngroupi, ...] = dd - np.median(dd)
+
+        # Calculate times of all pixels in the input integration, then use that to calculate
+        # phase of all pixels. Times here is in integer numbers of 10us pixels starting from
+        # the first data pixel in the input image. Times can be a very large integer, so use
+        # a big datatype. Phaseall (below) is just 0-1.0.
+
+        # A safer option is to calculate times_per_integration and calculate the phase at each
+        # int separately. That way times array will have a smaller number of elements at each
+        # int, with less risk of datatype overflow. Still, use the largest datatype available
+        # for the time_this_int array.
+
+        times_this_int = np.zeros((ngroups, ny, nx4), dtype='ulonglong')
+        phaseall = np.zeros((nints, ngroups, ny, nx4))
+
+        # non-roi rowclocks between subarray frames (this will be 0 for fullframe)
+        extra_rowclocks = (1024. - ny) * (4 + 3.)
+        # e.g. ((1./390.625) / 10e-6) = 256.0 pix and ((1./218.52055) / 10e-6) = 457.62287 pix
+        period_in_pixels = (1./frequency) / 10.0e-6
+
+        start_time, ref_pix_sample = 0, 3
+
+        # Need colstop for phase calculation in case of last refpixel in a row. Technically,
+        # this number comes from the subarray definition (see subarray_cases dict above), but
+        # calculate it from the input image header here just in case the subarray definitions
+        # are not available to this routine.
+        colstop = int( xsize/4 + xstart - 1 )
+        log.info('Phase calculation per integration')
+        for l in range(nints_to_phase):
+            log.debug('  Working on integration: {}'.format(l+1))
+            for k in range(ngroups):   # frames
+                for j in range(ny):    # rows
+                    # nsamples= 1 for fast, 9 for slow (from metadata)
+                    times_this_int[k, j, :] = np.arange(nx4, dtype='ulonglong') * nsamples + start_time
+
+                    # If the last pixel in a row is a reference pixel, need to push it out
+                    # by ref_pix_sample sample times. The same thing happens for the first
+                    # ref pix in each row, but that gets absorbed into the inter-row pad and
+                    # can be ignored here. Since none of the current subarrays hit the
+                    # right-hand reference pixel, this correction is not in play, but for
+                    # fast and slow fullframe (e.g. 10Hz) it should be applied. And even
+                    # then, leaving this out adds just a *tiny* phase error on the last ref
+                    # pix in a row (only) - it does not affect the phase of the other pixels.
+
+                    if colstop == 258:
+                        ulonglong_ref_pix_sample = ref_pix_sample + 2**32
+                        times_this_int[k, j, nx4-1] = times_this_int[k, j, nx4-1] + ulonglong_ref_pix_sample
+
+                    # point to the first pixel of the next row (add "end-of-row" pad)
+                    start_time += rowclocks
+
+                # point to the first pixel of the next frame (add "end-of-frame" pad)
+                start_time += extra_rowclocks
+
+            # Convert "times" to phase each integration. Nothe times has units of
+            # number of 10us from the first data pixel in this integration, so to
+            # convert to phase, divide by the waveform *period* in float pixels
+            phase_this_int = times_this_int / period_in_pixels
+            phaseall[l, ...] = phase_this_int - phase_this_int.astype('ulonglong')
+
+            # add a frame time to account for the extra frame reset between MIRI integrations
+            start_time += frameclocks
+
+        # use phaseall vs dd_all
+
+        # Define the sizew of 1 wave of the phased waveform vector, then bin the whole
+        # dataset at this interval. This is essentially the desired number of bins along
+        # the waveform. This can be any number that is at least 1 less than the period in
+        # pixels. If larger, some bins could end up sparsely sampled. Fewer bins results
+        # in a smoother waveform but lower resolution. Can always be smoothed and/or
+        # streched to a different dimension later. By default, use period_in_pixels/2.0
+
+        # number of bins in one phased wave (e.g. 255, 218, etc)
+        if nbins is None:
+            # the IDL code sets nbins as ulong type (ulonglong in python)
+            nbins = int(period_in_pixels/2.0)
+
+        # bin the whole set
+        log.info('Calculating the phase amplitude for {} bins'.format(nbins))
+        # Define the binned waveform amplitude (pa = phase amplitude)
+        pa = np.arange(nbins, dtype=float)
+        # keep track of n per bin to check for low n
+        nu = np.arange(nbins)
+        for nb in range(nbins):
+            u = np.where((phaseall > nb/nbins) & (phaseall <= (nb + 1)/nbins) & (np.isfinite(dd_all)))
+            nu[nb] = phaseall[u].size
+            # calculate the sigma-clipped mean
+            dmean, _, _, _ = iter_stat_sig_clip(dd_all[u])
+            pa[nb] = dmean   # amplitude in this bin
+
+        pa -= np.median(pa)
+
+        # pa_phase is the phase corresponding to each bin in pa. The +0.5 here is to center in
+        # the bin. Together, (pa, pa_phase) define the waveform to phase against the
+        # reference_wave, or to use directly for a self correction (likely noiser than the
+        # reference_wave, but perfectly matched in phase by definition so no need for a
+        # crosscor or other phase-finder against a reference_wave). This may not be
+        # necessary. Can just assume pa and the reference_wave are exactly 1.0 waves long.
+        #pa_phase = (np.arange(nbins, dtype=np.float32) + 0.5) / nbins
+
+        # Correction using direct look-up from either pa (self-correction) or
+        # the reference wave (if provided)
+        # Make a vector to hold a version of pa that is close to an integer number of
+        # pixels long. This will be the noise amplitude look-up table.
+        lut = rebin(pa, [period_in_pixels])   # shift and resample reference_wave at pa's phase
+
+        # If a reference_wave_filename was provided, use it to measure phase shift between
+        # pa (binned phased wave) from the input image vs. the reference_wave using xcorr.
+        # These two methods give the same integer-reference_wave-element resolution results.
+        if emicorr_model is not None:
+            log.info('Using reference file to measure phase shift')
+            reference_wave = reference_wave_list[fi]
+            reference_wave_size = np.size(reference_wave)
+            rebinned_pa = rebin(pa, [reference_wave_size])
+            cc = np.zeros(reference_wave_size)
+            for i in range(reference_wave_size):
+                shifted_ref_wave = np.roll(reference_wave, i)
+                pears_coeff = np.corrcoef(shifted_ref_wave, rebinned_pa)
+                cc[i] = pears_coeff[0, 1]
+
+            # Use the reference file (shifted to match the phase of pa,
+            # and optionally amplitude scaled)
+            # shift and resample reference_wave at pa's phase
+            # u[0] is the phase shift of reference_wave *to* pa
+            u = np.where(cc >= max(cc))
+            lut_reference = rebin(np.roll(reference_wave, u[0]), [period_in_pixels])
+
+            # Scale reference wave amplitude to match the pa amplitude from this dataset by
+            # fitting a line rather than taking a mean ratio so that any DC offset drops out
+            # in the intercept (and to avoid potential divide-by-zeros)
+            m, b = np.polyfit(lut_reference, lut, 1)   # the slope is the amplitude ratio
+
+            # check pa-to-reference_wave phase match and optional scale
+            if scale_reference:
+                lut_reference = lut_reference * m
+
+            lut = lut_reference
+
+        if save_intermediate_results and save_onthefly_reffile is not None:
+            freq_pa_dict['frequencies'][frequency_name] = {'frequency' : frequency,
+                                                        'phase_amplitudes' : pa}
+
+        log.info('Creating phased-matched noise model to subtract from data')
+        # This is the phase matched noise model to subtract from each pixel of the input image
+        dd_noise = lut[(phaseall * period_in_pixels).astype(int)]
+
+        # Interleave (straight copy) into 4 amps
+        noise = np.ones((nints, ngroups, ny, nx))   # same size as input data
+        for k in range(4):
+            noise_x = np.arange(nx4)*4 + k
+            noise[..., noise_x] = dd_noise
+
+        # Subtract EMI noise from the input data
+        log.info('Subtracting EMI noise from data')
+        corr_data = orig_data - noise
+        output_model.data = corr_data
+
+    if save_intermediate_results and save_onthefly_reffile is not None:
+        if readpatt == 'FASTR1':
+            freqs_dict = {readpatt: freqs2correct}
+        else:
+            freqs_dict = {readpatt: {detector: freqs2correct} }
+        on_the_fly_subarr_case = {}
+        on_the_fly_subarr_case[subarray] = {
+            'rowclocks': rowclocks,
+            'frameclocks': frameclocks,
+            'freqs': freqs_dict }
+        freq_pa_dict['subarray_cases'] = on_the_fly_subarr_case
+        mk_reffile(freq_pa_dict, save_onthefly_reffile)
+
+    return output_model
+
+
+def sloper(data):
+    """ Fit slopes to all pix of a ramp, using numerical recipies plane-adding
+     returning intercept image.
+
+    Parameters
+    ----------
+    data : numpy array
+        3-D integration data array
+
+    Returns
+    -------
+    outarray : numpy array
+        3-D integration data array
+
+    intercept: numpy array
+        slope intercept values
+    """
+    ngroups, ny, nx = np.shape(data)
+    frametime = 1.0   # use 1.0 for per-frame slopes, otherwise put a real time here
+    grouptimes = np.arange(ngroups) * frametime
+
+    sxy = np.zeros((ny, nx))
+    sx = np.zeros((ny, nx))
+    sy = np.zeros((ny, nx))
+    sxx = np.zeros((ny, nx))
+
+    for groupcount in range(ngroups):
+        # do the incremental slope calculation operations
+        sxy = sxy + grouptimes[groupcount] * data[groupcount, ...]
+        sx = sx + grouptimes[groupcount]
+        sy = sy + data[groupcount, ...]
+        sxx = sxx + grouptimes[groupcount]**2
+
+    # calculate the final per-pixel slope values
+    outarray = (ngroups * sxy - sx * sy)/(ngroups * sxx - sx * sx)
+
+    # calculate the final per-pixel intercept values
+    intercept = (sxx * sy - sx * sxy)/(ngroups * sxx - sx * sx)
+
+    return outarray, intercept
+
+
+def minmed(data, minval=False, avgval=False, maxval=False):
+    """ Returns the median image of a stack of images, or if there are 2 or less
+    non-zero frames, returns the minimum (or this minimum can be forced).
+
+    Parameters
+    ----------
+    data : numpy array
+        3-D integration data array
+
+    minval : bool
+        Returned image will be the minimum
+
+    avgval : bool
+        Returned image will be the mean
+
+    maxval : bool
+        Returned image will be the maximum
+
+    Returns
+    -------
+    medimg : numpy array
+        Median image of a stack of images
+    """
+
+    ngroups, ny, nx = np.shape(data)
+    medimg = np.zeros((ny, nx))
+    # use a mask to ignore nans for calculations
+    masked_data = np.ma.array(data, mask=np.isnan(data))
+
+    for i in range(nx):
+        for j in range(ny):
+            vec = masked_data[:, j, i]
+            u = np.where(vec != 0)
+            n = vec[u].size
+            if n > 0:
+                if n <= 2 or minval:
+                    medimg[j, i] = np.ma.min(vec[u])
+                if maxval:
+                    medimg[j, i] = np.ma.max(vec[u])
+                if not minval and not maxval and not avgval:
+                    medimg[j, i] = np.ma.median(vec[u])
+                if avgval:
+                    dmean , _, _, _ = iter_stat_sig_clip(vec[u])
+                    medimg[j, i] = dmean
+    return medimg
+
+
+def get_subarcase(subarray_cases, subarray, readpatt, detector):
+    """ Get the rowclocks and frameclocks values for the given configuration.
+
+    Parameters
+    ----------
+    subarray_cases : dict or model object
+        Either default corrections dictionary or datamodel
+
+    subarray : str
+        Keyword value
+
+    readpatt : str
+        Keyword value
+
+    detector : str
+        Keyword value
+
+    Returns
+    -------
+    subname : str
+        Modified subarray name
+
+    rowclocks : int
+
+    frameclocks : int
+
+    frequencies : list
+        List of frequencies to correct according to subarray name
+
+    freqs_names : list
+        List of frequency names to use
+    """
+    subname, rowclocks, frameclocks, frequencies = None, None, None, None
+    if isinstance(subarray_cases, dict):
+        for subname in subarray_cases:
+            if subarray not in subname:
+                continue
+            if subname == 'FULL':
+                subname = subname + '_' + readpatt
+            rowclocks = subarray_cases[subname]["rowclocks"]
+            frameclocks = subarray_cases[subname]["frameclocks"]
+            if readpatt == "SLOWR1":
+                frequencies = subarray_cases[subname]["freqs"]["SLOWR1"][detector]
+            else:
+                frequencies = subarray_cases[subname]["freqs"]["FASTR1"]
+            break
+    else:
+        frequencies = []
+        mdl_dict = subarray_cases.to_flat_dict()
+        for subname in subarray_cases.subarray_cases:
+            subname = subname.split(sep='.')[0]
+            if subarray not in subname:
+                continue
+            log.debug('Found subarray case {}!'.format(subname))
+            if 'FULL' in subname:
+                subname = subname + '_' + readpatt
+            for item, val in mdl_dict.items():
+                if subname in item:
+                    if "rowclocks" in item:
+                        rowclocks = val
+                    elif "frameclocks" in item:
+                        frameclocks = val
+                    else:
+                        if readpatt == "SLOWR1" and "SLOWR1" in item and detector in item:
+                            frequencies.append(val)
+                        elif readpatt == "FASTR1" and "FASTR1" in item:
+                            frequencies.append(val)
+            if subname is not None and rowclocks is not None and frameclocks is not None and frequencies is not None:
+                break
+    return subname, rowclocks, frameclocks, frequencies
+
+
+def get_frequency_info(freqs_names_vals, frequency_name):
+    """Get the frequency number from the given dictionary
+    Parameters
+    ----------
+    freqs_names_vals : dict or model object
+        Either default corrections dictionary or datamodel
+
+    frequency_name : str
+        Frequency of interest
+
+    Returns
+    -------
+    frequency_number : float
+        Frequency
+
+    phase_amplitudes : array
+        1-D array of the corresponding phase amplidues for this frequency
+    """
+    if isinstance(freqs_names_vals, dict):
+        for freq_nme, val in freqs_names_vals.items():
+            if freq_nme == frequency_name:
+                return val
+    else:
+        freq_number, phase_amplitudes = None, None
+        mdl_dict = freqs_names_vals.to_flat_dict()
+        for item, val in mdl_dict.items():
+            if frequency_name in item:
+                if 'frequency' in item:
+                    freq_number = val
+                if 'phase_amplitudes' in item:
+                    phase_amplitudes = val
+            if freq_number is not None and phase_amplitudes is not None:
+                break
+        return freq_number, phase_amplitudes
+
+
+def iter_stat_sig_clip(data, sigrej=3.0, maxiter=10):
+    """ Compute the mean, mediand and/or sigma of data with iterative sigma clipping.
+    This funtion is based on djs_iterstat.pro (authors therein)
+
+    Parameters
+    ----------
+    data : numpy array
+
+    sigrej : float
+        Sigma for rejection
+
+    maxiter: int
+        Maximum number of sigma rejection iterations
+
+    Returns
+    -------
+    dmean : float
+        Computed mean
+
+    dmedian : float
+        Computed median
+
+    dsigma : float
+        Computed sigma
+
+    dmask : numpy array
+        Mask set to 1 for good points, and 0 for rejected points
+    """
+    # special cases of 0 or 1 data points
+    ngood = np.size(data)
+    dmean, dmedian, dsigma, dmask = 0.0, 0.0, 0.0, np.zeros(np.shape(data))
+    if ngood == 0:
+        log.debug('No data points for sigma clipping')
+        return dmean, dmedian, dsigma, dmask
+    elif ngood == 1:
+        log.debug('Only 1 data point for sigma clipping')
+        return dmean, dmedian, dsigma, dmask
+
+    # Compute the mean + standard deviation of the entire data array,
+    # these values will be returned if there are fewer than 2 good points.
+    dmask = np.ones(ngood, dtype='b') + 1
+    dmean = sum(data * dmask) / ngood
+    dsigma = np.sqrt(sum((data - dmean)**2) / (ngood - 1))
+    dsigma = dsigma
+    iiter = 1
+
+    # Iteratively compute the mean + stdev, updating the sigma-rejection thresholds
+    # each iteration
+    nlast = -1
+    while iiter < maxiter and nlast != ngood and ngood >= 2:
+        loval = dmean - sigrej * dsigma
+        hival = dmean + sigrej * dsigma
+        nlast = ngood
+
+        dmask[data < loval] = 0
+        dmask[data > hival] = 0
+        ngood = sum(dmask)
+
+        if ngood >= 2:
+            dmean = sum(data*dmask) / ngood
+            dsigma = np.sqrt( sum((data - dmean)**2 * dmask) / (ngood - 1) )
+            dsigma = dsigma
+
+        iiter += 1
+
+    return dmean, dmedian, dsigma, dmask
+
+
+def rebin(arr, newshape):
+    """Rebin an array to a new shape.
+
+    Parameters
+    ----------
+    arr : numpy array
+        array to rebin
+
+    newshape : list
+        New shape for the input array
+
+    Returns
+    -------
+    arr : numpy array
+        rebinned array
+    """
+    assert len(arr.shape) == len(newshape)
+
+    slices = [ slice(0,old, float(old)/new) for old,new in zip(arr.shape,newshape) ]
+    coordinates = np.mgrid[slices]
+    indices = coordinates.astype('i')   # choose the biggest smaller integer index
+
+    return arr[tuple(indices)]
+
+
+def mk_reffile(freq_pa_dict, emicorr_ref_filename):
+    """ Create the reference file hdulist object.
+
+    Parameters
+    ----------
+    freq_pa_dict : dictionary
+        Dictionary containing the phase amplitude (pa) array
+        corresponding to the appropriate frequency
+
+    emicorr_ref_filename : str
+        Full path and root name of on-the-fly reference file
+
+    Returns
+    -------
+    Nothing
+    """
+    # save the reference file if save_intermediate_results=True and no CRDS file exists
+    emicorr_model = datamodels.EmiModel(freq_pa_dict)
+    if 'fits' in emicorr_ref_filename:
+        emicorr_ref_filename = emicorr_ref_filename.replace('fits', 'asdf')
+    emicorr_model.save(emicorr_ref_filename)
+    emicorr_model.close()
+    log.info('On-the-fly reference file written as: %s', emicorr_ref_filename)
+
+
diff --git a/jwst/emicorr/emicorr_step.py b/jwst/emicorr/emicorr_step.py
new file mode 100755
index 000000000..c9a6e90bd
--- /dev/null
+++ b/jwst/emicorr/emicorr_step.py
@@ -0,0 +1,86 @@
+#! /usr/bin/env python
+
+
+from stdatamodels.jwst import datamodels
+from ..stpipe import Step
+from . import emicorr
+
+
+__all__ = ["EmiCorrStep"]
+
+
+class EmiCorrStep(Step):
+    """
+    EmiCorrStep: Apply MIRI EMI correction to a science image.
+    """
+
+    class_alias = "emicorr"
+
+    spec = """
+        save_intermediate_results = boolean(default=False)
+        user_supplied_reffile = string(default=None)  # ASDF user-supplied reference file
+        nints_to_phase = integer(default=None)  # Number of integrations to phase
+        nbins = integer(default=None)  # Number of bins in one phased wave
+        scale_reference = boolean(default=True)  # If True, the reference wavelength will be scaled to the data's phase amplitude
+    """
+
+    reference_file_types = ['emicorr']
+
+    def process(self, input):
+        with datamodels.open(input) as input_model:
+
+            # Catch the cases to skip, i.e. all instruments other than MIRI
+            instrument = input_model.meta.instrument.name
+            if instrument != 'MIRI':
+                self.log.warning('EMI correction not implemented for instrument: {}'.format(instrument))
+                input_model.meta.cal_step.emicorr = 'SKIPPED'
+                return input_model
+
+            # Setup parameters
+            pars = {
+                'save_intermediate_results': self.save_intermediate_results,
+                'user_supplied_reffile': self.user_supplied_reffile,
+                'nints_to_phase': self.nints_to_phase,
+                'nbins': self.nbins,
+                'scale_reference': self.scale_reference
+            }
+
+            # Get the reference file
+            save_onthefly_reffile, emicorr_ref_filename, emicorr_model = None, None, None
+            if self.user_supplied_reffile is None:
+                try:
+                    emicorr_ref_filename = self.get_reference_file(input_model, 'emicorr')
+                    # Create the reference file only if told to save outputs, else correct on-the-fly
+                    if emicorr_ref_filename == 'N/A':
+                        emicorr_ref_filename = None
+                    else:
+                        self.log.info('Using CRDS reference file: {}'.format(emicorr_ref_filename))
+                        emicorr_model = datamodels.EmiModel(emicorr_ref_filename)
+                except Exception:
+                    # No reference file in CRDS
+                    self.log.info('No CRDS emicorr reference file found. Creating on-the-fly reference file.')
+
+            else:
+                self.log.info('Using user-supplied reference file: {}'.format(self.user_supplied_reffile))
+                emicorr_model = datamodels.EmiModel(self.user_supplied_reffile)
+
+            # Do the correction
+            if self.save_intermediate_results:
+                if emicorr_ref_filename is None and self.user_supplied_reffile is None:
+                    # get the same full path as input file to save the on-the-fly reference file
+                    emicorr_ref_filename = Step._make_output_path(self, suffix='emi_ref_waves')
+                    save_onthefly_reffile = emicorr_ref_filename
+                else:
+                    save_onthefly_reffile = None
+            output_model = emicorr.do_correction(input_model, emicorr_model, save_onthefly_reffile, **pars)
+            if isinstance(output_model, str) or output_model is None:
+                # in this case output_model=subarray_readpatt configuration
+                self.log.warning('No correction match for this configuration: {}'.format(output_model))
+                self.log.warning('Step skipped')
+                input_model.meta.cal_step.emicorr = 'SKIPPED'
+                return input_model
+
+            # close and remove the reference file created on-the-fly
+            output_model.meta.cal_step.emicorr = 'COMPLETE'
+
+        return output_model
diff --git a/jwst/extract_2d/nirspec.py b/jwst/extract_2d/nirspec.py
index 571ad4256..aeccaea54 100644
--- a/jwst/extract_2d/nirspec.py
+++ b/jwst/extract_2d/nirspec.py
@@ -46,9 +46,12 @@ def nrs_extract2d(input_model, slit_name=None):
     open_slits = slit2msa.slits[:]
     if slit_name is not None:
         new_open_slits = []
-        slit_name = str(slit_name)
+        if not isinstance(slit_name, str):
+            slit_name = str(slit_name)
         for sub in open_slits:
-            if str(sub.name) == slit_name:
+            if not isinstance(sub.name, str):
+                str_subnme = str(sub.name)
+            if str_subnme == slit_name:
                 new_open_slits.append(sub)
                 break
         if len(new_open_slits) == 0:
diff --git a/jwst/lib/suffix.py b/jwst/lib/suffix.py
index 6881a221c..7f85215ec 100644
--- a/jwst/lib/suffix.py
+++ b/jwst/lib/suffix.py
@@ -190,6 +190,7 @@ _calculated_suffixes = {
  'persistence',
  'image2pipeline',
  'klip'
+ 'emicorrstep',
 }
 
 
diff --git a/jwst/pipeline/calwebb_dark.py b/jwst/pipeline/calwebb_dark.py
index 6d37dff72..9673ffd88 100644
--- a/jwst/pipeline/calwebb_dark.py
+++ b/jwst/pipeline/calwebb_dark.py
@@ -8,6 +8,7 @@ from ..stpipe import Pipeline
 # step imports
 from ..group_scale import group_scale_step
 from ..dq_init import dq_init_step
+from ..emicorr import emicorr_step
 from ..saturation import saturation_step
 from ..ipc import ipc_step
 from ..superbias import superbias_step
@@ -40,6 +41,7 @@ class DarkPipeline(Pipeline):
     # Define aliases to steps
     step_defs = {'group_scale': group_scale_step.GroupScaleStep,
                  'dq_init': dq_init_step.DQInitStep,
+                 'emicorr': emicorr_step.EmiCorrStep,
                  'saturation': saturation_step.SaturationStep,
                  'ipc': ipc_step.IPCStep,
                  'superbias': superbias_step.SuperBiasStep,
@@ -67,6 +69,7 @@ class DarkPipeline(Pipeline):
 
             input = self.group_scale(input)
             input = self.dq_init(input)
+            input = self.emicorr(input)
             input = self.saturation(input)
             input = self.ipc(input)
             input = self.firstframe(input)
diff --git a/jwst/pipeline/calwebb_detector1.py b/jwst/pipeline/calwebb_detector1.py
index ed85ba213..16087ee73 100644
--- a/jwst/pipeline/calwebb_detector1.py
+++ b/jwst/pipeline/calwebb_detector1.py
@@ -8,6 +8,7 @@ from ..stpipe import Pipeline
 # step imports
 from ..group_scale import group_scale_step
 from ..dq_init import dq_init_step
+from ..emicorr import emicorr_step
 from ..saturation import saturation_step
 from ..ipc import ipc_step
 from ..superbias import superbias_step
@@ -49,6 +50,7 @@ class Detector1Pipeline(Pipeline):
     # Define aliases to steps
     step_defs = {'group_scale': group_scale_step.GroupScaleStep,
                  'dq_init': dq_init_step.DQInitStep,
+                 'emicorr': emicorr_step.EmiCorrStep,
                  'saturation': saturation_step.SaturationStep,
                  'ipc': ipc_step.IPCStep,
                  'superbias': superbias_step.SuperBiasStep,
@@ -87,6 +89,7 @@ class Detector1Pipeline(Pipeline):
 
             input = self.group_scale(input)
             input = self.dq_init(input)
+            input = self.emicorr(input)
             input = self.saturation(input)
             input = self.ipc(input)
             input = self.firstframe(input)
diff --git a/jwst/pipeline/emicorr.cfg b/jwst/pipeline/emicorr.cfg
new file mode 100644
index 000000000..6f24969c9
--- /dev/null
+++ b/jwst/pipeline/emicorr.cfg
@@ -0,0 +1,2 @@
+name = "emicorr"
+class = "jwst.emicorr.EmiCorrStep"


Test patch:
diff --git a/jwst/emicorr/tests/__init__.py b/jwst/emicorr/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/jwst/emicorr/tests/test_emicorr.py b/jwst/emicorr/tests/test_emicorr.py
new file mode 100644
index 000000000..80e8cb01c
--- /dev/null
+++ b/jwst/emicorr/tests/test_emicorr.py
@@ -0,0 +1,144 @@
+"""
+
+Unit tests for EMI correction
+
+"""
+
+
+import numpy as np
+from jwst.emicorr import emicorr, emicorr_step
+from stdatamodels.jwst.datamodels import Level1bModel, EmiModel
+
+default_subarray_cases = emicorr.default_subarray_cases
+
+def mk_data_mdl(data, subarray, readpatt, detector):
+    # create input_model
+    input_model = Level1bModel(data=data)
+    input_model.meta.instrument.name = 'MIRI'
+    input_model.meta.instrument.detector = detector
+    input_model.meta.exposure.type = 'MIR_4QPM'
+    input_model.meta.subarray.name = subarray
+    input_model.meta.exposure.readpatt = readpatt
+    input_model.meta.subarray.xsize = 288
+    input_model.meta.subarray.xstart = 1
+    input_model.meta.exposure.nsamples = 1
+    return input_model
+
+
+def test_EmiCorrStep():
+    data = np.ones((1, 5, 20, 20))
+    input_model = mk_data_mdl(data, 'MASK1550', 'FASTR1', 'MIRIMAGE')
+    expected_outmdl = input_model.copy()
+
+    nirmdl = input_model.copy()
+    nirmdl.meta.instrument.name = 'NIRISS'
+
+    stp = emicorr_step.EmiCorrStep()
+    nir_result = stp.call(nirmdl)
+
+    # expect no change because instrument is not MIRI
+    assert expected_outmdl.data.all() == nir_result.data.all()
+
+
+def test_do_correction():
+    data = np.ones((1, 5, 20, 20))
+    input_model = mk_data_mdl(data, 'MASK1550', 'FASTR1', 'MIRIMAGE')
+    pars = {
+        'save_intermediate_results': False,
+        'user_supplied_reffile': None,
+        'nints_to_phase': None,
+        'nbins': None,
+        'scale_reference': True
+    }
+    emimdl = {'frequencies': {
+                "Hz390": {'frequency': 390.625,
+                    'phase_amplitudes': np.ones(20)+0.1},
+                "Hz10": { 'frequency': 10.039216,
+                    'phase_amplitudes': np.ones(20)+0.5}},
+            'subarray_cases': default_subarray_cases
+    }
+    save_onthefly_reffile = None
+    emicorr_model = EmiModel(emimdl)
+    outmdl = emicorr.do_correction(input_model, emicorr_model, save_onthefly_reffile, **pars)
+
+    assert outmdl is not None
+
+
+def test_apply_emicorr():
+    data = np.ones((1, 5, 20, 20))
+    input_model = mk_data_mdl(data, 'MASK1550', 'FASTR1', 'MIRIMAGE')
+    emicorr_model, save_onthefly_reffile = None, None
+    outmdl = emicorr.apply_emicorr(input_model, emicorr_model, save_onthefly_reffile,
+                save_intermediate_results=False, user_supplied_reffile=None,
+                nints_to_phase=None, nbins=None, scale_reference=True)
+
+    assert outmdl is not None
+
+
+def test_get_subarcase():
+    subarray, readpatt, detector = 'MASK1550', 'FASTR1', 'MIRIMAGE'
+    subarray_info_r = emicorr.get_subarcase(default_subarray_cases, subarray, readpatt, detector)
+    subname_r, rowclocks_r, frameclocks_r, freqs2correct_r = subarray_info_r
+
+    # set up a fake configuration
+    subarray_info_f = emicorr.get_subarcase(default_subarray_cases, 'FULL', 'kkkkk', detector)
+    subname_f, rowclocks_f, frameclocks_f, freqs2correct_f = subarray_info_f
+
+    # test if we get the right configuration
+    compare_real = ['MASK1550', 82, 23968, ["Hz390", "Hz10"]]
+    subname_real, rowclocks_real, frameclocks_real, freqs2correct_real = compare_real
+
+    assert subname_real == subname_r
+    assert rowclocks_real == rowclocks_r
+    assert frameclocks_real == frameclocks_r
+    assert freqs2correct_real[0] == freqs2correct_r[0]
+
+
+def test_sloper():
+    data = np.ones((5, 5, 5))
+    outarray, intercept = emicorr.sloper(data)
+
+    compare_arr, compare_intercept = np.zeros((5, 5)), np.ones((5, 5))
+
+    assert compare_arr.all() == outarray.all()
+    assert compare_intercept.all() == intercept.all()
+
+
+def test_minmed():
+    data = np.ones((5, 5, 5))
+    compare_arr = data.copy()
+
+    medimg = emicorr.minmed(data)
+
+    assert compare_arr.all() == medimg.all()
+
+
+def test_iter_stat_sig_clip():
+    data = np.ones((5, 5, 5, 5))
+    data[0, 0, 0, 0] = 0.55
+    data[1, 1, 1, 1] = 1.11
+    data[2, 2, 2, 2] = 2.22
+    data[3, 3, 3, 3] = 3.33
+    data[4, 4, 4, 4] = 4.44
+    u = np.where(data > 1)
+    dmean, dmedian, dsigma, dmask = emicorr.iter_stat_sig_clip(data[u])
+
+    compare_mean, compare_median = 2.7750000000000004, 0.0
+    compare_sigma, compare_mask = 1.326703756361177, np.ones(np.size(data), dtype='b')+1
+
+    assert compare_mean == dmean
+    assert compare_median == dmedian
+    assert compare_sigma == dsigma
+    assert compare_mask.all() == dmask.all()
+
+
+def test_rebin():
+    data = np.ones(10)
+    data[1] = 0.55
+    data[5] = 1.55
+    data[9] = 2.0
+    rebinned_data = emicorr.rebin(data, [7])
+
+    compare_arr = np.array([1., 0.55, 1. , 1., 1.55, 1., 1.])
+
+    assert compare_arr.all() == rebinned_data.all()