Add files via upload

Author: gsteranka

src/norm_diff_class.py

@@ -0,0 +1,257 @@
+#!/usr/bin/env python
+"""
+
+norm_diff_class.py
+
+Purpose: Apply previously made calibration file to calibrate raw data
+         (frequency correct and normalize), which produces the diffraction
+         pattern given to Fresnel Inversion code.
+
+Revisions:
+      gjs_produce_normalized_diffraction_pattern.py
+   Mar 07 2018 - gsteranka - Original version. Tentative rough draft, since
+                             haven't yet chosen format of calibration file, and
+                             also don't have Ryan's code to find window width
+                             from a given window type and resolution.
+   Mar 19 2018 - gsteranka - Edited to take new tentative final form of cal
+                             file, which has the columns as 1) SPM 2) Predicted
+                             sky frequency 3) Residual frequency fit
+                             4) Freespace power spline fit 5) Frequency offset
+                             fit. Also edited to save an obs file and return
+                             everything in the obs file too
+      produce_normalized_diffraction_pattern.py
+   Mar 20 2018 - gsteranka - Eliminate debug steps, and switched obs_file to be
+                             optional inputs
+      norm_diff_class.py
+   Mar 22 2018 - gsteranka - Edit previous code to be a class, with Fresnel
+                            inversion inputs as attributes
+   Mar 30 2018 - gsteranka - Edited reading in of geometry file to read in csv
+                             format
+   Apr 03 2018 - gsteranka - When first defining p_norm_vals, do it by
+                             normalizing IQ_c_desired, not IQ_m_desired
+   Apr 09 2018 - gsteranka - Obs files now save absolute value of rho dot,
+                             which eliminates a lot of error checking in
+                             Fresnel Inversion step. Added
+                             "fill_value='extrapolate'" to rho_interp_func,
+                             p_free_interp_func, spm_desired_func, and
+                             f_sky_pred_func
+"""
+
+import numpy as np
+import pandas as pd
+from scipy.interpolate import interp1d
+import spiceypy as spice
+import sys
+
+from rsr_reader import RSRReader
+from resample_IQ import resample_IQ
+
+
+class NormDiff(object):
+    """Produce a normalized diffraction pattern using a previously made
+    calibration file (cal_file). Contains everything needed for Fresnel
+    Inversion as attributes.
+
+    Attributes:
+        rho_km_desired (np.ndarray): Ring intercept radius values, in km,
+            at final spacing specified in dr_km
+        spm_desired (np.ndarray): Set of SPM values corresponding to
+            rho_km_desired
+        p_norm_vals (np.ndarray): Power normalized to 1. This is the diffraction
+            pattern that is input to the Fresnel inversion step
+        phase_rad_vals (np.ndarray): Phase of the complex signal, in radians.
+            This is the other part of the diffraction pattern that is input to
+            the Fresnel Inversion step
+        B_rad_vals (np.ndarray): Ring opening angle of Saturn
+        D_km_vals (np.ndarray): Spacecraft-RIP (Ring Intercept Point) distance
+        F_km_vals (np.ndarray): Fresnel scale, in km
+        f_sky_pred_vals (np.ndarray): Predicted sky frequency, in Hz.
+        phi_rad_vals (np.ndarray): Observed ring azimuth, in radians
+        rho_dot_kms_vals (np.ndarray): Ring intercept point velocity
+    """
+
+    def __init__(self, rsr_file, dr_km, rho_km_range, res_km, window_type,
+             geo_file, cal_file, dr_km_tol=0.01,
+             decimate_16khz_to_1khz=False):
+        """
+        Instantiation defines all attributes of instance.
+
+        Args:
+            rsr_file (str): Full path name of RSR file to process
+            dr_km (float): Desired final radial spacing of processed data, in km
+            rho_km_range (list): Radius range to process over, in km
+            res_km (float): Resolution to later do Fresnel Inversion at, in km
+            window_type (str): Window type to later do Fresnel Inversion at
+            geo_file (str): Full path name of a geoemtry file corresponding to
+            rsr_file
+            cal_file (str): Full path name of a calibration file corresponding to
+                rsr_file
+            dr_km_tol (float): Optional keyword argument, in km, that specifies the
+                maximum distance the starting point of the final set of radius
+                values can be away from an integer number of dr_km. For example, if
+                you say dr_km_tol=0.01 and dr_km=0.25, the starting point of the
+                final set of rho values could be 70000.26 km, but not 70000.261 km
+            decimate_16khz_to_1khz (bool): Optional boolean keyword argument that
+                specifies whether or not to decimate a 16khz sample rate file down
+                to 1khz sample rate. Can only specify when rsr_file is 16khz. Value
+                of this should match what you said in RSRfile.get_IQ method when you
+                initially made the calibration file. Defualt is False
+        """
+
+        # Magical code to find radius range necessary for Fresnel Inversion.
+        # For now just use the input radius range
+        rho_km_range_inversion = rho_km_range
+
+        # Read geometry file for radius-SPM conversion
+        geo = pd.read_csv(geo_file, header=None)
+        spm_geo = np.asarray(geo[0][:])
+        rho_km_geo = np.asarray(geo[3][:])
+        B_rad_geo = np.asarray(geo[6][:]) * spice.rpd()
+        D_km_geo = np.asarray(geo[7][:])
+        phi_rad_geo = np.asarray(geo[5][:]) * spice.rpd()
+        rho_dot_kms_geo = np.asarray(geo[8][:])
+        F_km_geo = np.asarray(geo[10][:])
+
+        # Function to convert SPM to rho
+        rho_interp_func = interp1d(spm_geo, rho_km_geo, fill_value='extrapolate')
+
+        # Find SPM range corresponding to rho range. Bottom-most line of this
+        # block is to handle ingress files
+        spm_start = max(spm_geo[rho_km_geo <= rho_km_range[0]])
+        spm_end = min(spm_geo[rho_km_geo >= rho_km_range[1]])
+        spm_range = [spm_start, spm_end]
+        spm_range = [min(spm_range), max(spm_range)]
+
+        rsr_inst = RSRReader(rsr_file)
+        (spm_vals, IQ_m) = rsr_inst.get_IQ(spm_range=spm_range,
+                                           decimate_16khz_to_1khz=
+                                           decimate_16khz_to_1khz)
+        rho_km_vals = rho_interp_func(spm_vals)
+
+        cal = pd.read_csv(cal_file, delim_whitespace=True, header=None)
+
+        spm_cal = np.asarray(cal[0][:])
+        f_sky_pred_cal = np.asarray(cal[1][:])
+        p_free_cal = np.asarray(cal[3][:])
+        f_offset_fit_cal = np.asarray(cal[4][:])
+        rho_km_cal = rho_interp_func(spm_cal)
+
+        # If specified rho range that cal file can't cover
+        if (spm_range[0]<min(spm_cal)) | (spm_range[1]>max(spm_cal)):
+            print('ERROR (NormDiff): Specified cal file is missing \n'
+                  +'points required to pre-process points in rho_km_range.\n'
+                  +'Calibration file rho range is '+str(min(rho_km_cal))+'km, '
+                  +str(max(rho_km_cal))+'km, while specified rho_km_range is \n'
+                  +str(rho_km_range[0])+'km, '+str(rho_km_range[1])+'km')
+            sys.exit()
+
+        # Inteprolate frequency offset to finer spacing in preparation
+        # for integration
+        dt = 0.1
+        n_pts = round((spm_cal[-1] - spm_cal[0])/dt)
+        spm_cal_interp = spm_cal[0] + dt*np.arange(n_pts)
+        f_offset_fit_func = interp1d(spm_cal, f_offset_fit_cal,
+                                     fill_value='extrapolate')
+        f_offset_fit_interp = f_offset_fit_func(spm_cal_interp)
+
+        # Integrate frequency offset to detrend IQ_m
+        f_detrend_interp = np.cumsum(f_offset_fit_interp)*dt
+        f_detrend_interp_rad = f_detrend_interp * spice.twopi()
+        f_detrend_rad_func = interp1d(spm_cal_interp, f_detrend_interp_rad,
+                                      fill_value='extrapolate')
+        f_detrend_rad = f_detrend_rad_func(spm_vals)
+
+        # Detrend raw measured I and Q
+        IQ_c = IQ_m*np.exp(-1j*f_detrend_rad)
+
+        (rho_km_desired, IQ_c_desired) = resample_IQ(rho_km_vals, IQ_c, dr_km,
+                                                     dr_km_tol=dr_km_tol)
+
+        # Interpolate freespace power (spline fit) to rho_km_desired
+        p_free_interp_func = interp1d(rho_km_cal, p_free_cal,
+                                      fill_value='extrapolate')
+        p_free = p_free_interp_func(rho_km_desired)
+
+        # Construct power from IQ_c_desired and not IQ_m_desired or resampled
+        # power, because p_free was constructed from downsampled IQ_c. It's
+        # important to not use IQ_m_desired because I and Q each vary
+        # sinusoidally, meaning that if you resample too much, they get
+        # averaged to near zero.
+        p_norm_vals = (abs(IQ_c_desired)**2)/p_free
+        phase_rad_vals = np.arctan2(np.imag(IQ_c_desired),
+                                    np.real(IQ_c_desired))
+
+        spm_desired_func = interp1d(rho_km_geo, spm_geo,
+                                    fill_value='extrapolate')
+        spm_desired = spm_desired_func(rho_km_desired)
+
+        self.__set_attributes(rho_km_desired, spm_desired, p_norm_vals,
+                              phase_rad_vals,
+                              rho_km_geo, B_rad_geo, D_km_geo, F_km_geo,
+                              phi_rad_geo, rho_dot_kms_geo,
+                              rho_km_cal, f_sky_pred_cal)
+
+
+    def __set_attributes(self, rho_km_desired, spm_desired, p_norm_vals,
+                         phase_rad_vals,
+                         rho_km_geo, B_rad_geo, D_km_geo, F_km_geo, phi_rad_geo,
+                         rho_dot_kms_geo,
+                         rho_km_cal, f_sky_pred_cal):
+        """Private method called by __init__ to set attributes of the
+        instance"""
+
+        self.rho_km_desired = rho_km_desired
+        self.spm_desired = spm_desired
+        self.p_norm_vals = p_norm_vals
+        self.phase_rad_vals = phase_rad_vals
+
+        # Ring opening angle at final spacing
+        B_rad_func = interp1d(rho_km_geo, B_rad_geo, fill_value='extrapolate')
+        self.B_rad_vals = B_rad_func(rho_km_desired)
+
+        # Spacecraft - RIP distance at final spacing
+        D_km_func = interp1d(rho_km_geo, D_km_geo, fill_value='extrapolate')
+        self.D_km_vals = D_km_func(rho_km_desired)
+
+        # Fresnel scale at final spacing
+        F_km_func = interp1d(rho_km_geo, F_km_geo, fill_value='extrapolate')
+        self.F_km_vals = F_km_func(rho_km_desired)
+
+        # Sky frequency at final spacing
+        f_sky_pred_func = interp1d(rho_km_cal, f_sky_pred_cal,
+                                   fill_value='extrapolate')
+        self.f_sky_pred_vals = f_sky_pred_func(rho_km_desired)
+
+        # Observed Ring Azimuth at final spacing
+        phi_rad_func = interp1d(rho_km_geo, phi_rad_geo, fill_value='extrapolate')
+        self.phi_rad_vals = phi_rad_func(rho_km_desired)
+
+        # RIP velocity at final spacing
+        rho_dot_kms_func = interp1d(rho_km_geo, rho_dot_kms_geo,
+                                    fill_value='extrapolate')
+        self.rho_dot_kms_vals = rho_dot_kms_func(rho_km_desired)
+
+
+    def save_obs(self, obs_file):
+        """Save an obs file, which contains all of the attributes in the
+        instance, and therefore has everything needed for Fresnel Inversion
+
+        Args:
+            obs_file (str): Full path name of obs file to be created
+        """
+
+        # Save obs file with absolute value of rho dot. With this  and
+        # increasing rho_km_desired, ingress files can be put into
+        # inversion code without any error checks to differentiate
+        # from  egress
+        np.savetxt(obs_file,
+                   np.c_[self.rho_km_desired, self.spm_desired,
+                         self.p_norm_vals, self.phase_rad_vals,
+                         self.B_rad_vals, self.D_km_vals, self.F_km_vals,
+                         self.f_sky_pred_vals, self.phi_rad_vals,
+                         abs(self.rho_dot_kms_vals)],
+                   fmt='%32.16f '*10)
+
+
+if __name__ == '__main__':
+    pass

