Merge pull request #290 from SuperDARN/ehn/tdiff-correction

EHN: Tdiff Correction for Elevation Data

pydarn/__init__.py

@@ -56,6 +56,7 @@
 from .utils.scan import build_scan
 from .utils.geo import geocentric_coordinates
 from .utils.coordinates import Coords
+from .utils.recalculate_elevation import recalculate_elevation
 from .utils.filters import Boxcar
 
 # import plotting

pydarn/utils/recalculate_elevation.py

@@ -0,0 +1,142 @@
+# (C) Copyright SuperDARN Canada, University of Saskatchewan 2022
+# Author: Carley Martin 20221101
+#
+# Disclaimer:
+# pyDARN is under the LGPL v3 license found in the root directory LICENSE.md
+# Everyone is permitted to copy and distribute verbatim copies of this license
+# document, but changing it is not allowed.
+#
+# This version of the GNU Lesser General Public License incorporates the terms
+# and conditions of version 3 of the GNU General Public License,
+# supplemented by the additional permissions listed below.
+#
+# About:
+# Modified from the RST elevation_v2.c code (Simon Shepherd's code).
+# Method: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017RS006348
+# THIS CODE IS DESIGNED TO AMEND ELEVATION
+# AFTER FITACF FITTING (AT THE VISUALIZATION STEP TO QUICKLY CHECK THE EFFECTS
+# OF DIFFERENT TDIFF VALUES) NOT TO REPLACE THE FITTING ITSELF
+#
+# Modifications:
+# 20221221 - Bharat Kunduri: Updated to RST elevation code
+
+import numpy as np
+
+from typing import List
+from copy import deepcopy
+
+from pydarn import (SuperDARNRadars, C)
+
+
+def recalculate_elevation(dmap_data: List[dict], tdiff: float,
+                          overwrite: bool = False,
+                          interferometer_offset: list = None):
+    """
+    Recalculates elevation values for a given tdiff Value
+
+    Parameters
+    -----------
+    dmap_data: list of dictionaries
+        fitacf data
+    tdiff: float
+        propagation time from interferometer array to phasing matrix
+        input minus propagation time from main array antenna, microseconds
+    overwrite: bool
+        If true then return a new dmap_data with new elevation to plot with
+        if false then return dictionary of new elevations for further use
+    interferometer_offset: list
+        select position of interferometer array wrt the main array
+        needs to be list of [X, Y, Z] e.g. [0.0, 100.0, 1.0]
+
+    Raises
+    ------
+
+    Returns
+    -------
+    elv_amended: dictionary of lists
+        amended elevation values for each record given
+    """
+    if not overwrite:
+        # Make phi1 output into dictionary
+        elv_amended = {}
+    else:
+        dmap_amended = deepcopy(dmap_data)
+
+    # Hardware config for radar
+    # Doesn't have to be in the loop, since we're accessing only the first rec
+    radar_hdw = SuperDARNRadars.radars[dmap_data[0]['stid']].hardware_info
+    if interferometer_offset is not None:
+        int_pos = interferometer_offset
+    else:
+        int_pos = radar_hdw.interferometer_offset
+
+    # If inteferometer is infront (+1) or behind (-1) main array
+    if int_pos[1] < 0:
+        sgn = -1.0
+    else:
+        sgn = 1.0
+
+    boff = (radar_hdw.beams / 2.0) - 0.5
+
+    # For each record in data, recalculate the elevation
+    for ind in range(0, len(dmap_data)):
+        # If there is no elevation data, append empty list and skip
+        if 'phi0' not in dmap_data[ind]:
+            print("No elevation data. 'phi0' parameter missing"
+                  " from the record")
+            if not overwrite:
+                elv_amended[ind] = []
+            continue
+
+        # Beam direction off boresight in RADIANS
+        phi0 = np.radians(radar_hdw.beam_separation
+                          * (dmap_data[ind]['bmnum'] - boff))
+        # Cos and Sin of phi in shape of phi0
+        cp0 = np.cos(phi0)
+        sp0 = np.sin(phi0)
+
+        # Phase delay [radians] due to electrical path difference.
+        psi_ele = (-2.0 * np.pi * (dmap_data[ind]['tfreq'])
+                   * 1000.0 * tdiff * 1.0e-6)
+        # Elevation angle (a0) where psi (phase difference) is maximum
+        a0 = np.arcsin(sgn * int_pos[2] * cp0
+                       / np.sqrt(int_pos[1]**2 + int_pos[2]**2))
+        if a0 < 0:
+            a0 = 0
+        ca0 = np.cos(a0)
+        sa0 = np.sin(a0)
+
+        # maximum phase = psi_ele + psi_geo(a0)
+        psi_max = psi_ele + 2.0 * np.pi * (dmap_data[ind]['tfreq']) *\
+            (1e3 / C) * (int_pos[0] * sp0 + int_pos[1]
+                         * np.sqrt(ca0*ca0 - sp0*sp0) + int_pos[2] * sa0)
+
+        # compute the number of 2pi factors necessary to map to correct region
+        dpsi = (psi_max - dmap_data[ind]['phi0'])
+        if int_pos[1] > 0:
+            n2pi = np.floor(dpsi / (2.0 * np.pi))
+        else:
+            n2pi = np.ceil(dpsi / (2.0 * np.pi))
+        d2pi = n2pi * 2.0 * np.pi
+        # map observed phase to correct extended phase
+        psi_obs = dmap_data[ind]['phi0'] + d2pi
+        # solve for the elevation angle
+        E = (psi_obs / (2.0*np.pi*dmap_data[ind]['tfreq']*1.0e3)
+             + tdiff*1e-6) * C - int_pos[0] * sp0
+
+        alpha = np.arcsin((E*int_pos[2]
+                           + np.sqrt(E*E * int_pos[2]**2
+                           - (int_pos[1]**2 + int_pos[2]**2)
+                           * (E*E - int_pos[1]*int_pos[1]*cp0*cp0)))
+                          / (int_pos[1]*int_pos[1] + int_pos[2]*int_pos[2]))
+
+        # Convert theta back to degrees
+        if overwrite:
+            dmap_amended[ind].update({'elv': np.degrees(alpha)})
+        else:
+            elv_amended[ind] = np.degrees(alpha)
+
+    if overwrite:
+        return dmap_amended
+    else:
+        return elv_amended