meteortools.utils.Math

  1# Various maths functions scraped from WMPL
  2
  3import numpy as np
  4import math
  5from datetime import datetime, timedelta, MINYEAR
  6# Copyright (C) 2018-2023 Mark McIntyre
  7
  8# Define Julian epoch
  9JULIAN_EPOCH = datetime(2000, 1, 1, 12) # J2000.0 noon
 10J2000_JD = timedelta(2451545) # J2000.0 epoch in julian days
 11
 12
 13def jd2Date(jd, UT_corr=0, dt_obj=False):
 14    """ Converts the given Julian date to (year, month, day, hour, minute, second, millisecond) tuple or a datetime.  
 15
 16    Arguments:  
 17        jd: [float] Julian date  
 18
 19    Keyword arguments:  
 20        UT_corr: [float] UT correction in hours (difference from local time to UT)  
 21        dt_obj: [bool] returns a datetime object if True. False by default.  
 22
 23    Return:  
 24        (year, month, day, hour, minute, second, millisecond) or a datetime
 25
 26    """
 27
 28    try:
 29
 30        dt = timedelta(days=jd)
 31        
 32        date = dt + JULIAN_EPOCH - J2000_JD + timedelta(hours=UT_corr) 
 33
 34    # If the date is out of range (i.e. before year 1) use year 1. This is the limitation in the datetime
 35    # library. Time handling should be switched to astropy.time
 36    except OverflowError:
 37        date = datetime(MINYEAR, 1, 1, 0, 0, 0)
 38
 39
 40    # Return a datetime object if dt_obj == True
 41    if dt_obj:
 42        return date
 43
 44    return date.year, date.month, date.day, date.hour, date.minute, date.second, date.microsecond/1000.0
 45
 46
 47def datetime2JD(dt):
 48    """ Convert a datetime to a julian date  
 49    """
 50
 51    return date2JD(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, dt.microsecond/1000.0)
 52
 53
 54def date2JD(year, month, day, hour, minute, second, millisecond=0, UT_corr=0.0):
 55    """ Convert date and time to Julian Date in J2000.0.  
 56    
 57    Arguments:  
 58        year: [int] year  
 59        month: [int] month  
 60        day: [int] day of the date  
 61        hour: [int] hours  
 62        minute: [int] minutes  
 63        second: [int] seconds  
 64
 65    Kwargs:  
 66        millisecond: [int] milliseconds (optional)  
 67        UT_corr: [float] UT correction in hours (difference from local time to UT)  
 68    
 69    Return:  
 70        [float] julian date, J2000.0 epoch  
 71    """
 72
 73    # Convert all input arguments to integer (except milliseconds)
 74    year, month, day, hour, minute, second = map(int, (year, month, day, hour, minute, second))
 75
 76    # Create datetime object of current time
 77    dt = datetime(year, month, day, hour, minute, second, int(millisecond*1000))
 78
 79    # Calculate Julian date
 80    julian = dt - JULIAN_EPOCH + J2000_JD - timedelta(hours=UT_corr)
 81    
 82    # Convert seconds to day fractions
 83    return julian.days + (julian.seconds + julian.microseconds/1000000.0)/86400.0
 84
 85
 86def jd2LST(julian_date, lon):
 87    """ Convert Julian date to Local Sidereal Time and Greenwich Sidereal Time. The times used are apparent 
 88        times, not mean times.  
 89
 90    Source: J. Meeus: Astronomical Algorithms  
 91
 92    Arguments:  
 93        julian_date: [float] decimal julian date, epoch J2000.0  
 94        lon: [float] longitude of the observer in degrees  
 95    
 96    Return:  
 97        (LST, GST): [tuple of floats] a tuple of Local Sidereal Time and Greenwich Sidereal Time  
 98    """
 99
