From 0fff77f868d5efee1c5d0a142bbce1b8167020e9 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Sun, 25 Jan 2026 16:42:34 -0800
Subject: [PATCH 1/8] Remove weird latitudes from legacy EMTF z-files

---
 data/synthetic/emtf_results/from_matlab_256_26.zss | 2 +-
 data/synthetic/emtf_results/test2.zss              | 2 +-
 data/synthetic/emtf_results/test2r1.zrr            | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/data/synthetic/emtf_results/from_matlab_256_26.zss b/data/synthetic/emtf_results/from_matlab_256_26.zss
index e7e40edb..da98c7b0 100644
--- a/data/synthetic/emtf_results/from_matlab_256_26.zss
+++ b/data/synthetic/emtf_results/from_matlab_256_26.zss
@@ -2,7 +2,7 @@
  ********** WITH FULL ERROR COVARINCE**********
 Robust Single station                                                           
 station    :test1               
-coordinate  1007.996  0.0  declination  0.0
+coordinate  37.996  0.0  declination  0.0
 number of channels   5   number of frequencies   26
  orientations and tilts of each channel 
     1     0.00     0.00 tes  Hx
diff --git a/data/synthetic/emtf_results/test2.zss b/data/synthetic/emtf_results/test2.zss
index dd9e65da..6978cf02 100644
--- a/data/synthetic/emtf_results/test2.zss
+++ b/data/synthetic/emtf_results/test2.zss
@@ -2,7 +2,7 @@
  ********** WITH FULL ERROR COVARINCE**********
 Robust Single station                                                           
 station    :test2               
-coordinate  2951.999     0.000  declination     0.00
+coordinate  31.999     0.000  declination     0.00
 number of channels   5   number of frequencies   25
  orientations and tilts of each channel 
     1     0.00     0.00 tes  Hx
diff --git a/data/synthetic/emtf_results/test2r1.zrr b/data/synthetic/emtf_results/test2r1.zrr
index 1bed707e..a51781d7 100644
--- a/data/synthetic/emtf_results/test2r1.zrr
+++ b/data/synthetic/emtf_results/test2r1.zrr
@@ -2,7 +2,7 @@
  ********** WITH FULL ERROR COVARINCE**********
 Robust Remote Reference EMAP                                                    
 station    :test2r1             
-coordinate  1007.996   102.190  declination     0.00
+coordinate  37.996   102.190  declination     0.00
 number of channels   5   number of frequencies   25
  orientations and tilts of each channel 
     1     0.00     0.00 tes  Hx

From 1ce80c021475a1650c82c15e58f1410d03535fd8 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Sun, 25 Jan 2026 16:45:35 -0800
Subject: [PATCH 2/8] add mt_metadata z-file reader to test

- aslo kludge in old z-file reader to verify that old comparison plots still look correct
---
 aurora/sandbox/io_helpers/zfile_murphy.py     | 20 ++++++++++++++++++-
 ..._compare_aurora_vs_archived_emtf_pytest.py |  7 ++++++-
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/aurora/sandbox/io_helpers/zfile_murphy.py b/aurora/sandbox/io_helpers/zfile_murphy.py
index 8d397e9b..3d0ebdf6 100644
--- a/aurora/sandbox/io_helpers/zfile_murphy.py
+++ b/aurora/sandbox/io_helpers/zfile_murphy.py
@@ -243,7 +243,25 @@ def impedance(self, angle: Optional[float] = 0.0):
         u[hy_index, hy_index] = np.sin(
             (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
         )
-        u = np.linalg.inv(u)  # Identity if angle=0
+        try:
+            u = np.linalg.inv(u)
+        except np.linalg.LinAlgError:
+            u[hy_index, hx_index] = np.cos(
+                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
+            )
+        u[hy_index, hy_index] = np.sin(
+            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
+        )
+        try:
+            u = np.linalg.inv(u)
+        except np.linalg.LinAlgError:
+            u[hy_index, hx_index] = np.cos(
+                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
+            )
+            u[hy_index, hy_index] = np.sin(
+                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
+            )
+            u = np.linalg.inv(u)  # Identity if angle=0
 
         # build transformation matrix for predicted channels (electric fields)
         ex_index = self.channels.index("Ex")
diff --git a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
index 3eb85dc2..1604d614 100644
--- a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
+++ b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 from loguru import logger
 from mth5.helpers import close_open_files
@@ -13,6 +14,7 @@
     compute_rms,
     get_expected_rms_misfit,
 )
+from aurora.transfer_function.compare import CompareTF
 from aurora.transfer_function.emtf_z_file_helpers import (
     merge_tf_collection_to_match_z_file,
 )
@@ -69,7 +71,10 @@ def aurora_vs_emtf(
         z_file_path=z_file_path,
         return_collection=True,
     )
-
+    comparator = CompareTF(tf_01=z_file_path, tf_02=auxilliary_z_file)
+    result = comparator.compare_transfer_functions()
+    assert np.isclose(result["impedance_ratio"]["Z_10"], 1.0, atol=1e-2)
+    assert np.isclose(result["impedance_ratio"]["Z_01"], 1.0, atol=1e-2)
     aux_data = read_z_file(auxilliary_z_file)
     aurora_rho_phi = merge_tf_collection_to_match_z_file(aux_data, tf_collection)
     data_dict = {}

From f20385d9584ec34e350a9dcbe62addc442beae11 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Sun, 25 Jan 2026 21:11:55 -0800
Subject: [PATCH 3/8] remove ad-hoc z-file comparison from
 test_compare_aurora_vs_archived_emtf_pytest

