aacgmv2.utils.igrf_dipole_axis() - Code Metrics - Inspection of "Merge branch 'deprecation' into develop" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — develop ( 885684...7d5043 )

by Angeline

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

aburrell / aacgmv2

Push — develop ( 885684...7d5043 )

aacgmv2.utils.igrf_dipole_axis() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like