100    # t = (julian_date - J2000_JD.days)/36525.0
101
102    # Greenwich Sidereal Time
103    #GST = 280.46061837 + 360.98564736629*(julian_date - J2000_JD.days) + 0.000387933*t**2 - (t**3)/38710000
104    #GST = (GST + 360)%360
105
106    GST = np.degrees(calcApparentSiderealEarthRotation(julian_date))
107
108    # Local Sidereal Time
109    LST = (GST + lon + 360) % 360
110    
111    return LST, GST
112
113
114def jd2DynamicalTimeJD(jd):
115    """ Converts the given Julian date to dynamical time (i.e. Terrestrial Time, TT) Julian date. The 
116        conversion takes care of leap seconds.   
117
118    Arguments:  
119        jd: [float] Julian date.  
120
121    Return:  
122        [float] Dynamical time Julian date.  
123    """
124
125    # Leap seconds as of 2017 (default)
126    leap_secs = 37.0
127
128
129    # Get the relevant number of leap seconds for the given JD
130    # COMMENTED OUT as tai-utc.dat is empty
131    #for jd_leap, ls in config.leap_seconds:
132    #    
133    #    if jd > jd_leap:
134    #        leap_secs = ls
135            
136
137    # Calculate the dynamical JD
138    jd_dyn = jd + (leap_secs + 32.184)/86400.0
139
140
141    return jd_dyn
142
143
144def greatCircleDistance(lat1, lon1, lat2, lon2):
145    """ Calculate the great circle distance in kilometers between two points on the Earth. 
146        Source: https://gis.stackexchange.com/a/56589/15183  
147
148    Arguments:  
149        lat1: [float] Latitude 1 (radians).  
150        lon1: [float] Longitude 1 (radians).  
151        lat2: [float] Latitude 2 (radians).  
152        lon2: [float] Longitude 2 (radians).  
153
154    Return:  
155        [float]: Distance in kilometers.  
156    """
157    
158    # Haversine formula
159    dlon = lon2 - lon1 
160    dlat = lat2 - lat1 
161
162    a = np.sin(dlat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dlon/2)**2
163    c = 2*np.arcsin(np.sqrt(a))
164
165    # Distance in kilometers.
166    dist = 6371*c
167
168    return dist
169
170
171def angleBetweenSphericalCoords(phi1, lambda1, phi2, lambda2):
172    """ Calculates the angle between two points on a sphere.  
173    
174    Arguments:  
175        phi1: [float] Latitude 1 (radians).  
176        lambda1: [float] Longitude 1 (radians).  
177        phi2: [float] Latitude 2 (radians).  
178        lambda2: [float] Longitude 2 (radians).  
179
180    Return:  
181        [float] Angle between two coordinates (radians).  
182    """
183
184    return np.arccos(np.sin(phi1)*np.sin(phi2) + np.cos(phi1)*np.cos(phi2)*np.cos(lambda2 - lambda1))
185
186
187def equatorialCoordPrecession(start_epoch, final_epoch, ra, dec):
188    """ Precess right Ascension and declination from one epoch to another, taking only precession into 
189        account.  
190
191        Implemented from: Jean Meeus - Astronomical Algorithms, 2nd edition, pages 134-135  
192    
193    Arguments:  
194        start_epoch: [float] Julian date of the starting epoch  
195        final_epoch: [float] Julian date of the final epoch  
196        ra: [float] non-corrected right ascension in radians  
197        dec: [float] non-corrected declination in radians  
198    
199    Return:  
200        (ra, dec): [tuple of floats] precessed equatorial coordinates in radians  
201
202    """
203
204
205    T = (start_epoch - J2000_JD.days)/36525.0
206    t = (final_epoch - start_epoch)/36525.0
207
208    # Calculate correction parameters
209    zeta = ((2306.2181 + 1.39656*T - 0.000139*T**2)*t + (0.30188 - 0.000344*T)*t**2 + 0.017998*t**3)/3600
210    z = ((2306.2181 + 1.39656*T - 0.000139*T**2)*t + (1.09468 + 0.000066*T)*t**2 + 0.018203*t**3)/3600
211    theta = ((2004.3109 - 0.85330*T - 0.000217*T**2)*t - (0.42665 + 0.000217*T)*t**2 - 0.041833*t**3)/3600
212
213    # Convert parameters to radians
214    zeta, z, theta = map(math.radians, (zeta, z, theta))
215
216    # Calculate the next set of parameters
217    A = np.cos(dec) * np.sin(ra + zeta)
218    B = np.cos(theta)* np.cos(dec) * np.cos(ra + zeta) - np.sin(theta) * np.sin(dec)
219    C = np.sin(theta)* np.cos(dec) * np.cos(ra + zeta) + np.cos(theta) * np.sin(dec)
220
221    # Calculate right ascension
222    ra_corr = np.arctan2(A, B) + z
223
224    # Calculate declination (apply a different equation if close to the pole, closer then 0.5 degrees)
225    if (np.pi/2 - np.abs(dec)) < np.radians(0.5):
226        dec_corr = np.sign(C)*np.arccos(np.sqrt(A**2 + B**2))
227    else:
228        dec_corr = np.arcsin(C)
229
230    # Wrap right ascension to [0, 2*pi] range
231    ra_corr = ra_corr % (2*np.pi)
232
233    # Wrap declination to [-pi/2, pi/2] range
234    dec_corr = (dec_corr + np.pi/2) % np.pi - np.pi/2
235
236    return ra_corr, dec_corr
237
238
239equatorialCoordPrecession_vect = np.vectorize(equatorialCoordPrecession, excluded=['start_epoch'])
240"""Vectorize the equatorialCoordPrecession, so ra and dec can be passed as numpy arrays"""
241
242
243def raDec2AltAz(ra, dec, jd, lat, lon):
244    """ Convert right ascension and declination to azimuth (+east of sue north) and altitude.  
245
246    Arguments:  
247        ra: [float] right ascension in radians  
248        dec: [float] declination in radians  
249        jd: [float] Julian date  
250        lat: [float] latitude in radians  
251        lon: [float] longitude in radians  
252
253    Return:  
254        (azim, elev): [tuple]  
255            azim: [float] azimuth (+east of due north) in radians  
256            elev: [float] elevation above horizon in radians  
257
258        """
259
260    # Calculate Local Sidereal Time
261    lst = np.radians(jd2LST(jd, np.degrees(lon))[0])
262
263    # Calculate the hour angle
264    ha = lst - ra
265
266    # Constrain the hour angle to [-pi, pi] range
267    ha = (ha + np.pi) % (2*np.pi) - np.pi
268
269    # Calculate the azimuth
270    azim = np.pi + np.arctan2(np.sin(ha), np.cos(ha)*np.sin(lat) - np.tan(dec)*np.cos(lat))
271
272    # Calculate the sine of elevation
273    sin_elev = np.sin(lat)*np.sin(dec) + np.cos(lat)*np.cos(dec)*np.cos(ha)
274
275    # Wrap the sine of elevation in the [-1, +1] range
276    sin_elev = (sin_elev + 1) % 2 - 1
277
278    elev = np.arcsin(sin_elev)
279
280    return azim, elev
281
282
283raDec2AltAz_vect = np.vectorize(raDec2AltAz, excluded=['lat', 'lon'])
284"""Vectorize the raDec2AltAz function so it can take numpy arrays for: ra, dec, jd"""
285
286
287def altAz2RADec(azim, elev, jd, lat, lon):
288    """ Convert azimuth and altitude in a given time and position on Earth to right ascension and 
289        declination.   
290
291    Arguments:  
292        azim: [float] azimuth (+east of due north) in radians  
293        elev: [float] elevation above horizon in radians  
294        jd: [float] Julian date  
295        lat: [float] latitude of the observer in radians  
296        lon: [float] longitde of the observer in radians  
297
298    Return:  
299        (RA, dec): [tuple]  
300            RA: [float] right ascension (radians)  
301            dec: [float] declination (radians)  
302    """
303    
304    # Calculate hour angle
305    ha = np.arctan2(-np.sin(azim), np.tan(elev)*np.cos(lat) - np.cos(azim)*np.sin(lat))
306
307    # Calculate Local Sidereal Time
308    lst = np.radians(jd2LST(jd, np.degrees(lon))[0])
309    
310    # Calculate right ascension
311    ra = (lst - ha) % (2*np.pi)
312
313    # Calculate declination
314    dec = np.arcsin(np.sin(lat)*np.sin(elev) + np.cos(lat)*np.cos(elev)*np.cos(azim))
315
316    return ra, dec
317
318
319altAz2RADec_vect = np.vectorize(altAz2RADec, excluded=['lat', 'lon'])
320""" Vectorize the altAz2RADec function so it can take numpy arrays for: azim, elev, jd """
321
322
323def calcApparentSiderealEarthRotation(julian_date):
324    """ Calculate apparent sidereal rotation GST of the Earth.  
325        
326        Calculated according to:  
327        Clark, D. L. (2010). Searching for fireball pre-detections in sky surveys. The School of Graduate and 
328        Postdoctoral Studies. University of Western Ontario, London, Ontario, Canada, MSc Thesis.  
329
330    """
331
332    t = (julian_date - J2000_JD.days)/36525.0
333
334    # Calculate the Mean sidereal rotation of the Earth in radians (Greenwich Sidereal Time)
335    GST = 280.46061837 + 360.98564736629*(julian_date - J2000_JD.days) + 0.000387933*t**2 - (t**3)/38710000
336    GST = (GST + 360) % 360
337    GST = math.radians(GST)
338
339    # print('GST:', np.degrees(GST), 'deg')
340
341    # Calculate the dynamical time JD
342    jd_dyn = jd2DynamicalTimeJD(julian_date)
343
344
345    # Calculate Earth's nutation components
346    delta_psi, delta_eps = calcNutationComponents(jd_dyn)
347
348    # print('Delta Psi:', np.degrees(delta_psi), 'deg')
349    # print('Delta Epsilon:', np.degrees(delta_eps), 'deg')
350
351
352    # Calculate the mean obliquity (in arcsec)
353    u = (jd_dyn - 2451545.0)/3652500.0
354    eps0 = 84381.448 - 4680.93*u - 1.55*u**2 + 1999.25*u**3 - 51.38*u**4 - 249.67*u**5 - 39.05*u**6 \
355        + 7.12*u**7 + 27.87*u**8 + 5.79*u**9 + 2.45*u**10
356
357    # Convert to radians
358    eps0 /= 3600
359    eps0 = np.radians(eps0)
360
361    # print('Mean obliquity:', np.degrees(eps0), 'deg')
362
363    # Calculate apparent sidereal Earth's rotation
364    app_sid_rot = (GST + delta_psi*math.cos(eps0 + delta_eps)) % (2*math.pi)
365
366    return app_sid_rot
367
368
369def calcNutationComponents(jd_dyn):
370    """ Calculate Earth's nutation components from the given Julian date.  
371
372    Source: Meeus (1998) Astronomical algorithms, 2nd edition, chapter 22.  
373    
374    The precision is limited to 0.5" in nutation in longitude and 0.1" in nutation in obliquity. The errata 
375    for the 2nd edition was used to correct the equation for delta_psi.
376    
377    Arguments:  
378        jd_dyn: [float] Dynamical Julian date. See wmpl.Utils.TrajConversions.jd2DynamicalTimeJD function.  
379
380    Return:  
381        (delta_psi, delta_eps): [tuple of floats] Differences from mean nutation due to the influence of
382            the Moon and minor effects (radians).  
383    """
384
385
386    T = (jd_dyn - J2000_JD.days)/36525.0
387
388    # # Mean Elongation of the Moon from the Sun
389    # D = 297.85036 + 445267.11148*T - 0.0019142*T**2 + (T**3)/189474
390
391    # # Mean anomaly of the Earth with respect to the Sun
392    # M = 357.52772 + 35999.05034*T - 0.0001603*T**2 - (T**3)/300000
393
394    # # Mean anomaly of the Moon
395    # Mm = 134.96298 + 477198.867398*T + 0.0086972*T**2 + (T**3)/56250
396
397    # # Argument of latitude of the Moon
398    # F = 93.27191  + 483202.017538*T - 0.0036825*T**2 + (T**3)/327270
399
400
401    # Longitude of the ascending node of the Moon's mean orbit on the ecliptic, measured from the mean equinox
402    # of the date
403    omega = 125.04452 - 1934.136261*T
404
405
406    # Mean longitude of the Sun (deg)
407    L = 280.4665 + 36000.7698*T
408
409    # Mean longitude of the Moon (deg)
410    Ll = 218.3165 + 481267.8813*T
411
412
413    # Nutation in longitude
414    delta_psi = -17.2*math.sin(math.radians(omega)) - 1.32*math.sin(np.radians(2*L)) \
415        - 0.23*math.sin(math.radians(2*Ll)) + 0.21*math.sin(math.radians(2*omega))
416
417    # Nutation in obliquity
418    delta_eps = 9.2*math.cos(math.radians(omega)) + 0.57*math.cos(math.radians(2*L)) \
419        + 0.1*math.cos(math.radians(2*Ll)) - 0.09*math.cos(math.radians(2*omega))
420
421
422    # Convert to radians
423    delta_psi = np.radians(delta_psi/3600)
424    delta_eps = np.radians(delta_eps/3600)
425
426    return delta_psi, delta_eps