src/resample_IQ.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+"""
+
+resample_IQ.py
+
+Purpose: Resample I and Q from uniformly spaced time to uniformly spaced
+         radius. Normally, you resample from raw resolution data.
+
+Revisions:
+      gjs_resample_IQ.py
+   Feb 27 2018 - gsteranka - Original version
+   Mar 01 2018 - gsteranka - Took pre_resample function outside of resample_IQ
+                             function, since I found out it doesn't need to be
+                             inside to be recognized in a call to resample_IQ
+                             from another file
+   Mar 07 2018 - gsteranka - Input is a complex number IQ_c, instead of
+                             separate real and imaginary parts
+   Mar 08 2018 - gsteranka - Change input arrays so rho starts within dr_km_tol
+                             of an integer number of dr_km. Makes it so after
+                             resampling, everything is at a more even number
+      resample_IQ.py
+   Mar 20 2018 - gsteranka - Copy to official version and remove debug steps
+   Apr 09 2018 - gsteranka - Edited pre_resample function to handle ingerss
+                             files as well as egress
+"""
+
+import numpy as np
+from scipy import signal
+from scipy.interpolate import interp1d
+
+def pre_resample(rho_km, vec, freq):
+    """Sub-function inside of resample_IQ to put vectors to uniform spaced
+    radius at a similar spacing as raw resolution.  For ingress occultations,
+    this step implicitly reverses the radius scale when interpolating"""
+
+    # Average radius spacing over region
+    ts_avg = abs(rho_km[-1]  - rho_km[0])/float(len(rho_km) - 1)
+    p = 1
+    q = int(round(1.0/(ts_avg*freq)))
+    dr_grid = float(p)/(q*freq)
+
+    # Uniform radius grid at near-raw resolution to interpolate to
+    # For ingress, this implicitly reverses radius scale!
+    n_pts = round(abs(rho_km[-1] - rho_km[0])/dr_grid)
+    rho_grid = min(rho_km) + dr_grid*np.arange(n_pts)
+
+    # Interpolate to near-raw resolution
+    # For ingress, this implicitly reverses radius scale!
+    vec_grid_interp = interp1d(rho_km, vec, kind='linear',
+                               fill_value='extrapolate')
+    vec_grid = vec_grid_interp(rho_grid)
+
+    return rho_grid, vec_grid, p, q
+
+
+def resample_IQ(rho_km, IQ_c, dr_desired, dr_km_tol=0.01, TEST=False):
+    """Function to resample I and Q to uniformly spaced radius. Based off of
+    Matlab's resample function
+
+    Example:
+        >>> fit_inst = FreqOffsetFit(RSR_inst, geo_inst, f_spm, f_offset, \
+                f_USO, kernels, k=k, rho_exclude=rho_exclude)
+        >>> (spm_raw, IQ_c_raw) = fit_inst.get_IQ_c()
+        >>> (rho_km_desired, IQ_c_resampled) = resample_IQ(rho_km_raw, \
+                IQ_c_raw, dr_desired, dr_km_tol=dr_km_tol, TEST=TEST)
+
+    Args:
+        rho_km (np.ndarray): Set of ring intercept point values at initial
+            resolution before resampling
+        IQ_c (np.ndarray): Frequency corrected complex signal at initial
+            resolution before resampling
+        dr_desired (float): Desired final radial spacing to resample to
+        dr_km_tol (float): Maximum distance from an integer number of
+            dr_desired that the first rho value will be at. Makes the final
+            set of rho values more regular-looking. For example, if you say
+            dr_km_tol=0.01 with a dr_desired of 0.25, your final set of rho
+            values might look something like [70000.26, 70000.51, ...]
+        TEST (bool): Testing variable to print out the first few resampled
+            results
+    """
+
+    # Begin raw resolution rho within dr_km_tol of integer number of dr_desired
+    rho_km_remainder = rho_km % dr_desired
+    rho_km_start_ind = int((np.argwhere(rho_km_remainder < dr_km_tol))[0])
+    rho_km = rho_km[rho_km_start_ind:-1]
+    IQ_c = IQ_c[rho_km_start_ind:-1]
+
+    I_c = np.real(IQ_c)
+    Q_c = np.imag(IQ_c)
+
+    # Pre-resampling steps. Interpolates to uniform radius at near-raw spacing
+    (rho_km_uniform, I_c_uniform, p, q) = pre_resample(rho_km, I_c,
+                                                       1.0/dr_desired)
+    (rho_km_uniform, Q_c_uniform, p, q) = pre_resample(rho_km, Q_c,
+                                                       1.0/dr_desired)
+
+    # Downsample by factor q to desired final spacing
+    I_c_desired = signal.resample_poly(I_c_uniform, p, q)
+    Q_c_desired = signal.resample_poly(Q_c_uniform, p, q)
+
+    rho_km_desired = rho_km_uniform[0] + dr_desired*np.arange(len(I_c_desired))
+
+    if TEST:
+        print('First 10 rho, I_c, Q_c:')
+        for i in range(10):
+            print('%24.16f %24.16f %24.16f' %
+                  (rho_km_desired[i], I_c_desired[i], Q_c_desired[i]))
+
+    return rho_km_desired, I_c_desired + 1j*Q_c_desired
+
+
+if __name__ == '__main__':
+    pass