aacgmv2.deprecated.igrf_dipole_axis() - Code Metrics - Inspection of "V2.6 update" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — develop (#47)

by Angeline

created 2019-12-30 22:35 UTC

aacgmv2.deprecated.igrf_dipole_axis() A

↳ Parent: aacgmv2.deprecated

Complexity

Conditions

Size

Total Lines	70
Code Lines	30

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	30
dl	0
loc	70
rs	9.16
c	0
b	0
f	0
cc	3
nop	1

How to fix Long Method

# Copyright (C) 2019 NRL 
# Author: Angeline Burrell
# Disclaimer: This code is under the MIT license, whose details can be found at
# the root in the LICENSE file
#
# -*- coding: utf-8 -*-
"""Pythonic wrappers for AACGM-V2 C functions that were depricated in the
change from version 2.0.0 to version 2.0.2

References
-------------------------------------------------------------------------------
Laundal, K. M. and A. D. Richmond (2016), Magnetic Coordinate Systems, Space
 Sci. Rev., doi:10.1007/s11214-016-0275-y.

"""

from __future__ import division, absolute_import, unicode_literals
import datetime as dt
import numpy as np
import warnings
import aacgmv2

dep_str = "".join(["Deprecated routine will be removed in version 2.6.1 ",
                   "unless users express interest in keeping it"])

def subsol(year, doy, utime):
    """Finds subsolar geocentric longitude and latitude.

    Parameters
    ------------
    year : (int)
        Calendar year between 1601 and 2100
    doy : (int)
        Day of year between 1-365/366
    utime : (float)
        Seconds since midnight on the specified day

    Returns
    ---------
    sbsllon : (float)
        Subsolar longitude in degrees E for the given date/time
    sbsllat : (float)
        Subsolar latitude in degrees N for the given date/time

    Notes
    --------
    Based on formulas in Astronomical Almanac for the year 1996, p. C24.
    (U.S. Government Printing Office, 1994). Usable for years 1601-2100,
    inclusive. According to the Almanac, results are good to at least 0.01
    degree latitude and 0.025 degrees longitude between years 1950 and 2050.
    Accuracy for other years has not been tested. Every day is assumed to have
    exactly 86400 seconds; thus leap seconds that sometimes occur on December
    31 are ignored (their effect is below the accuracy threshold of the
    algorithm).
    After Fortran code by A. D. Richmond, NCAR. Translated from IDL
    by K. Laundal.

    """
    warnings.warn(dep_str, category=FutureWarning)

    # Convert from 4 digit year to 2 digit year
    yr2 = year - 2000

    if year >= 2101:
        aacgmv2.logger.error('subsol invalid after 2100. Input year is:', year)

    # Determine if this year is a leap year
    nleap = np.floor((year - 1601) / 4)
    nleap = nleap - 99
    if year <= 1900:
        if year <= 1600:
            print('subsol.py: subsol invalid before 1601. Input year is:', year)
        ncent = np.floor((year - 1601) / 100)
        ncent = 3 - ncent
        nleap = nleap + ncent

    # Calculate some of the coefficients needed to deterimine the mean longitude
    # of the sun and the mean anomaly
    l_0 = -79.549 + (-0.238699 * (yr2 - 4 * nleap) + 3.08514e-2 * nleap)
    g_0 = -2.472 + (-0.2558905 * (yr2 - 4 * nleap) - 3.79617e-2 * nleap)

    # Days (including fraction) since 12 UT on January 1 of IYR2:
    dfrac = (utime / 86400 - 1.5) + doy

    # Mean longitude of Sun:
    l_sun = l_0 + 0.9856474 * dfrac

    # Mean anomaly:
    grad = np.radians(g_0 + 0.9856003 * dfrac)

    # Ecliptic longitude:
    lmrad = np.radians(l_sun + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
    sinlm = np.sin(lmrad)

    # Days (including fraction) since 12 UT on January 1 of 2000:
    epoch_day = dfrac + 365.0 * yr2 + nleap

    # Obliquity of ecliptic:
    epsrad = np.radians(23.439 - 4.0e-7 * epoch_day)

    # Right ascension:
    alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))

    # Declination, which is the subsolar latitude:
    sbsllat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))

    # Equation of time (degrees):
    etdeg = l_sun - alpha
    etdeg = etdeg - 360.0 * np.round(etdeg / 360.0)

    # Apparent time (degrees):
    aptime = utime / 240.0 + etdeg    # Earth rotates one degree every 240 s.

    # Subsolar longitude:
    sbsllon = 180.0 - aptime
    sbsllon = sbsllon - 360.0 * np.round(sbsllon / 360.0)

    return sbsllon, sbsllat

def gc2gd_lat(gc_lat):
    """Convert geocentric latitude to geodetic latitude using WGS84.

    Parameters
    -----------
    gc_lat : (array_like or float)
        Geocentric latitude in degrees N

    Returns
    ---------
    gd_lat : (same as input)
        Geodetic latitude in degrees N
    """
    warnings.warn(dep_str, category=FutureWarning)

    
    wgs84_e2 = 0.006694379990141317 - 1.0
    return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / wgs84_e2))

