aacgmv2.utils.igrf_dipole_axis() - Code Metrics - Inspection of "Un-deprecation of utilities" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — develop (#57)

by Angeline

created 2020-08-31 13:59 UTC

aacgmv2.utils.igrf_dipole_axis() A

↳ Parent: aacgmv2.utils

Complexity

Conditions

Size

Total Lines	70
Code Lines	29

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	29
dl	0
loc	70
rs	9.184
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 -*-
"""utilities that support the AACGM-V2 C functions.

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 aacgmv2


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

    """

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

    return gd_lat


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).

    References
    ----------
    After Fortran code by A. D. Richmond, NCAR. Translated from IDL
    by K. Laundal.

    """

    # 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 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

    """

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

aburrell / aacgmv2

Pull Request — develop (#57)

aacgmv2.utils.igrf_dipole_axis() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like