apexpy.helpers.subsol() - Code Metrics - Inspection of "RC v2.1.0" - aburrell/apexpy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — main (#147)

by Angeline

created 2024-12-30 23:46 UTC

apexpy.helpers.subsol() B

↳ Parent: apexpy.helpers

Complexity

Conditions

Size

Total Lines	101
Code Lines	45

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	45
dl	0
loc	101
rs	7.4
c	0
b	0
f	0
cc	7
nop	1

How to fix Long Method

# -*- coding: utf-8 -*-

"""This module contains helper functions used by :class:`~apexpy.Apex`."""

import datetime as dt
import numpy as np
import time


def set_array_float(in_val):
    """Set array data type to float.

    Parameters
    ----------
    in_val : any
        Input value, only modified if it is a np.ndarray

    Returns
    -------
    out_val : any
        Output value, if `in_val` was an array, `out_val` will be an array of
        type `np.float64`.

    """

    if isinstance(in_val, np.ndarray):
        out_val = in_val.astype(np.float64)
    else:
        out_val = in_val

    return out_val


def checklat(lat, name='lat'):
    """Makes sure the latitude is inside [-90, 90], clipping close values
    (tolerance 1e-4).

    Parameters
    ----------
    lat : array-like
        latitude
    name : str, optional
        parameter name to use in the exception message

    Returns
    -------
    lat : ndarray or float
        Same as input where values just outside the range have been
        clipped to [-90, 90]

    Raises
    ------
    ValueError
        if any values are too far outside the range [-90, 90]

    """
    if np.any(np.abs(lat) > 90 + 1e-5):
        raise ValueError('{:} must be in [-90, 90]'.format(name))

    return np.clip(lat, -90.0, 90.0)


def getsinIm(alat):
    """Computes sinIm from modified apex latitude.

    Parameters
    ----------
    alat : array-like
        Modified apex latitude

    Returns
    -------
    sinIm : ndarray or float

    """

    alat = np.float64(alat)

    return 2 * np.sin(np.radians(alat)) / np.sqrt(4 - 3
                                                  * np.cos(np.radians(alat))**2)


def getcosIm(alat):
    """Computes cosIm from modified apex latitude.

    Parameters
    ----------
    alat : array-like
        Modified apex latitude

    Returns
    -------
    cosIm : ndarray or float

    """

    alat = np.float64(alat)

    return np.cos(np.radians(alat)) / np.sqrt(4 - 3
                                              * np.cos(np.radians(alat))**2)


def toYearFraction(date):
    """Converts :class:`datetime.date` or :class:`datetime.datetime` to decimal
    year.

    Parameters
    ----------
    date : :class:`datetime.date` or :class:`datetime.datetime`
        Input date or datetime object

    Returns
    -------
    year : float
        Decimal year

    Notes
    -----
    The algorithm is taken from http://stackoverflow.com/a/6451892/2978652

    """

    def sinceEpoch(date):
        """returns seconds since epoch"""
        return time.mktime(date.timetuple())

    year = date.year
    startOfThisYear = dt.datetime(year=year, month=1, day=1)
    startOfNextYear = dt.datetime(year=year + 1, month=1, day=1)

    yearElapsed = sinceEpoch(date) - sinceEpoch(startOfThisYear)
    yearDuration = sinceEpoch(startOfNextYear) - sinceEpoch(startOfThisYear)
    fraction = yearElapsed / yearDuration

    year += fraction

    return year


def gc2gdlat(gclat):
    """Converts geocentric latitude to geodetic latitude using WGS84.

    Parameters
    ---------
    gclat : array-like
        Geocentric latitude

    Returns
    -------
    gdlat : ndarray or float
        Geodetic latitude

    """
    WGS84_e2 = 0.006694379990141317  # WGS84 first eccentricity squared
    return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gclat)) / (WGS84_e2 - 1)))