def igrf_dipole_axis(date):
    """Get Cartesian unit vector pointing at dipole pole in the north,
    according to IGRF

    Parameters
    -------------
    date : (dt.datetime)
        Date and time

    Returns
    ----------
    m_0: (np.ndarray)
        Cartesian 3 element unit vector pointing at dipole pole in the north
        (geocentric coords)

    Notes
    ----------
    IGRF coefficients are read from the igrf12coeffs.txt file. It should also
    work after IGRF updates.  The dipole coefficients are interpolated to the
    date, or extrapolated if date > latest IGRF model
    """
    warnings.warn(dep_str, category=FutureWarning)

    # get time in years, as float:
    year = date.year
    doy = date.timetuple().tm_yday
    year_days = int(dt.date(date.year, 12, 31).strftime("%j"))
    year = year + doy / year_days

    # read the IGRF coefficients
    with open(aacgmv2.IGRF_COEFFS, 'r') as f_igrf:
        lines = f_igrf.readlines()

    years = lines[3].split()[3:][:-1]
    years = np.array(years, dtype=float)  # time array

    g10 = lines[4].split()[3:]
    g11 = lines[5].split()[3:]
    h11 = lines[6].split()[3:]

    # secular variation coefficients (for extrapolation)
    g10sv = np.float32(g10[-1])
    g11sv = np.float32(g11[-1])
    h11sv = np.float32(h11[-1])

    # model coefficients:
    g10 = np.array(g10[:-1], dtype=float)
    g11 = np.array(g11[:-1], dtype=float)
    h11 = np.array(h11[:-1], dtype=float)

    # get the gauss coefficient at given time:
    if year <= years[-1]:
        # regular interpolation
        g10 = np.interp(year, years, g10)
        g11 = np.interp(year, years, g11)
        h11 = np.interp(year, years, h11)
    else:
        # extrapolation
        dyear = year - years[-1]
        g10 = g10[-1] + g10sv * dyear
        g11 = g11[-1] + g11sv * dyear
        h11 = h11[-1] + h11sv * dyear

    # calculate pole position
    B_0 = np.sqrt(g10**2 + g11**2 + h11**2)

    # Calculate output
    m_0 = -np.array([g11, h11, g10]) / B_0

    return m_0


