aacgmv2.utils.igrf_dipole_axis() - Code Metrics - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

aacgmv2.utils.igrf_dipole_axis() A
last analyzed 2024-12-27 21:19 UTC

↳ Parent: aacgmv2.utils

Complexity

Conditions

Size

Total Lines	68
Code Lines	29

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	29
dl	0
loc	68
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.

"""

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 : array-like or float
        Geodetic latitude in degrees N, same type as input `gc_lat`

    """
    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):
    """Find 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 (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) 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			import datetime as dt
17			import numpy as np
18
19			import aacgmv2
20
21
22			def gc2gd_lat(gc_lat):
23			"""Convert geocentric latitude to geodetic latitude using WGS84.
24
25			Parameters
26			-----------
27			gc_lat : array-like or float
28			Geocentric latitude in degrees N
29
30			Returns
31			---------
32			gd_lat : array-like or float
33			Geodetic latitude in degrees N, same type as input `gc_lat`
34
35			"""
36			wgs84_e2 = 0.006694379990141317 - 1.0
37			gd_lat = np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / wgs84_e2))
38
39			return gd_lat
40
41
42			def subsol(year, doy, utime):
43			"""Find subsolar geocentric longitude and latitude.
44
45			Parameters
46			----------
47			year : int
48			Calendar year between 1601 and 2100
49			doy : int
50			Day of year between 1-365/366
51			utime : float
52			Seconds since midnight on the specified day
53
54			Returns
55			-------
56			sbsllon : float
57			Subsolar longitude in degrees E for the given date/time
58			sbsllat : float
59			Subsolar latitude in degrees N for the given date/time
60
61			Raises
62			------
63			ValueError
64			If year is out of range
65
66			Notes
67			-----
68			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
69			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
70			inclusive. According to the Almanac, results are good to at least 0.01
71			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
72			Accuracy for other years has not been tested. Every day is assumed to have
73			exactly 86400 seconds; thus leap seconds that sometimes occur on December
74			31 are ignored (their effect is below the accuracy threshold of the
75			algorithm).
76
77			References
78			----------
79			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
80			by K. Laundal.
81
82			"""
83			# Convert from 4 digit year to 2 digit year
84			yr2 = year - 2000
85
86			if year >= 2101 or year <= 1600:
87			raise ValueError('subsol valid between 1601-2100. Input year is:', year)
88
89			# Determine if this year is a leap year
90			nleap = np.floor((year - 1601) / 4)
91			nleap = nleap - 99
92			if year <= 1900:
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 (IGRF).
143
144			Parameters
145			----------
146			date : dt.datetime
147			Date and time
148
149			Returns
150			-------
151			m_0 : np.ndarray
152			Cartesian 3 element unit vector pointing at dipole pole in the north
153			(geocentric coords)
154
155			Notes
156			-----
157			IGRF coefficients are read from the igrf12coeffs.txt file. It should also
158			work after IGRF updates. The dipole coefficients are interpolated to the
159			date, or extrapolated if date > latest IGRF model
160
161			"""
162			# Get time in years, as float
163			year = date.year
164			doy = date.timetuple().tm_yday
165			year_days = dt.date(date.year, 12, 31).timetuple().tm_yday
166			year = year + doy / year_days
167
168			# Read the IGRF coefficients
169			with open(aacgmv2.IGRF_COEFFS) 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] and year >= years[0]:
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

aacgmv2.utils.igrf_dipole_axis() A last analyzed 2024-12-27 21:19 UTC

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

aacgmv2.utils.igrf_dipole_axis() A
last analyzed 2024-12-27 21:19 UTC