def subsol(datetime):
    """Finds subsolar geocentric latitude and longitude.

    Parameters
    ----------
    datetime : :class:`datetime.datetime` or :class:`numpy.ndarray[datetime64]`
        Date and time in UTC (naive objects are treated as UTC)

    Returns
    -------
    sbsllat : float
        Latitude of subsolar point
    sbsllon : float
        Longitude of subsolar point

    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.

    """
    # Convert to year, day of year and seconds since midnight
    if isinstance(datetime, dt.datetime):
        year = np.asanyarray([datetime.year])
        doy = np.asanyarray([datetime.timetuple().tm_yday])
        ut = np.asanyarray([datetime.hour * 3600 + datetime.minute * 60
                            + datetime.second + datetime.microsecond / 1.0e6])
    elif isinstance(datetime, np.ndarray):
        # This conversion works for datetime of wrong precision or unit epoch
        times = datetime.astype('datetime64[us]')
        year_floor = times.astype('datetime64[Y]')
        day_floor = times.astype('datetime64[D]')
        year = year_floor.astype(int) + 1970
        doy = (day_floor - year_floor).astype(int) + 1
        ut = (times.astype('datetime64[us]') - day_floor).astype(float)
        ut /= 1e6
    else:
        raise ValueError("input must be datetime.datetime or numpy array")

    if not (np.all(1601 <= year) and np.all(year <= 2100)):
        raise ValueError('Year must be in [1601, 2100]')

    yr = year - 2000

    nleap = np.floor((year - 1601.0) / 4.0).astype(int)
    nleap -= 99
    mask_1900 = year <= 1900
    if np.any(mask_1900):
        ncent = np.floor((year[mask_1900] - 1601.0) / 100.0).astype(int)
        ncent = 3 - ncent
        nleap[mask_1900] = nleap[mask_1900] + ncent

    l0 = -79.549 + (-0.238699 * (yr - 4.0 * nleap) + 3.08514e-2 * nleap)
    g0 = -2.472 + (-0.2558905 * (yr - 4.0 * nleap) - 3.79617e-2 * nleap)

    # Days (including fraction) since 12 UT on January 1 of IYR:
    df = (ut / 86400.0 - 1.5) + doy

    # Mean longitude of Sun:
    lmean = l0 + 0.9856474 * df

    # Mean anomaly in radians:
    grad = np.radians(g0 + 0.9856003 * df)

    # Ecliptic longitude:
    lmrad = np.radians(lmean + 1.915 * np.sin(grad)
                       + 0.020 * np.sin(2.0 * grad))
    sinlm = np.sin(lmrad)

    # Obliquity of ecliptic in radians:
    epsrad = np.radians(23.439 - 4e-7 * (df + 365 * yr + nleap))

    # Right ascension:
    alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))

    # Declination, which is also the subsolar latitude:
    sslat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))

    # Equation of time (degrees):
    etdeg = lmean - alpha
    nrot = np.round(etdeg / 360.0)
    etdeg = etdeg - 360.0 * nrot

    # Subsolar longitude calculation. Earth rotates one degree every 240 s.
    sslon = 180.0 - (ut / 240.0 + etdeg)
    nrot = np.round(sslon / 360.0)
    sslon = sslon - 360.0 * nrot

    # Return a single value from the output if the input was a single value
    if isinstance(datetime, dt.datetime):
        return sslat[0], sslon[0]
    return sslat, sslon