- now use mt_metadata and aurora/transfer_function/compare.py
---
 aurora/sandbox/io_helpers/zfile_murphy.py     | 15 +---
 aurora/transfer_function/compare.py           |  5 +-
 ..._compare_aurora_vs_archived_emtf_pytest.py | 74 ++++++-------------
 3 files changed, 31 insertions(+), 63 deletions(-)

diff --git a/aurora/sandbox/io_helpers/zfile_murphy.py b/aurora/sandbox/io_helpers/zfile_murphy.py
index 3d0ebdf6..d4bd27a3 100644
--- a/aurora/sandbox/io_helpers/zfile_murphy.py
+++ b/aurora/sandbox/io_helpers/zfile_murphy.py
@@ -237,10 +237,10 @@ def impedance(self, angle: Optional[float] = 0.0):
         u[hx_index, hy_index] = np.sin(
             (self.orientation[hx_index, 0] - angle) * np.pi / 180.0
         )
-        u[hy_index, hx_index] = np.cos(
+        u[hy_index, hx_index] = np.sin(
             (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
         )
-        u[hy_index, hy_index] = np.sin(
+        u[hy_index, hy_index] = np.cos(
             (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
         )
         try:
@@ -252,16 +252,7 @@ def impedance(self, angle: Optional[float] = 0.0):
         u[hy_index, hy_index] = np.sin(
             (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
         )
-        try:
-            u = np.linalg.inv(u)
-        except np.linalg.LinAlgError:
-            u[hy_index, hx_index] = np.cos(
-                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-            )
-            u[hy_index, hy_index] = np.sin(
-                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-            )
-            u = np.linalg.inv(u)  # Identity if angle=0
+        u = np.linalg.inv(u)
 
         # build transformation matrix for predicted channels (electric fields)
         ex_index = self.channels.index("Ex")
diff --git a/aurora/transfer_function/compare.py b/aurora/transfer_function/compare.py
index 5e6cb814..a01483d4 100644
--- a/aurora/transfer_function/compare.py
+++ b/aurora/transfer_function/compare.py
@@ -286,6 +286,7 @@ def compare_transfer_functions(
         self,
         rtol: float = 1,
         atol: float = 1,
+        atol_phase: float = 4.0,
     ) -> dict:
         """
         Compare transfer functions between two transfer_functions objects.
@@ -299,6 +300,8 @@ def compare_transfer_functions(
             Relative tolerance for np.allclose, defaults to 1e-2
         atol: float
             Absolute tolerance for np.allclose, defaults to 1e-2
+        atol_phase: float
+            Absolute tolerance for phase comparison, defaults to 4.0 degrees
 
         Returns
         -------
@@ -357,7 +360,7 @@ def compare_transfer_functions(
                         np.angle(z1[:, ii, jj]),
                         np.angle(z2[:, ii, jj]),
                         rtol=rtol,
-                        atol=atol,
+                        atol=atol_phase,
                     )
 
                     result["impedance_error_close"] = np.allclose(
diff --git a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
index 1604d614..91c1dede 100644
--- a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
+++ b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
@@ -1,3 +1,18 @@
+"""
+Tests comparing aurora processing results against archived EMTF results.
+
+Development Notes:
+- In the early days of development these tests were useful to check that
+the same results were obtained from processing an mth5 with both stations in it
+and two separate mth5 files.  This could probably be made into a simpler test in mth5
+that checks that the data are the same in the two files.
+- This used to make a homebrew resistivity and phase plot for comparison between
+archived emtf z-files, but has been replaced with mt_metadata methods.
+- TODO: Check phases in these plots -- they are off by 180 so there may be a sign error
+in the data, or maybe the emtf results are using a different convention.  Need to investigate.
+- The comparison with the matlab emtf results uses a slighly different windowing method.
+
+"""
 import numpy as np
 import pytest
 from loguru import logger
@@ -6,18 +21,8 @@
 
 from aurora.general_helper_functions import DATA_PATH
 from aurora.pipelines.process_mth5 import process_mth5
-from aurora.sandbox.io_helpers.zfile_murphy import read_z_file
 from aurora.test_utils.synthetic.make_processing_configs import create_test_run_config
-from aurora.test_utils.synthetic.plot_helpers_synthetic import plot_rho_phi
-from aurora.test_utils.synthetic.rms_helpers import (
-    assert_rms_misfit_ok,
-    compute_rms,
-    get_expected_rms_misfit,
-)
 from aurora.transfer_function.compare import CompareTF
-from aurora.transfer_function.emtf_z_file_helpers import (
-    merge_tf_collection_to_match_z_file,
-)
 
 
 # Path to baseline EMTF results in source tree
@@ -31,9 +36,8 @@ def aurora_vs_emtf(
     auxilliary_z_file,
     z_file_base,
     tfk_dataset,
+    atol_phase=4.0,
     make_rho_phi_plot=True,
-    show_rho_phi_plot=False,
-    use_subtitle=True,
 ):
     """
     Compare aurora processing results against EMTF baseline.
@@ -62,55 +66,25 @@ def aurora_vs_emtf(
         test_case_id, tfk_dataset, matlab_or_fortran=emtf_version
     )
 
-    expected_rms_misfit = get_expected_rms_misfit(test_case_id, emtf_version)
     z_file_path = AURORA_RESULTS_PATH.joinpath(z_file_base)
 
-    tf_collection = process_mth5(
+    process_mth5(
         processing_config,
         tfk_dataset=tfk_dataset,
         z_file_path=z_file_path,
         return_collection=True,
     )
     comparator = CompareTF(tf_01=z_file_path, tf_02=auxilliary_z_file)
-    result = comparator.compare_transfer_functions()
+    result = comparator.compare_transfer_functions(atol_phase=4.0)
     assert np.isclose(result["impedance_ratio"]["Z_10"], 1.0, atol=1e-2)
     assert np.isclose(result["impedance_ratio"]["Z_01"], 1.0, atol=1e-2)
-    aux_data = read_z_file(auxilliary_z_file)
-    aurora_rho_phi = merge_tf_collection_to_match_z_file(aux_data, tf_collection)
-    data_dict = {}
-    data_dict["period"] = aux_data.periods
-    data_dict["emtf_rho_xy"] = aux_data.rxy
-    data_dict["emtf_phi_xy"] = aux_data.pxy
-    for xy_or_yx in ["xy", "yx"]:
-        aurora_rho = aurora_rho_phi["rho"][xy_or_yx]
-        aurora_phi = aurora_rho_phi["phi"][xy_or_yx]
-        aux_rho = aux_data.rho(xy_or_yx)
-        aux_phi = aux_data.phi(xy_or_yx)
-        rho_rms_aurora, phi_rms_aurora = compute_rms(
-            aurora_rho, aurora_phi, verbose=True
+    assert result["impedance_phase_close"]
+    if make_rho_phi_plot:
+        comparator.plot_two_transfer_functions(
+            save_plot_path=AURORA_RESULTS_PATH.joinpath(
+                f"{test_case_id}_aurora_vs_emtf_{emtf_version}_tf_compare.png"
+            ),
         )
-        rho_rms_emtf, phi_rms_emtf = compute_rms(aux_rho, aux_phi)
-        data_dict["aurora_rho_xy"] = aurora_rho
-        data_dict["aurora_phi_xy"] = aurora_phi
-        if expected_rms_misfit is not None:
-            assert_rms_misfit_ok(
-                expected_rms_misfit, xy_or_yx, rho_rms_aurora, phi_rms_aurora
-            )
-
-        if make_rho_phi_plot:
-            plot_rho_phi(
-                xy_or_yx,
-                tf_collection,
-                rho_rms_aurora,
-                rho_rms_emtf,
-                phi_rms_aurora,
-                phi_rms_emtf,
-                emtf_version,
-                aux_data=aux_data,
-                use_subtitle=use_subtitle,
-                show_plot=show_rho_phi_plot,
-                output_path=AURORA_RESULTS_PATH,
-            )
 
 
 @pytest.mark.slow

From 25e37c14067097d89b23dda1a4b3668a9cefc826 Mon Sep 17 00:00:00 2001
From: JP <peacock.jared@gmail.com>
Date: Sun, 25 Jan 2026 21:13:08 -0800
Subject: [PATCH 4/8] Update pip install commands in tests.yaml

update to pip install mt_metadata main and mth5 master
---
 .github/workflows/tests.yaml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index 1255ed2e..13189467 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -47,9 +47,9 @@ jobs:
           source .venv/bin/activate
           uv pip install --upgrade pip
           uv pip install -e ".[dev,test]"
-          # uv pip install mt_metadata[obspy]
-          uv pip install "mt_metadata[obspy] @ git+https://github.com/kujaku11/mt_metadata.git@patches"
-          uv pip install git+https://github.com/kujaku11/mth5.git@patches
+          uv pip install mt_metadata[obspy]
+          # uv pip install "mt_metadata[obspy] @ git+https://github.com/kujaku11/mt_metadata.git@patches"
+          uv pip install git+https://github.com/kujaku11/mth5.git@master
           
           # uv pip install mth5
           uv pip install git+https://github.com/kujaku11/mth5_test_data.git

From 48a0aa4968a032bdb3d036858cbe8e68cfc67e44 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Sun, 25 Jan 2026 21:36:20 -0800
Subject: [PATCH 5/8] remove ad-hoc rms checks - replaced with ratio checks and
 absolute angle checks

---
 .../synthetic/plot_helpers_synthetic.py       |   3 +
 aurora/test_utils/synthetic/rms_helpers.py    | 125 ------------------
 2 files changed, 3 insertions(+), 125 deletions(-)
 delete mode 100644 aurora/test_utils/synthetic/rms_helpers.py

diff --git a/aurora/test_utils/synthetic/plot_helpers_synthetic.py b/aurora/test_utils/synthetic/plot_helpers_synthetic.py
index e8a9a44a..c12ccaf8 100644
--- a/aurora/test_utils/synthetic/plot_helpers_synthetic.py
+++ b/aurora/test_utils/synthetic/plot_helpers_synthetic.py
@@ -8,6 +8,9 @@ def make_subtitle(
 ):
     """
 
+    Development Notes:
+    compute_rms is now deprecated.  If that tool is desired it should be reimplemented
+    in aurora/transfer_function/compare.py
     Parameters
     ----------
     rho_rms_aurora: float
diff --git a/aurora/test_utils/synthetic/rms_helpers.py b/aurora/test_utils/synthetic/rms_helpers.py
deleted file mode 100644
index 7a8f26b8..00000000
--- a/aurora/test_utils/synthetic/rms_helpers.py
+++ /dev/null
@@ -1,125 +0,0 @@
-"""
-    This module contains methods associated with RMS calculations that are used in testing
-    aurora processing on synthetic data.
-
-"""
-import numpy as np
-from loguru import logger
-
-
-def compute_rms(rho, phi, model_rho_a=100.0, model_phi=45.0, verbose=False):
-    """
-    Computes the RMS between processing results (rho, phi) and model (rho, phi).
-
-    It is used to make annotations for comparative plots for synthetic data. Could be
-    used in general to compare different processing results. For example by replacing
-    model_rho_a and model_phi with other processing results, or other (non-uniform)
-    model results.
-
-    Parameters
-    ----------
-    rho: numpy.ndarray
-        1D array of computed apparent resistivities (expected in Ohmm)
-    phi: numpy.ndarrayx
-        1D array of computed phases (expected in degrees)
-    model_rho_a: float or numpy array
-        if numpy array must be the same shape as rho
-    model_phi: float or numpy array
-        if numpy array must be the same shape as phi.
-    Returns
-    -------
-    rho_rms: float
-        rms misfit between the model apparent resistivity and the computed resistivity
-    phi_rms: float
-        rms misfit between the model phase (or phases) and the computed phase
-    """
-    rho_rms = np.sqrt(np.mean((rho - model_rho_a) ** 2))
-    phi_rms = np.sqrt(np.mean((phi - model_phi) ** 2))
-    if verbose:
-        logger.info(f"rho_rms = {rho_rms}")
-        logger.info(f"phi_rms = {phi_rms}")
-    return rho_rms, phi_rms
-
-
-def get_expected_rms_misfit(test_case_id: str, emtf_version=None) -> dict:
-    """
-    Returns hard-coded expected results from synthetic data processing.
-    These results are a benchmark against which test results are compared on push to
-    github.
-
-    Parameters
-    ----------
-    test_case_id
-    emtf_version
-
-    Returns
-    -------
-
-    """
-    expected_rms_misfit = {}
-    expected_rms_misfit["rho"] = {}
-    expected_rms_misfit["phi"] = {}
-    if test_case_id == "test1":
-        if emtf_version == "fortran":
-            expected_rms_misfit["rho"]["xy"] = 4.433905
-            expected_rms_misfit["phi"]["xy"] = 0.910484
-            expected_rms_misfit["rho"]["yx"] = 3.658614
-            expected_rms_misfit["phi"]["yx"] = 0.844645
-        elif emtf_version == "matlab":
-            expected_rms_misfit["rho"]["xy"] = 2.706098
-            expected_rms_misfit["phi"]["xy"] = 0.784229
-            expected_rms_misfit["rho"]["yx"] = 3.745280
-            expected_rms_misfit["phi"]["yx"] = 1.374938
-
-    elif test_case_id == "test2r1":
-        expected_rms_misfit["rho"]["xy"] = 3.971313
-        expected_rms_misfit["phi"]["xy"] = 0.982613
-        expected_rms_misfit["rho"]["yx"] = 3.967259
-        expected_rms_misfit["phi"]["yx"] = 1.62881
-    return expected_rms_misfit
-
-
-def assert_rms_misfit_ok(
-    expected_rms_misfit,
-    xy_or_yx,
-    rho_rms_aurora,
-    phi_rms_aurora,
-    rho_tol=1e-4,
-    phi_tol=1e-4,
-) -> None:
-    """
-    Compares actual RMS misfit from processing against expected values.
-    Raises Assertion errors if test processing results different from expected.
-
-    Parameters
-    ----------
-    expected_rms_misfit: dictionary
-        precomputed RMS misfits for test data in rho and phi
-    xy_or_yx: str
-        mode
-    rho_rms_aurora: float
-    phi_rms_aurora: float
-    """
-    expected_rms_rho = expected_rms_misfit["rho"][xy_or_yx]
-    expected_rms_phi = expected_rms_misfit["phi"][xy_or_yx]
-    logger.info(f"expected_rms_rho_{xy_or_yx} {expected_rms_rho}")
-    logger.info(f"expected_rms_phi_{xy_or_yx} {expected_rms_phi}")
-    if not np.isclose(rho_rms_aurora - expected_rms_rho, 0, atol=rho_tol):
-        logger.error("==== AURORA ====\n")
-        logger.error(rho_rms_aurora)
-        logger.error("==== EXPECTED ====\n")
-        logger.error(expected_rms_rho)
-        logger.error("==== DIFFERENCE ====\n")
-        logger.error(rho_rms_aurora - expected_rms_rho)
-        raise AssertionError("Expected misfit for resistivity is not correct")
-
-    if not np.isclose(phi_rms_aurora - expected_rms_phi, 0, atol=phi_tol):
-        logger.error("==== AURORA ====\n")
-        logger.error(phi_rms_aurora)
-        logger.error("==== EXPECTED ====\n")
-        logger.error(expected_rms_phi)
-        logger.error("==== DIFFERENCE ====\n")
-        logger.error(phi_rms_aurora - expected_rms_phi)
-        raise AssertionError("Expected misfit for phase is not correct")
-
-    return

From cf43dbcc2bf655de11628cd80de0e12a8b800f16 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Mon, 26 Jan 2026 08:48:08 -0800
Subject: [PATCH 6/8] adjust resisitivity axes limits in comparison plot

---
 aurora/transfer_function/compare.py              | 16 ++++++++++++++++
 ...est_compare_aurora_vs_archived_emtf_pytest.py |  2 ++
 2 files changed, 18 insertions(+)

diff --git a/aurora/transfer_function/compare.py b/aurora/transfer_function/compare.py
index a01483d4..0e81ad2a 100644
--- a/aurora/transfer_function/compare.py
+++ b/aurora/transfer_function/compare.py
@@ -77,6 +77,11 @@ def plot_two_transfer_functions(
         label_01="emtf",
         label_02="aurora",
         save_plot_path=None,
+        rho_xx_ylims=None,
+        rho_xy_ylims=None,
+        rho_yx_ylims=None,
+        rho_yy_ylims=None,
+        phi_ylims=None,
     ):
         """
         Plots two transfer functions for comparison.
@@ -94,10 +99,12 @@ def plot_two_transfer_functions(
         -------
 
         """
+        xy = "xy"
         fig = plt.figure(figsize=(12, 6))
 
         for ii in range(2):
             for jj in range(2):
+                component = f"{xy[ii]}{xy[jj]}"
                 plot_num_res = 1 + ii * 2 + jj
                 plot_num_phase = 5 + ii * 2 + jj
                 ax = fig.add_subplot(2, 4, plot_num_res)
@@ -127,6 +134,14 @@ def plot_two_transfer_functions(
                     ax.set_ylabel("Apparent Resistivity ($\Omega \cdot m$)")
                     ax.legend()
                 ax.grid(True, which="both", ls="--", lw=0.5, color="gray")
+                if component == "xx" and rho_xx_ylims is not None:
+                    ax.set_ylim(rho_xx_ylims)
+                if component == "xy" and rho_xy_ylims is not None:
+                    ax.set_ylim(rho_xy_ylims)
+                if component == "yx" and rho_yx_ylims is not None:
+                    ax.set_ylim(rho_yx_ylims)
+                if component == "yy" and rho_yy_ylims is not None:
+                    ax.set_ylim(rho_yy_ylims)
 
                 ax2 = fig.add_subplot(2, 4, plot_num_phase)
                 ax2.semilogx(
@@ -149,6 +164,7 @@ def plot_two_transfer_functions(
                 if plot_num_phase == 5:
                     ax2.set_ylabel("Phase (degrees)")
                     ax2.legend()
+
                 ax2.grid(True, which="both", ls="--", lw=0.5, color="gray")
 
         fig.tight_layout()
diff --git a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
index 91c1dede..65b41d4b 100644
--- a/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
+++ b/tests/synthetic/test_compare_aurora_vs_archived_emtf_pytest.py
@@ -84,6 +84,8 @@ def aurora_vs_emtf(
             save_plot_path=AURORA_RESULTS_PATH.joinpath(
                 f"{test_case_id}_aurora_vs_emtf_{emtf_version}_tf_compare.png"
             ),
+            rho_xy_ylims=(10.0, 1000.0),
+            rho_yx_ylims=(10.0, 1000.0),
         )
 
 

From 1bb252b66764ad0bb5236ab315269d8db15714d3 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Mon, 26 Jan 2026 08:51:46 -0800
Subject: [PATCH 7/8] remove deprecated z-file / tf_collection handling

---
 .../transfer_function/emtf_z_file_helpers.py  | 67 +------------------
 1 file changed, 2 insertions(+), 65 deletions(-)

diff --git a/aurora/transfer_function/emtf_z_file_helpers.py b/aurora/transfer_function/emtf_z_file_helpers.py
index a4fca770..031abdde 100644
--- a/aurora/transfer_function/emtf_z_file_helpers.py
+++ b/aurora/transfer_function/emtf_z_file_helpers.py
@@ -8,14 +8,10 @@
 
 """
 import pathlib
+from typing import Optional, Union
 
-import numpy as np
-from aurora.transfer_function.transfer_function_collection import (
-    TransferFunctionCollection,
-)
-from aurora.sandbox.io_helpers.zfile_murphy import ZFile
 from loguru import logger
-from typing import Optional, Union
+
 
 EMTF_CHANNEL_ORDER = ["hx", "hy", "hz", "ex", "ey"]
 
@@ -47,65 +43,6 @@ def get_default_orientation_block(n_ch: int = 5) -> list:
     return orientation_strs
 
 
-def merge_tf_collection_to_match_z_file(
-    aux_data: ZFile, tf_collection: TransferFunctionCollection
-) -> dict:
-    """
-    method to merge tf data from a tf_collection with a Z-file when there are potentially
-    multiple estimates of TF at the same periods for different decimation levels.
-
-    Development Notes:
-    Currently this is only used for the the synthetic test where aurora results
-    are compared against a stored legacy Z-file.  Given data from a z_file, and a
-    tf_collection, the tf_collection may have several TF estimates at the same
-    frequency from multiple decimation levels.  This tries to make a single array as
-    a function of period for all rho and phi.
-
-    Parameters
-    ----------
-    aux_data: aurora.sandbox.io_helpers.zfile_murphy.ZFile
-        Object representing a z-file
-    tf_collection: aurora.transfer_function.transfer_function_collection
-    .TransferFunctionCollection
-        Object representing the transfer function returned from the aurora processing
-
-
-    Returns
-    -------
-    result: dict of dicts
-        Keyed by ["rho", "phi"], below each of these is an ["xy", "yx",] entry.  The
-        lowest level entries are numpy arrays.
-    """
-    rxy = np.full(len(aux_data.decimation_levels), np.nan)
-    ryx = np.full(len(aux_data.decimation_levels), np.nan)
-    pxy = np.full(len(aux_data.decimation_levels), np.nan)
-    pyx = np.full(len(aux_data.decimation_levels), np.nan)
-    dec_levels = list(set(aux_data.decimation_levels))
-    dec_levels = [int(x) for x in dec_levels]
-    dec_levels.sort()
-
-    for dec_level in dec_levels:
-        aurora_tf = tf_collection.tf_dict[dec_level - 1]
-        indices = np.where(aux_data.decimation_levels == dec_level)[0]
-        for ndx in indices:
-            period = aux_data.periods[ndx]
-            # find the nearest period in aurora_tf
-            aurora_ndx = np.argmin(np.abs(aurora_tf.periods - period))
-            rxy[ndx] = aurora_tf.rho[aurora_ndx, 0]
-            ryx[ndx] = aurora_tf.rho[aurora_ndx, 1]
-            pxy[ndx] = aurora_tf.phi[aurora_ndx, 0]
-            pyx[ndx] = aurora_tf.phi[aurora_ndx, 1]
-
-    result = {}
-    result["rho"] = {}
-    result["phi"] = {}
-    result["rho"]["xy"] = rxy
-    result["phi"]["xy"] = pxy
-    result["rho"]["yx"] = ryx
-    result["phi"]["yx"] = pyx
-    return result
-
-
 def clip_bands_from_z_file(
     z_path: Union[str, pathlib.Path],
     n_bands_clip: int,

From bd8bf1d068982f928bdcd6e91fc46d1979664637 Mon Sep 17 00:00:00 2001
From: "Karl N. Kappler" <magnetotellurics@gmail.com>
Date: Mon, 26 Jan 2026 11:19:24 -0800
Subject: [PATCH 8/8] deprecate legacy z-file reader

---
 aurora/sandbox/io_helpers/zfile_murphy.py | 449 ----------------------
 1 file changed, 449 deletions(-)
 delete mode 100644 aurora/sandbox/io_helpers/zfile_murphy.py

diff --git a/aurora/sandbox/io_helpers/zfile_murphy.py b/aurora/sandbox/io_helpers/zfile_murphy.py
deleted file mode 100644
index d4bd27a3..00000000
--- a/aurora/sandbox/io_helpers/zfile_murphy.py
+++ /dev/null
@@ -1,449 +0,0 @@
-"""
-This module contains a class that was contributed by Ben Murphy for working with EMTF "Z-files"
-"""
-
-import pathlib
-import re
-from typing import Optional, Union
-
-import numpy as np
-
-
-class ZFile:
-    def __init__(self, filename: Union[str, pathlib.Path]):
-        """
-        Constructor
-
-        Parameters
-        ----------
-        filename: Union[str, pathlib.Path]
-            The path to the z-file.
-
-        """
-        self.filename = filename
-        self.station = ""
-        self.decimation_levels = None
-        self.lower_harmonic_indices = None
-        self.upper_harmonic_indices = None
-        self.f = None
-
-    def open_file(self) -> None:
-        """attempt to open file"""
-        try:
-            self.f = open(self.filename, "r")
-        except IOError:
-            raise IOError("File not found.")
-        return
-
-    def skip_header_lines(self) -> None:
-        """Skip over the header when reading"""
-        f = self.f
-        f.readline()
-        f.readline()
-        f.readline()
-        return
-
-    def get_station_id(self) -> None:
-        """get station ID from zfile"""
-        f = self.f
-        line = f.readline()
-        if line.lower().startswith("station"):
-            station = line.strip().split(":", 1)[1]
-        else:
-            station = line.strip()
-        self.station = station
-
-    def read_coordinates_and_declination(self) -> None:
-        """read coordinates and declination"""
-        f = self.f
-        line = f.readline().strip().lower()
-        match = re.match(
-            r"\s*coordinate\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+"
-            r"declination\s+(-?\d+\.?\d*)",
-            line,
-        )
-        self.coordinates = (float(match.group(1)), float(match.group(2)))
-        self.declination = float(match.group(3))
-        return
-
-    def read_number_of_channels_and_number_of_frequencies(self):
-        """read_number_of_channels_and_number_of_frequencies"""
-        f = self.f
-        line = f.readline().strip().lower()
-        match = re.match(
-            r"\s*number\s+of\s+channels\s+(\d+)\s+number\s+of"
-            r"\s+frequencies\s+(\d+)",
-            line,
-        )
-        nchannels = int(match.group(1))
-        nfreqs = int(match.group(2))
-        self.nchannels = nchannels
-        self.nfreqs = nfreqs
-
-    def read_channel_information(self):
-        """read_channel_information"""
-        f = self.f
-        self.orientation = np.zeros((self.nchannels, 2))
-        self.channels = []
-        for i in range(self.nchannels):
-            line = f.readline().strip()
-            match = re.match(
-                r"\s*\d+\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+\w*\s+" r"(\w+)", line
-            )
-            self.orientation[i, 0] = float(match.group(1))
-            self.orientation[i, 1] = float(match.group(2))
-            if len(match.group(3)) > 2:
-                # sometimes the channel ID comes out with extra stuff
-                self.channels.append(match.group(3)[:2].title())
-            else:
-                self.channels.append(match.group(3).title())
-        return
-
-    def load(self):
-        """load Z-file and populate attributes of class"""
-        self.open_file()
-        self.skip_header_lines()
-        self.get_station_id()
-        self.read_coordinates_and_declination()
-        self.read_number_of_channels_and_number_of_frequencies()
-
-        f = self.f
-
-        # skip line
-        f.readline()
-        self.read_channel_information()
-        f.readline()
-
-        # initialize empty arrays for transfer functions
-        # note that EMTF, and consequently this code, assumes two predictor
-        #    channels (horizontal magnetics)
-        # nchannels - 2 therefore is the number of predicted channels
-        self.decimation_levels = np.zeros(self.nfreqs)
-        self.periods = np.zeros(self.nfreqs)
-        self.lower_harmonic_indices = np.zeros(self.nfreqs)
-        self.upper_harmonic_indices = np.zeros(self.nfreqs)
-        self.transfer_functions = np.zeros(
-            (self.nfreqs, self.nchannels - 2, 2), dtype=np.complex64
-        )
-
-        # residual covariance -- square matrix with dimension as number of
-        # predicted channels
-        self.sigma_e = np.zeros(
-            (self.nfreqs, self.nchannels - 2, self.nchannels - 2), dtype=np.complex64
-        )
-
-        # inverse coherent signal power -- square matrix, with dimension as the
-        #    number of predictor channels
-        # since EMTF and this code assume N predictors is 2,
-        #    this dimension is hard-coded
-        self.sigma_s = np.zeros((self.nfreqs, 2, 2), dtype=np.complex64)
-
-        # now read data for each period
-        for i in range(self.nfreqs):
-            # extract period
-            line = f.readline().strip()
-            match = re.match(
-                r"\s*period\s*:\s+(\d+\.?\d*)\s+" r"decimation\s+level", line
-            )
-            self.periods[i] = float(match.group(1))
-
-            splitted_line1 = line.split("level")
-            splitted_line2 = splitted_line1[1].split("freq")
-            self.decimation_levels[i] = int(splitted_line2[0].strip())
-            splitted_line1 = line.split("from")
-            splitted_line2 = splitted_line1[1].split("to")
-            self.lower_harmonic_indices[i] = int(splitted_line2[0].strip())
-            self.upper_harmonic_indices[i] = int(splitted_line2[1].strip())
-            # skip two lines
-            f.readline()
-            f.readline()
-
-            # read transfer functions
-            for j in range(self.nchannels - 2):
-                comp1_r, comp1_i, comp2_r, comp2_i = f.readline().strip().split()
-                self.transfer_functions[i, j, 0] = float(comp1_r) + 1.0j * float(
-                    comp1_i
-                )
-                self.transfer_functions[i, j, 1] = float(comp2_r) + 1.0j * float(
-                    comp2_i
-                )
-
-            # skip label line
-            f.readline()
-
-            # read inverse coherent signal power matrix
-            val1_r, val1_i = f.readline().strip().split()
-            val2_r, val2_i, val3_r, val3_i = f.readline().strip().split()
-            self.sigma_s[i, 0, 0] = float(val1_r) + 1.0j * float(val1_i)
-            self.sigma_s[i, 1, 0] = float(val2_r) + 1.0j * float(val2_i)
-            self.sigma_s[i, 0, 1] = float(val2_r) - 1.0j * float(val2_i)
-            self.sigma_s[i, 1, 1] = float(val3_r) + 1.0j * float(val3_i)
-
-            # skip label line
-            f.readline()
-
-            # read residual covariance
-            for j in range(self.nchannels - 2):
-                values = f.readline().strip().split()
-                for k in range(j + 1):
-                    if j == k:
-                        self.sigma_e[i, j, k] = float(values[2 * k]) + 1.0j * float(
-                            values[2 * k + 1]
-                        )
-                    else:
-                        self.sigma_e[i, j, k] = float(values[2 * k]) + 1.0j * float(
-                            values[2 * k + 1]
-                        )
-                        self.sigma_e[i, k, j] = float(values[2 * k]) - 1.0j * float(
-                            values[2 * k + 1]
-                        )
-
-        f.close()
-
-    def impedance(self, angle: Optional[float] = 0.0):
-        """
-        Compute the Impedance and errors from the transfer function.
-        - note u,v are identity matrices if angle=0
-
-        Parameters
-        ----------
-        angle: float
-            The angle about the vertical axis by which to rotate the Z tensor.
-
-        Returns
-        -------
-        z: np.ndarray
-            The impedance tensor
-        error: np.ndarray
-            Errors for the impedance tensor
-        """
-        # check to see if there are actually electric fields in the TFs
-        if "Ex" not in self.channels and "Ey" not in self.channels:
-            raise ValueError(
-                "Cannot return apparent resistivity and phase "
-                "data because these TFs do not contain electric "
-                "fields as a predicted channel."
-            )
-
-        # transform the TFs first...
-        # build transformation matrix for predictor channels
-        #    (horizontal magnetic fields)
-        hx_index = self.channels.index("Hx")
-        hy_index = self.channels.index("Hy")
-        u = np.eye(2, 2)
-        u[hx_index, hx_index] = np.cos(
-            (self.orientation[hx_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hx_index, hy_index] = np.sin(
-            (self.orientation[hx_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hy_index, hx_index] = np.sin(
-            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hy_index, hy_index] = np.cos(
-            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-        )
-        try:
-            u = np.linalg.inv(u)
-        except np.linalg.LinAlgError:
-            u[hy_index, hx_index] = np.cos(
-                (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-            )
-        u[hy_index, hy_index] = np.sin(
-            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-        )
-        u = np.linalg.inv(u)
-
-        # build transformation matrix for predicted channels (electric fields)
-        ex_index = self.channels.index("Ex")
-        ey_index = self.channels.index("Ey")
-        v = np.eye(self.transfer_functions.shape[1], self.transfer_functions.shape[1])
-        v[ex_index - 2, ex_index - 2] = np.cos(
-            (self.orientation[ex_index, 0] - angle) * np.pi / 180.0
-        )
-        v[ey_index - 2, ex_index - 2] = np.sin(
-            (self.orientation[ex_index, 0] - angle) * np.pi / 180.0
-        )
-        v[ex_index - 2, ey_index - 2] = np.cos(
-            (self.orientation[ey_index, 0] - angle) * np.pi / 180.0
-        )
-        v[ey_index - 2, ey_index - 2] = np.sin(
-            (self.orientation[ey_index, 0] - angle) * np.pi / 180.0
-        )
-
-        # matrix multiplication...
-        rotated_transfer_functions = np.matmul(
-            v, np.matmul(self.transfer_functions, u.T)
-        )
-        rotated_sigma_s = np.matmul(u, np.matmul(self.sigma_s, u.T))
-        rotated_sigma_e = np.matmul(v, np.matmul(self.sigma_e, v.T))
-
-        # now pull out the impedance tensor
-        z = np.zeros((self.periods.size, 2, 2), dtype=np.complex64)
-        z[:, 0, 0] = rotated_transfer_functions[:, ex_index - 2, hx_index]  # Zxx
-        z[:, 0, 1] = rotated_transfer_functions[:, ex_index - 2, hy_index]  # Zxy
-        z[:, 1, 0] = rotated_transfer_functions[:, ey_index - 2, hx_index]  # Zyx
-        z[:, 1, 1] = rotated_transfer_functions[:, ey_index - 2, hy_index]  # Zyy
-
-        # and the variance information
-        var = np.zeros((self.periods.size, 2, 2))
-        var[:, 0, 0] = np.real(
-            rotated_sigma_e[:, ex_index - 2, ex_index - 2]
-            * rotated_sigma_s[:, hx_index, hx_index]
-        )
-
-        var[:, 0, 1] = np.real(
-            rotated_sigma_e[:, ex_index - 2, ex_index - 2]
-            * rotated_sigma_s[:, hy_index, hy_index]
-        )
-        var[:, 1, 0] = np.real(
-            rotated_sigma_e[:, ey_index - 2, ey_index - 2]
-            * rotated_sigma_s[:, hx_index, hx_index]
-        )
-        var[:, 1, 1] = np.real(
-            rotated_sigma_e[:, ey_index - 2, ey_index - 2]
-            * rotated_sigma_s[:, hy_index, hy_index]
-        )
-
-        error = np.sqrt(var)
-
-        return z, error
-
-    def tippers(self, angle=0.0):
-        """compute the tipper from transfer function"""
-
-        # check to see if there is a vertical magnetic field in the TFs
-        if "Hz" not in self.channels:
-            raise ValueError(
-                "Cannot return tipper data because the TFs do not "
-                "contain the vertical magnetic field as a "
-                "predicted channel."
-            )
-
-        # transform the TFs first...
-        # build transformation matrix for predictor channels
-        #    (horizontal magnetic fields)
-        hx_index = self.channels.index("Hx")
-        hy_index = self.channels.index("Hy")
-        u = np.eye(2, 2)
-        u[hx_index, hx_index] = np.cos(
-            (self.orientation[hx_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hx_index, hy_index] = np.sin(
-            (self.orientation[hx_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hy_index, hx_index] = np.cos(
-            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-        )
-        u[hy_index, hy_index] = np.sin(
-            (self.orientation[hy_index, 0] - angle) * np.pi / 180.0
-        )
-        u = np.linalg.inv(u)
-
-        # don't need to transform predicated channels (assuming no tilt in Hz)
-        hz_index = self.channels.index("Hz")
-        v = np.eye(self.transfer_functions.shape[1], self.transfer_functions.shape[1])
-
-        # matrix multiplication...
-        rotated_transfer_functions = np.matmul(
-            v, np.matmul(self.transfer_functions, u.T)
-        )
-        rotated_sigma_s = np.matmul(u, np.matmul(self.sigma_s, u.T))
-        rotated_sigma_e = np.matmul(v, np.matmul(self.sigma_e, v.T))
-
-        # now pull out tipper information
-        tipper = np.zeros((self.periods.size, 2), dtype=np.complex64)
-        tipper[:, 0] = rotated_transfer_functions[:, hz_index - 2, hx_index]  # Tx
-        tipper[:, 1] = rotated_transfer_functions[:, hz_index - 2, hy_index]  # Ty
-
-        # and the variance/error information
-        var = np.zeros((self.periods.size, 2))
-        var[:, 0] = np.real(
-            rotated_sigma_e[:, hz_index - 2, hz_index - 2]
-            * rotated_sigma_s[:, hx_index, hx_index]
-        )  # Tx
-        var[:, 1] = np.real(
-            rotated_sigma_e[:, hz_index - 2, hz_index - 2]
-            * rotated_sigma_s[:, hy_index, hy_index]
-        )  # Ty
-        error = np.sqrt(var)
-
-        return tipper, error
-
-    def apparent_resistivity(self, angle: float = 0.0):
-        """compute the apparent resistivity from the impedance."""
-        z_tensor, error = self.impedance(angle=angle)
-        Zxy = z_tensor[:, 0, 1]
-        Zyx = z_tensor[:, 1, 0]
-        T = self.periods
-        self.rxy = T * (abs(Zxy) ** 2) / 5.0
-        self.ryx = T * (abs(Zyx) ** 2) / 5.0
-        self.pxy = np.rad2deg(np.arctan(np.imag(Zxy) / np.real(Zxy)))
-        self.pyx = np.rad2deg(np.arctan(np.imag(Zyx) / np.real(Zyx)))
-        return
-
-    def rho(self, mode):
-        """
-        Return the apparent resistivity for the given mode.
-
-        Convenience function to help with streamlining synthetic tests - to be
-        eventually replaced by functionality in mt_metadata.tf
-
-        Parameters
-        ----------
-        mode: str
-            "xy" or "yx"
-
-        Returns
-        -------
-        rho
-        """
-        if mode == "xy":
-            return self.rxy
-        if mode == "yx":
-            return self.ryx
-
-    def phi(self, mode):
-        """
-        Return the phase for the given mode.
-
-        Convenience function to help with streamlining synthetic tests - to be
-        eventually replaced by functionality in mt_metadata.tf
-        Parameters
-        ----------
-        mode: str
-            "xy" or "yx"
-
-        Returns
-        -------
-        phi
-        """
-        if mode == "xy":
-            return self.pxy
-        if mode == "yx":
-            return self.pyx
-
-
-def read_z_file(z_file_path, angle=0.0) -> ZFile:
-    """
-    Reads a zFile and returns a ZFile object.
-
-    Parameters
-    ----------
-    z_file_path: string or pathlib.Path
-        The name of the EMTF-style z-file to operate on
-    angle: float
-        How much rotation to apply.  This is a kludge variable used to help compare
-        legacy SPUD results which are rotated onto a cardinal grid, vs aurora which
-        store the TF in the coordinate system of acquisition
-
-    Returns
-    -------
-    z_obj: ZFile
-        The zFile as an object.
-
-    """
-    z_obj = ZFile(z_file_path)
-    z_obj.load()
-    z_obj.apparent_resistivity(angle=angle)
-    return z_obj