Passed
Pull Request — main (#147)
by Angeline
01:39
created

apexpy.helpers.set_array_float()   A

Complexity

Conditions 2

Size

Total Lines 22
Code Lines 5

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 2
eloc 5
nop 1
dl 0
loc 22
rs 10
c 0
b 0
f 0
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