1			# -- coding: utf-8 --
2
3			"""This module contains helper functions used by :class:`~apexpy.Apex`."""
4
5			import datetime as dt
6			import numpy as np
7			import time
8
9
10			def set_array_float(in_val):
11			"""Set array data type to float.
12
13			Parameters
14			----------
15			in_val : any
16			Input value, only modified if it is a np.ndarray
17
18			Returns
19			-------
20			out_val : any
21			Output value, if `in_val` was an array, `out_val` will be an array of
22			type `np.float64`.
23
24			"""
25
26			if isinstance(in_val, np.ndarray):
27			out_val = in_val.astype(np.float64)
28			else:
29			out_val = in_val
30
31			return out_val
32
33
34			def checklat(lat, name='lat'):
35			"""Makes sure the latitude is inside [-90, 90], clipping close values
36			(tolerance 1e-4).
37
38			Parameters
39			----------
40			lat : array-like
41			latitude
42			name : str, optional
43			parameter name to use in the exception message
44
45			Returns
46			-------
47			lat : ndarray or float
48			Same as input where values just outside the range have been
49			clipped to [-90, 90]
50
51			Raises
52			------
53			ValueError
54			if any values are too far outside the range [-90, 90]
55
56			"""
57			if np.any(np.abs(lat) > 90 + 1e-5):
58			raise ValueError('{:} must be in [-90, 90]'.format(name))
59
60			return np.clip(lat, -90.0, 90.0)
61
62
63			def getsinIm(alat):
64			"""Computes sinIm from modified apex latitude.
65
66			Parameters
67			----------
68			alat : array-like
69			Modified apex latitude
70
71			Returns
72			-------
73			sinIm : ndarray or float
74
75			"""
76
77			alat = np.float64(alat)
78
79			return 2 * np.sin(np.radians(alat)) / np.sqrt(4 - 3
80			* np.cos(np.radians(alat))**2)
81
82
83			def getcosIm(alat):
84			"""Computes cosIm from modified apex latitude.
85
86			Parameters
87			----------
88			alat : array-like
89			Modified apex latitude
90
91			Returns
92			-------
93			cosIm : ndarray or float
94
95			"""
96
97			alat = np.float64(alat)
98
99			return np.cos(np.radians(alat)) / np.sqrt(4 - 3
100			* np.cos(np.radians(alat))**2)
101
102
103			def toYearFraction(date):
104			"""Converts :class:`datetime.date` or :class:`datetime.datetime` to decimal
105			year.
106
107			Parameters
108			----------
109			date : :class:`datetime.date` or :class:`datetime.datetime`
110			Input date or datetime object
111
112			Returns
113			-------
114			year : float
115			Decimal year
116
117			Notes
118			-----
119			The algorithm is taken from http://stackoverflow.com/a/6451892/2978652
120
121			"""
122
123			def sinceEpoch(date):
124			"""returns seconds since epoch"""
125			return time.mktime(date.timetuple())
126
127			year = date.year
128			startOfThisYear = dt.datetime(year=year, month=1, day=1)
129			startOfNextYear = dt.datetime(year=year + 1, month=1, day=1)
130
131			yearElapsed = sinceEpoch(date) - sinceEpoch(startOfThisYear)
132			yearDuration = sinceEpoch(startOfNextYear) - sinceEpoch(startOfThisYear)
133			fraction = yearElapsed / yearDuration
134
135			year += fraction
136
137			return year
138
139
140			def gc2gdlat(gclat):
141			"""Converts geocentric latitude to geodetic latitude using WGS84.
142
143			Parameters
144			---------
145			gclat : array-like
146			Geocentric latitude
147
148			Returns
149			-------
150			gdlat : ndarray or float
151			Geodetic latitude
152
153			"""
154			WGS84_e2 = 0.006694379990141317 # WGS84 first eccentricity squared
155			return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gclat)) / (WGS84_e2 - 1)))
156
157
158			def subsol(datetime):
159			"""Finds subsolar geocentric latitude and longitude.
160
161			Parameters
162			----------
163			datetime : :class:`datetime.datetime` or :class:`numpy.ndarray[datetime64]`
164			Date and time in UTC (naive objects are treated as UTC)
165
166			Returns
167			-------
168			sbsllat : float
169			Latitude of subsolar point
170			sbsllon : float
171			Longitude of subsolar point
172
173			Notes
174			-----
175			Based on formulas in Astronomical Almanac for the year 1996, p. C24.
176			(U.S. Government Printing Office, 1994). Usable for years 1601-2100,
177			inclusive. According to the Almanac, results are good to at least 0.01
178			degree latitude and 0.025 degrees longitude between years 1950 and 2050.
179			Accuracy for other years has not been tested. Every day is assumed to have
180			exactly 86400 seconds; thus leap seconds that sometimes occur on December
181			31 are ignored (their effect is below the accuracy threshold of the
182			algorithm).
183
184			After Fortran code by A. D. Richmond, NCAR. Translated from IDL
185			by K. Laundal.
186
187			"""
188			# Convert to year, day of year and seconds since midnight
189			if isinstance(datetime, dt.datetime):
190			year = np.asanyarray([datetime.year])
191			doy = np.asanyarray([datetime.timetuple().tm_yday])
192			ut = np.asanyarray([datetime.hour * 3600 + datetime.minute * 60
193			+ datetime.second + datetime.microsecond / 1.0e6])
194			elif isinstance(datetime, np.ndarray):
195			# This conversion works for datetime of wrong precision or unit epoch
196			times = datetime.astype('datetime64[us]')
197			year_floor = times.astype('datetime64[Y]')
198			day_floor = times.astype('datetime64[D]')
199			year = year_floor.astype(int) + 1970
200			doy = (day_floor - year_floor).astype(int) + 1
201			ut = (times.astype('datetime64[us]') - day_floor).astype(float)
202			ut /= 1e6
203			else:
204			raise ValueError("input must be datetime.datetime or numpy array")
205
206			if not (np.all(1601 <= year) and np.all(year <= 2100)):
207			raise ValueError('Year must be in [1601, 2100]')
208
209			yr = year - 2000
210
211			nleap = np.floor((year - 1601.0) / 4.0).astype(int)
212			nleap -= 99
213			mask_1900 = year <= 1900
214			if np.any(mask_1900):
215			ncent = np.floor((year[mask_1900] - 1601.0) / 100.0).astype(int)
216			ncent = 3 - ncent
217			nleap[mask_1900] = nleap[mask_1900] + ncent
218
219			l0 = -79.549 + (-0.238699 * (yr - 4.0 * nleap) + 3.08514e-2 * nleap)
220			g0 = -2.472 + (-0.2558905 * (yr - 4.0 * nleap) - 3.79617e-2 * nleap)
221
222			# Days (including fraction) since 12 UT on January 1 of IYR:
223			df = (ut / 86400.0 - 1.5) + doy
224
225			# Mean longitude of Sun:
226			lmean = l0 + 0.9856474 * df
227
228			# Mean anomaly in radians:
229			grad = np.radians(g0 + 0.9856003 * df)
230
231			# Ecliptic longitude:
232			lmrad = np.radians(lmean + 1.915 * np.sin(grad)
233			+ 0.020 * np.sin(2.0 * grad))
234			sinlm = np.sin(lmrad)
235
236			# Obliquity of ecliptic in radians:
237			epsrad = np.radians(23.439 - 4e-7 * (df + 365 * yr + nleap))
238
239			# Right ascension:
240			alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))
241
242			# Declination, which is also the subsolar latitude:
243			sslat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))
244
245			# Equation of time (degrees):
246			etdeg = lmean - alpha
247			nrot = np.round(etdeg / 360.0)
248			etdeg = etdeg - 360.0 * nrot
249
250			# Subsolar longitude calculation. Earth rotates one degree every 240 s.
251			sslon = 180.0 - (ut / 240.0 + etdeg)
252			nrot = np.round(sslon / 360.0)
253			sslon = sslon - 360.0 * nrot
254
255			# Return a single value from the output if the input was a single value
256			if isinstance(datetime, dt.datetime):
257			return sslat[0], sslon[0]
258			return sslat, sslon
259

aburrell / apexpy

Pull Request — main (#147)

apexpy.helpers.subsol() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like