Source code for ocbpy.ocb_scaling

#!/usr/bin/env python# -*- coding: utf-8 -*-
# Copyright (C) 2017, AGB & GC
# Full license can be found in License.md
#-----------------------------------------------------------------------------
"""Scale data affected by magnetic field direction or electric field

Routines
-------------------------------------------------------------------------------
normal_evar(evar, unscaled_r, scaled_r)
    Normalise a variable proportaional to the electric field (such as velocity)
normal_curl_evar(curl_evar, unscaled_r, scaled_r)
    Normalise a variable proportional to the curl of the electric field (such
    as vorticity)

Classes
-------------------------------------------------------------------------------
VectorData(object)
    Holds vector data in AACGM N-E-Z coordinates along with location
    information.  Converts vector from AACGM to OCB coordinates.

Moduleauthor
-------------------------------------------------------------------------------
Angeline G. Burrell (AGB), 12 May 2017, University of Texas, Dallas (UTDallas)

References
-------------------------------------------------------------------------------
Chisham, G. (2017), A new methodology for the development of high-latitude
 ionospheric climatologies and empirical models, Journal of Geophysical
 Research: Space Physics, 122, doi:10.1002/2016JA023235.
"""
import logbook as logging
import numpy as np

[docs]def normal_evar(evar, unscaled_r, scaled_r): """ Normalise a variable proportional to the electric field Parameters ----------- evar : (float) Variable related to electric field (e.g. velocity) unscaled_r : (float) Radius of polar cap in degrees scaled_r : (float) Radius of normalised OCB polar cap in degrees Returns -------- nvar : (float) Normalised variable Notes ------- Assumes that the cross polar cap potential is fixed across the polar cap regardless of the radius of the Open Closed field line Boundary. This is commonly assumed when looking at statistical patterns that control the IMF (which accounts for dayside reconnection) and assume that the nightside reconnection influence is averaged out over the averaged period. References ----------- Chisham, G. (2017), A new methodology for the development of high‐latitude ionospheric climatologies and empirical models, Journal of Geophysical Research: Space Physics, doi:10.1002/2016JA023235. """ nvar = evar * unscaled_r / scaled_r return nvar
[docs]def normal_curl_evar(curl_evar, unscaled_r, scaled_r): """ Normalise a variable proportional to the curl of the electric field Parameters ----------- curl_evar : (float) Variable related to electric field (e.g. vorticity) unscaled_r : (float) Radius of polar cap in degrees scaled_r : (float) Radius of normalised OCB polar cap in degrees Returns -------- nvar : (float) Normalised variable Notes ------- Assumes that the cross polar cap potential is fixed across the polar cap regardless of the radius of the Open Closed field line Boundary. This is commonly assumed when looking at statistical patterns that control the IMF (which accounts for dayside reconnection) and assume that the nightside reconnection influence is averaged out over the averaged period. References ----------- Chisham, G. (2017), A new methodology for the development of high‐latitude ionospheric climatologies and empirical models, Journal of Geophysical Research: Space Physics, doi:10.1002/2016JA023235. """ nvar = curl_evar * (unscaled_r / scaled_r)**2 return nvar
[docs]class VectorData(object): """ Object containing a vector data point Parameters ----------- dat_ind : (int) Data index (zero offset) ocb_ind : (int) OCBoundary record index matched to this data index (zero offset) aacgm_lat : (float) Vector AACGM latitude (degrees) aacgm_mlt : (float) Vector AACGM MLT (hours) ocb_lat : (float) Vector OCB latitude (degrees) (default=np.nan) ocb_mlt : (float) Vector OCB MLT (hours) (default=np.nan) aacgm_n : (float) AACGM North pointing vector (positive towards North) (default=0.0) aacgm_e : (float) AACGM East pointing vector (completes right-handed coordinate system (default = 0.0) aacgm_z : (float) AACGM Vertical pointing vector (positive down) (default=0.0) aacgm_mag : (float) Vector magnitude (default=np.nan) scale_func : (function) Function for scaling AACGM magnitude with arguements: [measurement value, mesurement AACGM latitude (degrees), mesurement OCB latitude (degrees)] (default=None) dat_name : (str) Data name (default=None) dat_units : (str) Data units (default=None) Attributes ------------ dat_name : (str or NoneType) Name of data dat_units : (str or NoneType) Units of data dat_ind : (int) Vector data index in external data array ocb_ind : (int) OCBoundary rec_ind value that matches dat_ind unscaled_r : (float) Radius of polar cap in degrees scaled_r : (float) Radius of normalised OCB polar cap in degrees aacgm_n : (float) AACGM north component of data vector (default=0.0) aacgm_e : (float) AACGM east component of data vector (default=0.0) aacgm_z : (float) AACGM vertical component of data vector (default=0.0) aacgm_mag : (float) Magnitude of data vector in AACGM coordinates (default=np.nan) aacgm_lat : (float) AACGM latitude of data vector in degrees aacgm_mlt : (float) AACGM MLT of data vector in hours ocb_n : (float) OCB north component of data vector (default=np.nan) ocb_e : (float) OCB east component of data vector (default=np.nan) ocb_z : (float) OCB vertical component of data vector (default=np.nan) ocb_mag : (float) OCB magnitude of data vector (default=np.nan) ocb_lat : (float) OCB latitude of data vector in degrees (default=np.nan) ocb_mlt : (float) OCB MLT of data vector in hours (default=np.nan) ocb_quad : (int) AACGM quadrant of OCB pole (default=0) vec_quad : (int) AACGM quadrant of Vector (default=0) pole_angle : (float) Angle at vector location appended by AACGM and OCB poles in degrees (default=np.nan) aacgm_naz : (float) AACGM north azimuth of data vector in degrees (default=np.nan) ocb_aacgm_lat : (float) AACGM latitude of OCB pole in degrees (default=np.nan) ocb_aacgm_mlt : (float) AACGM MLT of OCB pole in hours (default=np.nan) scale_func : (function or NoneType) Funciton that scales the magnitude of the data vector from AACGM polar cap coverage to OCB polar cap coverage Methods ----------- set_ocb(ocb, scale_func=None) Set the ocb coordinates and vector values define_quadrants() Define the OCB pole and vector AACGM quadrant scale_vector() Scale the data vector into OCB coordinates calc_ocb_polar_angle() calculate the OCB north azimuth angle calc_ocb_vec_sign(north=False, east=False, quads=dict()) calculate the signs of the OCB scaled vector components calc_vec_pole_angle(angular_res=1.0e-3) calculate the vector angle of the vector-poles triangle """ def __init__(self, dat_ind, ocb_ind, aacgm_lat, aacgm_mlt, ocb_lat=np.nan, ocb_mlt=np.nan, aacgm_n=0.0, aacgm_e=0.0, aacgm_z=0.0, aacgm_mag=np.nan, dat_name=None, dat_units=None, scale_func=None): """ Initialize VectorData object Parameters ----------- dat_ind : (int) Data index (zero offset) ocb_ind : (int) OCBoundary record index matched to this data index (zero offset) aacgm_lat : (float) Vector AACGM latitude (degrees) aacgm_mlt : (float) Vector AACGM MLT (hours) ocb_lat : (float) Vector OCB latitude (degrees) (default=np.nan) ocb_mlt : (float) Vector OCB MLT (hours) (default=np.nan) aacgm_n : (float) AACGM North pointing vector (positive towards North) (default=0.0) aacgm_e : (float) AACGM East pointing vector (completes right-handed coordinate system (default = 0.0) aacgm_z : (float) AACGM Vertical pointing vector (positive down) (default=0.0) aacgm_mag : (float) Vector magnitude (default = np.nan) dat_name : (str) Data name (default=None) dat_units : (str) Data units (default=None) scale_func : (function) Function for scaling AACGM magnitude with arguements: [measurement value, mesurement AACGM latitude (degrees), mesurement OCB latitude (degrees)] Not necessary if no magnitude scaling is needed. (default=None) Returns -------- self : Initialised VectorData class object by setting AACGM values """ # Assign the vector data name and units self.dat_name = dat_name self.dat_units = dat_units # Assign the data and OCB indices self.dat_ind = dat_ind self.ocb_ind = ocb_ind # Assign the AACGM vector values and location self.aacgm_n = aacgm_n self.aacgm_e = aacgm_e self.aacgm_z = aacgm_z self.aacgm_lat = aacgm_lat self.aacgm_mlt = aacgm_mlt if np.isnan(aacgm_mag): self.aacgm_mag = np.sqrt(aacgm_n**2 + aacgm_e**2 + aacgm_z**2) else: self.aacgm_mag = aacgm_mag # Assign the OCB vector default values and location self.ocb_n = np.nan self.ocb_e = np.nan self.ocb_z = np.nan self.ocb_mag = np.nan self.ocb_lat = ocb_lat self.ocb_mlt = ocb_mlt # Assign the default pole locations, relative angles, and quadrants self.ocb_quad = 0 self.vec_quad = 0 self.pole_angle = np.nan self.aacgm_naz = np.nan self.ocb_aacgm_lat = np.nan self.ocb_aacgm_mlt = np.nan # Assign the vector scaling function self.scale_func = scale_func return def __repr__(self): """ Provide readable representation of the DataVector object """ out = "Vector data:" if self.dat_name is not None: out = "{:s} {:s}".format(out, self.dat_name) if self.dat_units is not None: out = "{:s} ({:s})".format(out, self.dat_units) out = "{:s}\nData Index {:d}\tOCB Index {:d}\n".format(out, self.dat_ind, self.ocb_ind) out = "{:s}-------------------------------------------\n".format(out) out = "{:s}Location: [Mag. Lat. (degrees), MLT (hours)]\n".format(out) out = "{:s} AACGM: [{:.3f}, {:.3f}]\n".format(out, self.aacgm_lat, self.aacgm_mlt) if not np.isnan(self.ocb_lat) and not np.isnan(self.ocb_mlt): out = "{:s} OCB: [{:.3f}, {:.3f}]\n".format(out, self.ocb_lat, self.ocb_mlt) out = "\n{:s}-------------------------------------------\n".format(out) out = "{:s}Value: Magnitude [N, E, Z]\n".format(out) out = "{:s}AACGM: {:.3g} [".format(out, self.aacgm_mag) out = "{:s}{:.3g}, {:.3g}, {:.3g}]\n".format(out, self.aacgm_n, self.aacgm_e, self.aacgm_z) if not np.isnan(self.ocb_mag): out = "{:s} OCB: {:.3g} [".format(out, self.ocb_mag) out = "{:s}{:.3g}, {:.3g}, {:.3g}]\n".format(out, self.ocb_n, self.ocb_e, self.ocb_z) out = "\n{:s}-------------------------------------------\n".format(out) if self.scale_func is None: out = "{:s}No magnitude scaling function provided\n".format(out) else: out = "{:s}Scaling function: ".format(out) out = "{:s}{:s}\n".format(out, self.scale_func.__name__) return out def __str__(self): """ Provide readable representation of the DataVector object """ out = self.__repr__() return out
[docs] def set_ocb(self, ocb, scale_func=None): """ Set the OCBoundary values for this data point Parameters ------------ ocb : (OCBoundary) Open Closed Boundary class object scale_func : (function) Function for scaling AACGM magnitude with arguements: [measurement value, mesurement AACGM latitude (degrees), mesurement OCB latitude (degrees)] Not necessary if defined earlier or no scaling is needed. (default=None) Updates --------- self.unscaled_r : (float) Radius of polar cap in degrees self.scaled_r : (float) Radius of normalised OCB polar cap in degrees self.ocb_n : (float) Vector OCB North component self.ocb_e : (float) Vector OCB East component self.ocb_z : (float) Vector OCB vertical component (positive downward) self.ocb_mag : (float) Vector OCB magnitude self.ocb_lat : (float) Vector OCB latitude, if not updated already (degrees) self.ocb_mlt : (float) Vector OCB MLT, if not updated already (hours) self.ocb_quad : (int) OCB pole AACGM quadrant self.vec_quad : (int) Vector AACGM quadrant self.pole_angle : (float) Angle at the vector in the triangle formed by the poles and vector (degrees) self.aacgm_naz : (float) AACGM north azimuth angle (degrees) self.ocb_aacgm_lat : (float) AACGM latitude of the OCB pole (degrees) self.ocb_aacgm_mlt : (float) AACGM MLT of the OCB pole (hours) self.scale_func : (function) Function for scaling AACGM magnitude with arguements: [measurement value, unscaled polar cap radius (degrees), scaled polar cap radius (degrees)] Not necessary if defined earlier or if no scaling is needed. """ # Set the AACGM coordinates of the OCB pole self.unscaled_r = ocb.r[self.ocb_ind] self.scaled_r = 90.0 - abs(ocb.boundary_lat) self.ocb_aacgm_mlt = ocb.phi_cent[self.ocb_ind] / 15.0 self.ocb_aacgm_lat = 90.0 - ocb.r_cent[self.ocb_ind] # If the OCB vector coordinates weren't included in the initial info, # update them here if np.isnan(self.ocb_lat) or np.isnan(self.ocb_mlt): self.ocb_lat, self.ocb_mlt = ocb.normal_coord(self.aacgm_lat, self.aacgm_mlt) # Get the angle at the data vector appended by the AACGM and OCB poles self.calc_vec_pole_angle() # Set the OCB and Vector quadrants if not np.isnan(self.pole_angle): self.define_quadrants() # Set the scaling function if self.scale_func is None: if scale_func is None: # This is not necessarily a bad thing, if the value does not # need to be scaled. logging.info("no scaling function provided") else: self.scale_func = scale_func # Assign the OCB vector default values and location. Will also # update the AACGM north azimuth of the vector. self.scale_vector() return
[docs] def define_quadrants(self): """ Determine which quadrants (in AACGM coordinates) the OCB pole and data vector lie in Requires --------- self.ocb_aacgm_mlt : (float) OCB pole MLT in AACGM coordinates in hours self.aacgm_mlt : (float) Vector AACGM MLT in hours self.pole_angle : (float) vector angle in poles-vector triangle in degrees Updates -------- self.ocb_quad : (int) OCB pole quadrant self.vec_quad : (int) Vector quadrant Notes ------ North (N) and East (E) are defined by the AACGM directions centred on the data vector location, assuming vertical is positive downwards Quadrants: 1 [N, E]; 2 [N, W]; 3 [S, W]; 4 [S, E] """ # Test input assert(not np.isnan(self.ocb_aacgm_mlt)), \ logging.error("OCB pole location required") assert(not np.isnan(self.aacgm_mlt)), \ logging.error("Vector AACGM location required") assert(not np.isnan(self.pole_angle)), \ logging.error("vector angle in poles-vector triangle required") # Determine where the OCB pole is relative to the data vector ocb_adj_mlt = self.ocb_aacgm_mlt - self.aacgm_mlt while ocb_adj_mlt < 0.0: ocb_adj_mlt += 24.0 if abs(ocb_adj_mlt) >= 24.0: ocb_adj_mlt -= 24.0 * np.sign(ocb_adj_mlt) if self.pole_angle < 90.0: # OCB pole lies in quadrant 1 or 2 self.ocb_quad = 1 if ocb_adj_mlt < 12.0 else 2 else: # OCB poles lies in quadrant 3 or 4 self.ocb_quad = 3 if ocb_adj_mlt < 24.0 else 4 # Now determine which quadrant the vector is pointed into if self.aacgm_n >= 0.0: self.vec_quad = 1 if self.aacgm_e >= 0.0 else 2 else: self.vec_quad = 4 if self.aacgm_e >= 0.0 else 3 return
[docs] def scale_vector(self): """ Normalise a variable proportional to the curl of the electric field. Requires --------- self.ocb_lat : (float) OCB latitude in degrees self.ocb_mlt : (float) OCB MLT in hours self.ocb_aacgm_mlt : (float) OCB pole MLT in AACGM coordinates in hours self.pole_angle : (float) vector angle in poles-vector triangle Updates -------- ocb_n : (float) OCB scaled north component ocb_e : (float) OCB scaled east component ocb_z : (float) OCB scaled vertical component ocb_mag : (float) OCB scaled magnitude """ # Test input assert(not np.isnan(self.ocb_lat) and not np.isnan(self.ocb_mlt)), \ logging.error("OCB coordinates required") assert(not np.isnan(self.ocb_aacgm_mlt)), \ logging.error("OCB pole location required") assert(not np.isnan(self.pole_angle)), \ logging.error("vector angle in poles-vector triangle required") # Scale vertical component if self.scale_func is None or self.aacgm_z != 0.0: self.ocb_z = self.aacgm_z else: self.ocb_z = self.scale_func(self.aacgm_z, self.unscaled_r, self.scaled_r) if self.aacgm_n == 0.0 and self.aacgm_e == 0.0: # There's no magnitude, so nothing to adjust self.ocb_n = 0.0 self.ocb_e = 0.0 elif self.pole_angle == 0.0 or self.pole_angle == 180.0: # The measurement is aligned with the AACGM and OCB poles if self.scale_func is None: self.ocb_n = self.aacgm_n self.ocb_e = self.aacgm_e else: self.ocb_n = self.scale_func(self.aacgm_n, self.unscaled_r, self.scaled_r) self.ocb_e = self.scale_func(self.aacgm_e, self.unscaled_r, self.scaled_r) if self.pole_angle == 0.0 and self.aacgm_lat >= self.ocb_aacgm_lat: # The measurement is on or between the poles self.ocb_n *= -1.0 self.ocb_e *= -1.0 else: # If not defined, get the OCB and vector quadrants if(self.ocb_quad == 0 or self.vec_quad == 0): self.define_quadrants() if(self.ocb_quad == 0 or self.vec_quad == 0): logging.error("unable to define OCB and vector quadrants") return # Get the unscaled 2D vector magnitude vmag = np.sqrt(self.aacgm_n**2 + self.aacgm_e**2) # Calculate the AACGM north azimuth in degrees self.aacgm_naz = np.degrees(np.arccos(self.aacgm_n / vmag)) # Get the OCB north azimuth in radians ocb_angle = np.radians(self.calc_ocb_polar_angle()) # Get the sign of the North and East components vsigns = self.calc_ocb_vec_sign(north=True, east=True) # Scale the vector along the OCB north and account for # any changes associated with adjusting the size of the polar cap if self.scale_func is not None: vmag = self.scale_func(vmag, self.unscaled_r, self.scaled_r) self.ocb_n = vsigns['north'] * vmag * np.cos(ocb_angle) self.ocb_e = vsigns['east'] * vmag * np.sin(ocb_angle) # Calculate the scaled OCB vector magnitude self.ocb_mag = np.sqrt(self.ocb_n**2 + self.ocb_e**2 + self.ocb_z**2) return
[docs] def calc_ocb_polar_angle(self): """ Calculate the OCB north azimuth angle Requires --------- self.ocb_quad : (int) OCB quadrant self.vec_quad : (int) Vector quadrant self.aacgm_naz : (float) AACGM polar angle self.pole_angle : (float) Vector angle between AACGM pole, vector origin, and OCB pole Returns -------- ocb_naz : (float) Angle between measurement vector and OCB pole in degrees """ quad_range = np.arange(1, 5) assert self.ocb_quad in quad_range, \ logging.error("OCB quadrant undefined") assert self.vec_quad in quad_range, \ logging.error("Vector quadrant undefined") assert not np.isnan(self.aacgm_naz), \ logging.error("AACGM polar angle undefined") assert not np.isnan(self.pole_angle), \ logging.error("Vector angle undefined") # Initialise the output and set the quadrant dictionary ocb_naz = np.nan quads = {o:{v:True if self.ocb_quad == o and self.vec_quad == v else False for v in quad_range} for o in quad_range} # Calculate OCB polar angle based on quadrants and other angles if(((quads[2][4] or quads[2][2]) and self.aacgm_naz > self.pole_angle) or (self.aacgm_naz > self.pole_angle and quads[1][1]) or (quads[1][4] and self.aacgm_naz <= self.pole_angle + 90.0)): ocb_naz = self.aacgm_naz - self.pole_angle elif((self.aacgm_naz <= self.pole_angle and quads[1][1]) or (self.aacgm_naz <= self.pole_angle and (quads[2][4] or quads[2][2])) or (self.aacgm_naz > self.pole_angle - 90.0 and (quads[4][1] or quads[4][3] or quads[3][4] or quads[3][2]))): ocb_naz = self.pole_angle - self.aacgm_naz elif(self.aacgm_naz <= 90.0 - self.pole_angle and (quads[1][2] or quads[2][1] or quads[2][3])): ocb_naz = self.aacgm_naz + self.pole_angle elif((self.aacgm_naz > 90.0 - self.pole_angle and (quads[1][2] or quads[2][1] or quads[2][3])) or ((quads[4][4] or quads[4][2] or quads[3][1] or quads[3][3] or quads[1][3]) and self.aacgm_naz <= 180.0 - self.pole_angle)): ocb_naz = 180.0 - self.aacgm_naz - self.pole_angle elif(((quads[3][1] or quads[3][3] or quads[4][4] or quads [4][2] or quads[1][3]) and self.aacgm_naz > 180.0 - self.pole_angle) or (quads[1][4] and self.aacgm_naz > self.pole_angle + 90.0)): ocb_naz = self.aacgm_naz - 180.0 + self.pole_angle elif(self.aacgm_naz <= self.pole_angle - 90.0 and (quads[3][4] or quads[3][2] or quads[4][1] or quads[4][3])): ocb_naz = 180.0 - self.pole_angle + self.aacgm_naz return ocb_naz
[docs] def calc_ocb_vec_sign(self, north=False, east=False, quads=dict()): """ Get the sign of the North and East components Parameters ------------ north : (boolian) Get the sign of the north component (default=False) east : (boolian) Get the sign of the east component (default=False) quads : (dictionary) Dictionary of boolian values for OCB and Vector quadrants (default=dict()) Requires ---------- self.ocb_quad : (int) OCB pole quadrant self.vec_quad : (int) Vector quadrant self.aacgm_naz : (float) AACGM polar angle in degrees self.pole_angle : (float) Vector angle in degrees Returns --------- vsigns : (dict) Dictionary with keys 'north' and 'east' containing the desired signs """ quad_range = np.arange(1, 5) # Test input assert north or east, logging.warning("must set at least one direction") assert self.ocb_quad in quad_range, \ logging.error("OCB quadrant undefined") assert self.vec_quad in quad_range, \ logging.error("Vector quadrant undefined") assert not np.isnan(self.aacgm_naz), \ logging.error("AACGM polar angle undefined") assert not np.isnan(self.pole_angle), \ logging.error("Vector angle undefined") # Initialise output vsigns = {"north":0, "east":0} # If necessary, initialise quadrant dictionary if not np.all(quads.keys() == quad_range): quads = {o:{v:True if self.ocb_quad == o and self.vec_quad == v else False for v in quad_range} for o in quad_range} if north: pole_minus = self.pole_angle - 90.0 minus_pole = 90.0 - self.pole_angle pole_plus = self.pole_angle + 90.0 if(quads[1][1] or quads[2][2] or quads[3][3] or quads[4][4] or (quads[1][4] and self.aacgm_naz <= pole_plus) or (quads[1][2] and self.aacgm_naz > minus_pole) or (quads[2][1] and self.aacgm_naz <= minus_pole) or ((quads[3][4] or quads[4][3]) and self.aacgm_naz <= pole_minus) or ((quads[3][2] or quads[4][1]) and self.aacgm_naz > pole_minus) or (quads[2][3] and self.aacgm_naz > minus_pole)): vsigns["north"] = 1 elif((quads[1][2] and self.aacgm_naz > minus_pole) or (quads[1][4] and self.aacgm_naz > pole_plus) or (quads[2][1] and self.aacgm_naz > minus_pole) or ((quads[4][1] or quads[3][2]) and self.aacgm_naz <= pole_minus) or (quads[2][3] and self.aacgm_naz <= minus_pole) or ((quads[4][3] or quads[3][4]) and self.aacgm_naz > pole_minus) or quads[1][3] or quads[2][4] or quads[3][1] or quads[4][2]): vsigns["north"] = -1 if east: minus_pole = 180.0 - self.pole_angle if(quads[1][4] or quads[2][1] or quads[3][2] or quads[4][3] or (quads[1][1] and self.aacgm_naz > self.pole_angle) or (quads[1][3] and self.aacgm_naz > minus_pole) or ((quads[4][4] or quads[3][1]) and self.aacgm_naz <= minus_pole) or (quads[2][4] and self.aacgm_naz > self.pole_angle) or ((quads[4][2] or quads[3][3]) and self.aacgm_naz > minus_pole) or (quads[2][2] and self.aacgm_naz <= self.pole_angle)): vsigns["east"] = 1 elif(quads[1][2] or quads[2][3] or quads[3][4] or quads[4][1] or ((quads[4][4] or quads[3][1]) and self.aacgm_naz > minus_pole) or (quads[2][2] and self.aacgm_naz > self.pole_angle) or ((quads[4][2] or quads[3][3])and self.aacgm_naz <= minus_pole) or (quads[1][3] and self.aacgm_naz <= minus_pole) or ((quads[1][1] or quads[2][4]) and self.aacgm_naz <= self.pole_angle)): vsigns["east"] = -1 return vsigns
[docs] def calc_vec_pole_angle(self): """calculates the angle between the AACGM pole, a measurement, and the OCB pole using spherical triginometry Requires --------- self.aacgm_mlt : (float) AACGM MLT of vector origin in hours self.aacgm_lat : (float) AACGM latitude of vector origin in degrees self.ocb_aacgm_mlt : (float) AACGM MLT of the OCB pole in hours self.ocb_aacgm_lat : (float) AACGM latitude of the OCB pole in degrees Updates -------- self.pole_angle : (float) Angle in degrees between AACGM north, a measurement, and OCB north """ # Test input assert(not np.isnan(self.aacgm_mlt)), \ logging.error("AACGM MLT of Vector undefinded") assert(not np.isnan(self.ocb_aacgm_mlt)), \ logging.error("AACGM MLT of OCB pole is undefined") assert(not np.isnan(self.ocb_aacgm_lat)), \ logging.error("AACGM latitude of OCB pole is undefined") assert(not np.isnan(self.aacgm_lat)), \ logging.error("AACGM latitude of Vector undefined") # Convert the AACGM MLT of the observation and OCB pole to radians, # then calculate the difference between them. del_long = (self.ocb_aacgm_mlt - self.aacgm_mlt) * np.pi / 12.0 if del_long < 0.0: del_long += 2.0 * np.pi if del_long == 0.0: self.pole_angle = 0.0 return if del_long == np.pi: self.pole_angle = 180.0 return # Find the distance in radians between the two poles hemisphere = np.sign(self.ocb_aacgm_lat) rad_pole = hemisphere * np.pi * 0.5 del_pole = hemisphere * (rad_pole - np.radians(self.ocb_aacgm_lat)) # Get the distance in radians between the AACGM pole and the data point del_vect = hemisphere * (rad_pole - np.radians(self.aacgm_lat)) # Use the law of haversines, which goes to the spherical trigonometric # cosine rule for sides at large angles, but is more robust at small # angles, to find the length of the last side of the spherical triangle. del_ocb = archav(hav(del_pole - del_vect) + np.sin(del_pole) * np.sin(del_vect) * hav(del_long)) # Again use law of haversines, this time to find the polar angle hav_pole_angle = (hav(del_pole) - hav(del_vect - del_ocb)) \ / (np.sin(del_vect) * np.sin(del_ocb)) self.pole_angle = np.degrees(archav(hav_pole_angle)) return
[docs]def hav(alpha): """ Formula for haversine Parameters ---------- alpha : (float) Angle in radians Returns -------- hav_alpha : (float) Haversine of alpha, equal to the square of the sine of half-alpha """ hav_alpha = np.sin(alpha * 0.5)**2 return hav_alpha
[docs]def archav(hav): """ Formula for the inverse haversine Parameters ----------- hav : (float) Haversine of an angle Returns --------- alpha : (float) Angle in radians """ alpha = 2.0 * np.arcsin(np.sqrt(hav)) return alpha