aacgmv2.deprecated.subsol() - Code Metrics - Inspection of "Merge pull request #32 from aburrell/deprecation_w..." - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — develop ( 04a2cb...7e7468 )

by Angeline

created 2019-01-14 12:54 UTC

aacgmv2.deprecated.subsol() A

↳ Parent: aacgmv2.deprecated

Complexity

Conditions

Size

Total Lines	90
Code Lines	31

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	31
dl	0
loc	90
rs	9.1359
c	0
b	0
f	0
cc	4
nop	3

How to fix Long Method

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

Functions
-------------------------------------------------------------------------------
convert : Converts array location
subsol : finds subsolar geocentric longitude and latitude
gc2gd_lat : Convert between geocentric and geodetic coordinates
igrf_dipole_axis : Get Cartesian unit vector pointing at the IGRF north dipole
------------------------------------------------------------------------------

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 numpy as np
import logbook as logging
import warnings

def convert(lat, lon, alt, date=None, a2g=False, trace=False, allowtrace=False,
            badidea=False, geocentric=False):
    """Converts between geomagnetic coordinates and AACGM coordinates

    Parameters
    ------------
    lat : (float)
        Input latitude in degrees N (code specifies type of latitude)
    lon : (float)
        Input longitude in degrees E (code specifies type of longitude)
    alt : (float)
        Altitude above the surface of the earth in km
    date : (datetime)
        Datetime for magnetic field
    a2g : (bool)
        True for AACGM-v2 to geographic (geodetic), False otherwise
        (default=False)
    trace : (bool)
        If True, use field-line tracing, not coefficients (default=False)
    allowtrace : (bool)
        If True, use trace only above 2000 km (default=False)
    badidea : (bool)
        If True, use coefficients above 2000 km (default=False)
    geocentric : (bool)
        True for geodetic, False for geocentric w/RE=6371.2 (default=False)

    Returns
    -------
    lat_out : (float)
        Output latitude in degrees N
    lon_out : (float)
        Output longitude in degrees E
    """
    import aacgmv2

    dstr = "Deprecated routine, may be removed in future versions.  Recommend "
    dstr += "using convert_latlon or convert_latlon_arr"
    warnings.warn(dstr, category=FutureWarning)

    if(np.array(alt).max() > 2000 and not trace and not allowtrace and
       not badidea):
        estr = 'coefficients are not valid for altitudes above 2000 km. You'
        estr += ' must either use field-line tracing (trace=True '
        estr += 'or allowtrace=True) or indicate you know this is a bad idea'
        estr += ' (badidea=True)'
        logging.error(estr)
        raise ValueError(estr)

    # construct a code from the boolian flags
    bit_code = aacgmv2.convert_bool_to_bit(a2g=a2g, trace=trace,
                                           allowtrace=allowtrace,
                                           badidea=badidea,
                                           geocentric=geocentric)

    # convert location
    lat_out, lon_out, _ = aacgmv2.convert_latlon_arr(lat, lon, alt, date,
                                                     code=bit_code)

    return lat_out, lon_out

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.
    """
    dstr = "Deprecated routine, may be removed in future versions"
    warnings.warn(dstr, category=FutureWarning)

    
    yr2 = year - 2000

    if year >= 2101:
        logging.error('subsol invalid after 2100. Input year is:', 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

    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
    """
    dstr = "Deprecated routine, may be removed in future versions"
    warnings.warn(dstr, 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
    """
    import datetime as dt
    import aacgmv2

    dstr = "Deprecated routine, may be removed in future versions"
    warnings.warn(dstr, 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			# -- coding: utf-8 --
2			"""Pythonic wrappers for AACGM-V2 C functions that were depricated in the
3			change from version 2.0.0 to version 2.0.2
4
5			Functions
6			-------------------------------------------------------------------------------
7			convert : Converts array location
8			subsol : finds subsolar geocentric longitude and latitude
9			gc2gd_lat : Convert between geocentric and geodetic coordinates
10			igrf_dipole_axis : Get Cartesian unit vector pointing at the IGRF north dipole
11			------------------------------------------------------------------------------
12
13			References
14			-------------------------------------------------------------------------------
15			Laundal, K. M. and A. D. Richmond (2016), Magnetic Coordinate Systems, Space
16			Sci. Rev., doi:10.1007/s11214-016-0275-y.
17			-------------------------------------------------------------------------------
18			"""
19
20			from __future__ import division, absolute_import, unicode_literals
21			import numpy as np
22			import logbook as logging
23			import warnings
24
25			def convert(lat, lon, alt, date=None, a2g=False, trace=False, allowtrace=False,
26			badidea=False, geocentric=False):
27			"""Converts between geomagnetic coordinates and AACGM coordinates
28
29			Parameters
30			------------
31			lat : (float)
32			Input latitude in degrees N (code specifies type of latitude)
33			lon : (float)
34			Input longitude in degrees E (code specifies type of longitude)
35			alt : (float)
36			Altitude above the surface of the earth in km
37			date : (datetime)
38			Datetime for magnetic field
39			a2g : (bool)
40			True for AACGM-v2 to geographic (geodetic), False otherwise
41			(default=False)
42			trace : (bool)
43			If True, use field-line tracing, not coefficients (default=False)
44			allowtrace : (bool)
45			If True, use trace only above 2000 km (default=False)
46			badidea : (bool)
47			If True, use coefficients above 2000 km (default=False)
48			geocentric : (bool)
49			True for geodetic, False for geocentric w/RE=6371.2 (default=False)
50
51			Returns
52			-------
53			lat_out : (float)
54			Output latitude in degrees N
55			lon_out : (float)
56			Output longitude in degrees E
57			"""
58			import aacgmv2
59
60			dstr = "Deprecated routine, may be removed in future versions. Recommend "
61			dstr += "using convert_latlon or convert_latlon_arr"
62			warnings.warn(dstr, category=FutureWarning)
63
64			if(np.array(alt).max() > 2000 and not trace and not allowtrace and
65			not badidea):
66			estr = 'coefficients are not valid for altitudes above 2000 km. You'
67			estr += ' must either use field-line tracing (trace=True '
68			estr += 'or allowtrace=True) or indicate you know this is a bad idea'
69			estr += ' (badidea=True)'
70			logging.error(estr)
71			raise ValueError(estr)
72
73			# construct a code from the boolian flags
74			bit_code = aacgmv2.convert_bool_to_bit(a2g=a2g, trace=trace,
75			allowtrace=allowtrace,
76			badidea=badidea,
77			geocentric=geocentric)
78
79			# convert location
80			lat_out, lon_out, _ = aacgmv2.convert_latlon_arr(lat, lon, alt, date,
81			code=bit_code)
82
83			return lat_out, lon_out
84
85			def subsol(year, doy, utime):
86			"""Finds subsolar geocentric longitude and latitude.
87
88			Parameters
89			------------
90			year : (int)
91			Calendar year between 1601 and 2100
92			doy : (int)
93			Day of year between 1-365/366
94			utime : (float)
95			Seconds since midnight on the specified day
96
97			Returns
98			---------
99			sbsllon : (float)
100			Subsolar longitude in degrees E for the given date/time
101			sbsllat : (float)
102			Subsolar latitude in degrees N for the given date/time
103
104			Notes
105			--------
106			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
107			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
108			inclusive. According to the Almanac, results are good to at least 0.01
109			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
110			Accuracy for other years has not been tested. Every day is assumed to have
111			exactly 86400 seconds; thus leap seconds that sometimes occur on December
112			31 are ignored (their effect is below the accuracy threshold of the
113			algorithm).
114			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
115			by K. Laundal.
116			"""
117			dstr = "Deprecated routine, may be removed in future versions"
118			warnings.warn(dstr, category=FutureWarning)
119
120
121			yr2 = year - 2000
122
123			if year >= 2101:
124			logging.error('subsol invalid after 2100. Input year is:', year)
125
126			nleap = np.floor((year - 1601) / 4)
127			nleap = nleap - 99
128			if year <= 1900:
129			if year <= 1600:
130			print('subsol.py: subsol invalid before 1601. Input year is:', year)
131			ncent = np.floor((year - 1601) / 100)
132			ncent = 3 - ncent
133			nleap = nleap + ncent
134
135			l_0 = -79.549 + (-0.238699 * (yr2 - 4 * nleap) + 3.08514e-2 * nleap)
136			g_0 = -2.472 + (-0.2558905 * (yr2 - 4 * nleap) - 3.79617e-2 * nleap)
137
138			# Days (including fraction) since 12 UT on January 1 of IYR2:
139			dfrac = (utime / 86400 - 1.5) + doy
140
141			# Mean longitude of Sun:
142			l_sun = l_0 + 0.9856474 * dfrac
143
144			# Mean anomaly:
145			grad = np.radians(g_0 + 0.9856003 * dfrac)
146
147			# Ecliptic longitude:
148			lmrad = np.radians(l_sun + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
149			sinlm = np.sin(lmrad)
150
151			# Days (including fraction) since 12 UT on January 1 of 2000:
152			epoch_day = dfrac + 365.0 * yr2 + nleap
153
154			# Obliquity of ecliptic:
155			epsrad = np.radians(23.439 - 4.0e-7 * epoch_day)
156
157			# Right ascension:
158			alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))
159
160			# Declination, which is the subsolar latitude:
161			sbsllat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))
162
163			# Equation of time (degrees):
164			etdeg = l_sun - alpha
165			etdeg = etdeg - 360.0 * np.round(etdeg / 360.0)
166
167			# Apparent time (degrees):
168			aptime = utime / 240.0 + etdeg # Earth rotates one degree every 240 s.
169
170			# Subsolar longitude:
171			sbsllon = 180.0 - aptime
172			sbsllon = sbsllon - 360.0 * np.round(sbsllon / 360.0)
173
174			return sbsllon, sbsllat
175
176			def gc2gd_lat(gc_lat):
177			"""Convert geocentric latitude to geodetic latitude using WGS84.
178
179			Parameters
180			-----------
181			gc_lat : (array_like or float)
182			Geocentric latitude in degrees N
183
184			Returns
185			---------
186			gd_lat : (same as input)
187			Geodetic latitude in degrees N
188			"""
189			dstr = "Deprecated routine, may be removed in future versions"
190			warnings.warn(dstr, category=FutureWarning)
191
192
193			wgs84_e2 = 0.006694379990141317 - 1.0
194			return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / wgs84_e2))
195
196			def igrf_dipole_axis(date):
197			"""Get Cartesian unit vector pointing at dipole pole in the north,
198			according to IGRF
199
200			Parameters
201			-------------
202			date : (dt.datetime)
203			Date and time
204
205			Returns
206			----------
207			m_0: (np.ndarray)
208			Cartesian 3 element unit vector pointing at dipole pole in the north
209			(geocentric coords)
210
211			Notes
212			----------
213			IGRF coefficients are read from the igrf12coeffs.txt file. It should also
214			work after IGRF updates. The dipole coefficients are interpolated to the
215			date, or extrapolated if date > latest IGRF model
216			"""
217			import datetime as dt
218			import aacgmv2
219
220			dstr = "Deprecated routine, may be removed in future versions"
221			warnings.warn(dstr, category=FutureWarning)
222
223
224			# get time in years, as float:
225			year = date.year
226			doy = date.timetuple().tm_yday
227			year_days = int(dt.date(date.year, 12, 31).strftime("%j"))
228			year = year + doy / year_days
229
230			# read the IGRF coefficients
231			with open(aacgmv2.IGRF_COEFFS, 'r') as f_igrf:
232			lines = f_igrf.readlines()
233
234			years = lines[3].split()[3:][:-1]
235			years = np.array(years, dtype=float) # time array
236
237			g10 = lines[4].split()[3:]
238			g11 = lines[5].split()[3:]
239			h11 = lines[6].split()[3:]
240
241			# secular variation coefficients (for extrapolation)
242			g10sv = np.float32(g10[-1])
243			g11sv = np.float32(g11[-1])
244			h11sv = np.float32(h11[-1])
245
246			# model coefficients:
247			g10 = np.array(g10[:-1], dtype=float)
248			g11 = np.array(g11[:-1], dtype=float)
249			h11 = np.array(h11[:-1], dtype=float)
250
251			# get the gauss coefficient at given time:
252			if year <= years[-1]:
253			# regular interpolation
254			g10 = np.interp(year, years, g10)
255			g11 = np.interp(year, years, g11)
256			h11 = np.interp(year, years, h11)
257			else:
258			# extrapolation
259			dyear = year - years[-1]
260			g10 = g10[-1] + g10sv * dyear
261			g11 = g11[-1] + g11sv * dyear
262			h11 = h11[-1] + h11sv * dyear
263
264			# calculate pole position
265			B_0 = np.sqrt(g102 + g112 + h11**2)
266
267			# Calculate output
268			m_0 = -np.array([g11, h11, g10]) / B_0
269
270			return m_0
271

aburrell / aacgmv2

Push — develop ( 04a2cb...7e7468 )

aacgmv2.deprecated.subsol() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like