aacgmv2.utils.igrf_dipole_axis() - Code Metrics - Inspection of "V2.6.2 Release" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#60)

by Angeline

created 2020-09-09 17:06 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	4
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

    Raises
    ------
    ValueError if year is out of range

    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 or year <= 1600:
        raise ValueError('subsol valid between 1601-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:
        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 = dt.date(date.year, 12, 31).timetuple().tm_yday
    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] and year >= years[0]:
        # 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			Raises
64			------
65			ValueError if year is out of range
66
67			Notes
68			-----
69			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
70			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
71			inclusive. According to the Almanac, results are good to at least 0.01
72			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
73			Accuracy for other years has not been tested. Every day is assumed to have
74			exactly 86400 seconds; thus leap seconds that sometimes occur on December
75			31 are ignored (their effect is below the accuracy threshold of the
76			algorithm).
77
78			References
79			----------
80			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
81			by K. Laundal.
82
83			"""
84
85			# Convert from 4 digit year to 2 digit year
86			yr2 = year - 2000
87
88			if year >= 2101 or year <= 1600:
89			raise ValueError('subsol valid between 1601-2100. Input year is:', year)
90
91			# Determine if this year is a leap year
92			nleap = np.floor((year - 1601) / 4)
93			nleap = nleap - 99
94			if year <= 1900:
95			ncent = np.floor((year - 1601) / 100)
96			ncent = 3 - ncent
97			nleap = nleap + ncent
98
99			# Calculate some of the coefficients needed to deterimine the mean longitude
100			# of the sun and the mean anomaly
101			l_0 = -79.549 + (-0.238699 * (yr2 - 4 * nleap) + 3.08514e-2 * nleap)
102			g_0 = -2.472 + (-0.2558905 * (yr2 - 4 * nleap) - 3.79617e-2 * nleap)
103
104			# Days (including fraction) since 12 UT on January 1 of IYR2:
105			dfrac = (utime / 86400 - 1.5) + doy
106
107			# Mean longitude of Sun:
108			l_sun = l_0 + 0.9856474 * dfrac
109
110			# Mean anomaly:
111			grad = np.radians(g_0 + 0.9856003 * dfrac)
112
113			# Ecliptic longitude:
114			lmrad = np.radians(l_sun + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
115			sinlm = np.sin(lmrad)
116
117			# Days (including fraction) since 12 UT on January 1 of 2000:
118			epoch_day = dfrac + 365.0 * yr2 + nleap
119
120			# Obliquity of ecliptic:
121			epsrad = np.radians(23.439 - 4.0e-7 * epoch_day)
122
123			# Right ascension:
124			alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))
125
126			# Declination, which is the subsolar latitude:
127			sbsllat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))
128
129			# Equation of time (degrees):
130			etdeg = l_sun - alpha
131			etdeg = etdeg - 360.0 * np.round(etdeg / 360.0)
132
133			# Apparent time (degrees):
134			aptime = utime / 240.0 + etdeg # Earth rotates one degree every 240 s.
135
136			# Subsolar longitude:
137			sbsllon = 180.0 - aptime
138			sbsllon = sbsllon - 360.0 * np.round(sbsllon / 360.0)
139
140			return sbsllon, sbsllat
141
142
143			def igrf_dipole_axis(date):
144			"""Get Cartesian unit vector pointing at dipole pole in the north,
145			according to IGRF
146
147			Parameters
148			----------
149			date : (dt.datetime)
150			Date and time
151
152			Returns
153			-------
154			m_0: (np.ndarray)
155			Cartesian 3 element unit vector pointing at dipole pole in the north
156			(geocentric coords)
157
158			Notes
159			-----
160			IGRF coefficients are read from the igrf12coeffs.txt file. It should also
161			work after IGRF updates. The dipole coefficients are interpolated to the
162			date, or extrapolated if date > latest IGRF model
163
164			"""
165
166			# get time in years, as float:
167			year = date.year
168			doy = date.timetuple().tm_yday
169			year_days = dt.date(date.year, 12, 31).timetuple().tm_yday
170			year = year + doy / year_days
171
172			# read the IGRF coefficients
173			with open(aacgmv2.IGRF_COEFFS, 'r') as f_igrf:
174			lines = f_igrf.readlines()
175
176			years = lines[3].split()[3:][:-1]
177			years = np.array(years, dtype=float) # time array
178
179			g10 = lines[4].split()[3:]
180			g11 = lines[5].split()[3:]
181			h11 = lines[6].split()[3:]
182
183			# secular variation coefficients (for extrapolation)
184			g10sv = np.float32(g10[-1])
185			g11sv = np.float32(g11[-1])
186			h11sv = np.float32(h11[-1])
187
188			# model coefficients:
189			g10 = np.array(g10[:-1], dtype=float)
190			g11 = np.array(g11[:-1], dtype=float)
191			h11 = np.array(h11[:-1], dtype=float)
192
193			# get the gauss coefficient at given time:
194			if year <= years[-1] and year >= years[0]:
195			# regular interpolation
196			g10 = np.interp(year, years, g10)
197			g11 = np.interp(year, years, g11)
198			h11 = np.interp(year, years, h11)
199			else:
200			# extrapolation
201			dyear = year - years[-1]
202			g10 = g10[-1] + g10sv * dyear
203			g11 = g11[-1] + g11sv * dyear
204			h11 = h11[-1] + h11sv * dyear
205
206			# calculate pole position
207			B_0 = np.sqrt(g102 + g112 + h11**2)
208
209			# Calculate output
210			m_0 = -np.array([g11, h11, g10]) / B_0
211
212			return m_0
213

aburrell / aacgmv2

Pull Request — master (#60)

aacgmv2.utils.igrf_dipole_axis() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like