1			# Copyright (C) 2019 NRL
2			# Author: Angeline Burrell
3			# Disclaimer: This code is under the MIT license, whose details can be found at
4			# the root in the LICENSE file
5			#
6			# -- coding: utf-8 --
7			"""Pythonic wrappers for AACGM-V2 C functions that were depricated in the
8			change from version 2.0.0 to version 2.0.2
9
10			References
11			-------------------------------------------------------------------------------
12			Laundal, K. M. and A. D. Richmond (2016), Magnetic Coordinate Systems, Space
13			Sci. Rev., doi:10.1007/s11214-016-0275-y.
14
15			"""
16
17			from __future__ import division, absolute_import, unicode_literals
18			import datetime as dt
19			import numpy as np
20			import warnings
21			import aacgmv2
22
23			dep_str = "".join(["Deprecated routine will be removed in version 2.6.1 ",
24			"unless users express interest in keeping it"])
25
26			def subsol(year, doy, utime):
27			"""Finds subsolar geocentric longitude and latitude.
28
29			Parameters
30			------------
31			year : (int)
32			Calendar year between 1601 and 2100
33			doy : (int)
34			Day of year between 1-365/366
35			utime : (float)
36			Seconds since midnight on the specified day
37
38			Returns
39			---------
40			sbsllon : (float)
41			Subsolar longitude in degrees E for the given date/time
42			sbsllat : (float)
43			Subsolar latitude in degrees N for the given date/time
44
45			Notes
46			--------
47			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
48			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
49			inclusive. According to the Almanac, results are good to at least 0.01
50			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
51			Accuracy for other years has not been tested. Every day is assumed to have
52			exactly 86400 seconds; thus leap seconds that sometimes occur on December
53			31 are ignored (their effect is below the accuracy threshold of the
54			algorithm).
55			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
56			by K. Laundal.
57
58			"""
59			warnings.warn(dep_str, category=FutureWarning)
60
61			# Convert from 4 digit year to 2 digit year
62			yr2 = year - 2000
63
64			if year >= 2101:
65			aacgmv2.logger.error('subsol invalid after 2100. Input year is:', year)
66
67			# Determine if this year is a leap year
68			nleap = np.floor((year - 1601) / 4)
69			nleap = nleap - 99
70			if year <= 1900:
71			if year <= 1600:
72			print('subsol.py: subsol invalid before 1601. Input year is:', year)
73			ncent = np.floor((year - 1601) / 100)
74			ncent = 3 - ncent
75			nleap = nleap + ncent
76
77			# Calculate some of the coefficients needed to deterimine the mean longitude
78			# of the sun and the mean anomaly
79			l_0 = -79.549 + (-0.238699 * (yr2 - 4 * nleap) + 3.08514e-2 * nleap)
80			g_0 = -2.472 + (-0.2558905 * (yr2 - 4 * nleap) - 3.79617e-2 * nleap)
81
82			# Days (including fraction) since 12 UT on January 1 of IYR2:
83			dfrac = (utime / 86400 - 1.5) + doy
84
85			# Mean longitude of Sun:
86			l_sun = l_0 + 0.9856474 * dfrac
87
88			# Mean anomaly:
89			grad = np.radians(g_0 + 0.9856003 * dfrac)
90
91			# Ecliptic longitude:
92			lmrad = np.radians(l_sun + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
93			sinlm = np.sin(lmrad)
94
95			# Days (including fraction) since 12 UT on January 1 of 2000:
96			epoch_day = dfrac + 365.0 * yr2 + nleap
97
98			# Obliquity of ecliptic:
99			epsrad = np.radians(23.439 - 4.0e-7 * epoch_day)
100
101			# Right ascension:
102			alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))
103
104			# Declination, which is the subsolar latitude:
105			sbsllat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))
106
107			# Equation of time (degrees):
108			etdeg = l_sun - alpha
109			etdeg = etdeg - 360.0 * np.round(etdeg / 360.0)
110
111			# Apparent time (degrees):
112			aptime = utime / 240.0 + etdeg # Earth rotates one degree every 240 s.
113
114			# Subsolar longitude:
115			sbsllon = 180.0 - aptime
116			sbsllon = sbsllon - 360.0 * np.round(sbsllon / 360.0)
117
118			return sbsllon, sbsllat
119
120			def gc2gd_lat(gc_lat):
121			"""Convert geocentric latitude to geodetic latitude using WGS84.
122
123			Parameters
124			-----------
125			gc_lat : (array_like or float)
126			Geocentric latitude in degrees N
127
128			Returns
129			---------
130			gd_lat : (same as input)
131			Geodetic latitude in degrees N
132			"""
133			warnings.warn(dep_str, category=FutureWarning)
134
135
136			wgs84_e2 = 0.006694379990141317 - 1.0
137			return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / wgs84_e2))
138
139			def igrf_dipole_axis(date):
140			"""Get Cartesian unit vector pointing at dipole pole in the north,
141			according to IGRF
142
143			Parameters
144			-------------
145			date : (dt.datetime)
146			Date and time
147
148			Returns
149			----------
150			m_0: (np.ndarray)
151			Cartesian 3 element unit vector pointing at dipole pole in the north
152			(geocentric coords)
153
154			Notes
155			----------
156			IGRF coefficients are read from the igrf12coeffs.txt file. It should also
157			work after IGRF updates. The dipole coefficients are interpolated to the
158			date, or extrapolated if date > latest IGRF model
159			"""
160			warnings.warn(dep_str, category=FutureWarning)
161
162			# get time in years, as float:
163			year = date.year
164			doy = date.timetuple().tm_yday
165			year_days = int(dt.date(date.year, 12, 31).strftime("%j"))
166			year = year + doy / year_days
167
168			# read the IGRF coefficients
169			with open(aacgmv2.IGRF_COEFFS, 'r') as f_igrf:
170			lines = f_igrf.readlines()
171
172			years = lines[3].split()[3:][:-1]
173			years = np.array(years, dtype=float) # time array
174
175			g10 = lines[4].split()[3:]
176			g11 = lines[5].split()[3:]
177			h11 = lines[6].split()[3:]
178
179			# secular variation coefficients (for extrapolation)
180			g10sv = np.float32(g10[-1])
181			g11sv = np.float32(g11[-1])
182			h11sv = np.float32(h11[-1])
183
184			# model coefficients:
185			g10 = np.array(g10[:-1], dtype=float)
186			g11 = np.array(g11[:-1], dtype=float)
187			h11 = np.array(h11[:-1], dtype=float)
188
189			# get the gauss coefficient at given time:
190			if year <= years[-1]:
191			# regular interpolation
192			g10 = np.interp(year, years, g10)
193			g11 = np.interp(year, years, g11)
194			h11 = np.interp(year, years, h11)
195			else:
196			# extrapolation
197			dyear = year - years[-1]
198			g10 = g10[-1] + g10sv * dyear
199			g11 = g11[-1] + g11sv * dyear
200			h11 = h11[-1] + h11sv * dyear
201
202			# calculate pole position
203			B_0 = np.sqrt(g102 + g112 + h11**2)
204
205			# Calculate output
206			m_0 = -np.array([g11, h11, g10]) / B_0
207
208			return m_0
209

aburrell / aacgmv2

Pull Request — develop (#47)

aacgmv2.deprecated.igrf_dipole_axis() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like