-
Notifications
You must be signed in to change notification settings - Fork 0
/
saturnLs.py
executable file
·446 lines (381 loc) · 15.5 KB
/
saturnLs.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
#!/usr/bin/env python
from io import StringIO
import numpy as np
import pandas as pd
__doc__ = """
Functions to convert Saturn's solar longitude to and from SCET/UNIX timestamp, datetime, and Julian date.
Values:
* Ls: Solar longitude, apparent longitude of the sun seen from Saturn as defined in JPL Horizons.
* SY: Saturn year number where SY 1 is the Saturn year starting September 21, 1950.
* Ls2: Solar longitude since the beginning of SY 1 starting September 21, 1950.
* SCET: Time in seconds since January 1, 1970 (UNIX epoch).
* JD (or JDUT): Julian Date, UT.
* datetime: Python datetime object (including numpy datetime64, pandas Timestamp, and astropy Time).
* datestr (or datestring): Date and/or time as a string.
Written 2021 by LE Hanson.
"""
_HAS_SCIPY = True
try:
from scipy.interpolate import interp1d
except:
_HAS_SCIPY = False
__all__ = [
'datetime_to_Ls2', 'Ls2_to_datetime', 'datestr_to_Ls2', 'datestr_to_Ls',
'SCET_to_datetime', 'datetime_to_SCET',
'SCET_to_Ls2', 'Ls2_to_SCET',
'JD_to_SCET', 'SCET_to_JD',
'Ls2_to_Ls', 'Ls2_to_SYLs', 'SYLs_to_Ls2',
'dfa',
]
# calendar and JD date of start of several "Saturn years"
_sy_list = dict((
(1, {"date":"1774-01-03 02:22:37.246926848", "JDUT":2369002.59904221}),
(2, {"date":"1803-06-18 11:01:00.319862784", "JDUT":2379759.95903148}),
(3, {"date":"1832-12-03 18:51:45.533949952", "JDUT":2390521.28594368}),
(4, {"date":"1862-05-17 06:22:58.443740160", "JDUT":2401277.76595421}),
(5, {"date":"1891-10-30 07:07:31.777902603", "JDUT":2412035.7968955776}),
(6, {"date":"1921-04-12 00:26:22.458478928", "JDUT":2422791.5183154917}),
(7, {"date":"1950-09-21 19:02:40.934159994", "JDUT":2433546.2935293308}),
(8, {"date":"1980-03-03 16:12:43.076112986", "JDUT":2444302.175498566}),
(9, {"date":"2009-08-11 02:04:21.114835262", "JDUT":2455054.5863554957}),
(10, {"date":"2039-01-22 19:51:53.897613049", "JDUT":2465811.327707148}),
(11, {"date":"2068-06-29 08:45:54.395126820", "JDUT":2476561.8652129066}),
(12, {"date":"2097-12-13 07:54:47.293023586", "JDUT":2487320.82971404}),
(13, {"date":"2127-05-21 00:44:11.711929321", "JDUT":2498070.53069111}),
(14, {"date":"2156-11-05 03:52:56.824035645", "JDUT":2508831.661768797}),
(15, {"date":"2186-04-12 16:53:45.123068928", "JDUT":2519582.20399448}),
(16, {"date":"2215-10-02 16:21:44.526243840", "JDUT":2530346.18176535}),
(17, {"date":"2245-03-09 17:53:51.117792256", "JDUT":2541097.24573053}),
))
# DataFrame holding the date each "Saturn year" begins
dfsy = pd.DataFrame.from_dict(_sy_list, orient='index').rename_axis(index="SY")
dfsy['date'] = pd.to_datetime(dfsy['date'])
# ephemeris files from JPL Horizons
ephem_file_coarse = "saturn-Ls-1890_2160-sol_sparse.txt"
ephem_file_fine = "saturn-Ls-1970_2040-sol.txt"
# define default horizons output file to load
ephem_file = ephem_file_coarse
# csv ephemeris files, reduced from Horizons output
ephem_csv_long = "saturn-Ls-1890_2160-sol_sparse.csv"
ephem_csv_coarse = "saturn-Ls-sparse.csv"
ephem_csv_fine = "saturn-Ls.csv"
# define default csv file to load
ephem_csv = ephem_csv_fine
###################################
# Independent conversion functions
###################################
# these functions do not rely on the ephemeris data
# Ls2 and Ls, SY
def Ls2_to_Ls(ls2):
'''Convert Ls2 to Ls.'''
SY = ls2//360
return ls2 - 360*SY
def Ls2_to_SYLs(ls2):
'''Convert Ls2 to "Saturn year" SY and Ls.'''
SY = ls2//360 + 1
return SY, ls2 - 360*(SY-1)
def SYLs_to_Ls2(Ls, SY=3):
'''Convert "Saturn year" SY and Ls to Ls2.'''
return (SY-1)*360 + Ls
# convert SCET and datetime
def SCET_to_datetime(scet):
'''Convert SCET (seconds since UNIX epoch) to datetime'''
return pd.to_datetime(scet, unit='s')
def pddt_scet(pddt):
'''Convert python datetime to SCET.
This function strips any timezone information.'''
return (pddt - pd.to_datetime("1970")).total_seconds()
pddt_scet = np.frompyfunc(pddt_scet, nin=1, nout=1)
def datetime_to_SCET(pdt):
'''Convert datetime to SCET'''
from datetime import datetime
if isinstance(pdt, datetime):
pdt = pdt.replace(tzinfo=None)
return pddt_scet(pd.to_datetime(pdt))
datetime_to_SCET = np.frompyfunc(datetime_to_SCET, nin=1, nout=1)
###########################
# Ephemeris file functions
###########################
def scan_ephem(fname, header=False):
'''Scan a JPL Horizons ephemeris file (csv format) and return the csv data.'''
'''
expects ephemeris format returned by the following batch script,
send in email to horizons@ssd.jpl.nasa.gov with subject "BATCH-LONG"
!$$SOF
MAKE_EPHEM=YES
COMMAND=699
EPHEM_TYPE=OBSERVER
CENTER='500@399'
START_TIME='1750-01-01'
STOP_TIME='1800-01-01'
STEP_SIZE='12 HOURS'
QUANTITIES='41,44'
REF_SYSTEM='ICRF'
CAL_FORMAT='BOTH'
TIME_DIGITS='FRACSEC'
ANG_FORMAT='HMS'
APPARENT='AIRLESS'
RANGE_UNITS='AU'
SUPPRESS_RANGE_RATE='NO'
SKIP_DAYLT='NO'
SOLAR_ELONG='0,180'
EXTRA_PREC='NO'
RTS_ONLY='NO'
CSV_FORMAT='YES'
OBJ_DATA='YES'
'''
import re, gzip
if fname.endswith('.gz'):
with gzip.open(fname, "rb") as fin:
lines = []
for line in fin.readlines():
lines.append(line.decode())
else:
with open(fname, "r") as fin:
lines = fin.readlines()
soe = False
txt, hed = [], []
for ix in range(len(lines)):
ll = lines[ix]
if "$$SOE" in ll:
cols = lines[ix-2]
#cols = cols.replace(" ","").replace("n.a.","NaN")
cols = cols.strip().replace(r", , , ", r", ").replace("n.a.","NaN")
cols = re.sub(r",\s*$", "", cols)
cols = re.sub(",(\d[^,]+)", r",x\1", re.sub("[*\-]|_+", "_", re.sub("[%()./:]", "" ,cols)))
txt.append(cols.strip())
soe = True
elif "$$EOE" in ll:
soe = False
elif soe:
txt.append(ll.replace(", , , ", ", ").strip().strip(','))
elif header and ll[0] != '*':
hed.append(ll.strip().strip(','))
csv = "\n".join(txt).replace(", , , ", ", ")
if header:
return csv, "".join(hed)
return csv
def parse_ephem(s):
"""Parse the csv portion of JPL Horizon output (from scan_ephem) and return as a pandas DataFrame."""
df = pd.read_csv(StringIO(s), header=0, skipinitialspace=True, sep=',')
df = df.dropna("columns")
if "Calendar_Date_TDB" in df:
df["Calendar_Date_TDB"] = pd.to_datetime(df["Calendar_Date_TDB"].str[6:])
if "Date_UT_HRMNSCfff" in df:
df["Date_UT_HRMNSCfff"] = pd.to_datetime(df["Date_UT_HRMNSCfff"])
if "x1_way_down_LT" in df:
df["x1_way_down_LT"] = pd.to_timedelta(df["x1_way_down_LT"], 'min')
df = df.rename(columns={"x1_way_down_LT":"one_way_LT"})
return df[[c for c in df.columns if "Unnamed" not in c]]
def adjust_LT(df):
"""Account for light time.
Adjust for distance of observer to Saturn using the column "one_way_LT"."""
dfa = df.copy()
lt = dfa['one_way_LT'] # one-way light time
# calculate JDUT adjusted for light time
jdnew = dfa['JDUT'] - lt.dt.total_seconds()/86400 # correction to JDUT
ls2 = dfa["Ls2"]
# interpolate Ls2
if _HAS_SCIPY:
finterp = interp1d(jdnew, ls2, kind='linear', fill_value='extrapolate')
ls2new = finterp(dfa["JDUT"]) # calculate interpolated Ls2
else:
ls2new = np.interp(dfa["JDUT"], jdnew, ls2)
dls2 = ls2new - ls2 # difference in original Ls2 and new Ls2
dfa["Ls2"] = ls2new # replace Ls2 with corrected values
dfa["Ls"] += dls2 # correct Ls
# remove any nan's due to interpolation
dfa = dfa.dropna(subset=["Ls2", "Ls"])
# return without one_way_LT
return dfa.drop(columns=['one_way_LT'])
def make_df_Ls2SY(df):
"""Calculate Ls2, SY, and SCET, referenced from "Saturn year" zero starting in 1921.
Input is a DataFrame including fields ['date', 'JDUT', 'Ls']."""
_dfsy = dfsy.copy()
_dfsy["Ls"] = 0.0
dfa = pd.concat((df.copy(), _dfsy.reset_index())).sort_values('JDUT')
dfa["SY"] = dfa["SY"].ffill()
# if beyond SY table, estimate year
# could use Ls to find actual start of year but this is easy
if np.any(dfa.JDUT < _dfsy.JDUT.min()) or np.any(dfa.JDUT > _dfsy.JDUT.max()+10750):
jdyear = 10755.9 # julian days in saturn year
jdsy0 = 2369002.59904221-jdyear # SY 0 julian date
lx = (dfa.JDUT < _dfsy.JDUT.min()) | (dfa.JDUT > _dfsy.JDUT.max()+jdyear)
# estimate SY
dfa.loc[lx, "SY"] = np.floor((dfa.loc[lx, "JDUT"] - jdsy0)/jdyear)
dfa = dfa.set_index('JDUT').truncate(df.JDUT.min(), df.JDUT.max()).reset_index()
#dfa = dfa.dropna(thresh=2).reset_index(drop=True)
dfa["Ls2"] = dfa["Ls"] + 360*(dfa["SY"]-1)
dfa["SCET"] = (dfa.date - pd.to_datetime("1970")).dt.total_seconds()
# remove Ls = 0 rows
dfa = dfa.loc[~dfa["JDUT"].isin(_dfsy["JDUT"])]
# adjust for light time if column "one_way_LT" is present
if 'one_way_LT' in dfa:
dfa = adjust_LT(dfa) # adjust_LT removes "one_way_LT" column
return dfa[['date','JDUT','SCET','Ls2','Ls','SY']]
def make_Ls_df2(df):
"""Make the reference DataFrame from the ephemeris DataFrame."""
df = df.rename(columns={"Date_UT_HRMNSCfff":"date",
"App_Lon_Sun":"Ls",
"Date_JDUT":"JDUT"})
return make_df_Ls2SY(df)
def load_ephem(fname=ephem_file):
'''load JPL HORIZONS csv output as a DataFrame'''
return make_Ls_df2(parse_ephem(scan_ephem(fname)))
def save_csv_ephem(df, fname="saturn-Ls.csv", keep_date=False):
'''Save ephemeris DataFrame as a csv file.
This file loads significantly faster and is smaller than the raw Horizons output.'''
keepcols = ["JDUT", "Ls"]
if keep_date:
keepcols.insert(0, "date")
df[keepcols].to_csv(fname, index=False)
def load_csv_ephem(fname=ephem_csv):
'''Loads the ephemeris csv file saved by `save_csv_ephem`.'''
dfa = pd.read_csv(fname)
dfa["date"] = pd.to_datetime(dfa["JDUT"], unit='D', origin='julian')
return make_df_Ls2SY(dfa)
# load the ephemeris data (Horizons output or csv)
# this DataFrame is used in all of the Ls-time conversion functions by default
#dfa = load_ephem(ephem_file)
class dfa_store:
data = None
dfa = dfa_store()
dfa.data = load_csv_ephem()
#######################
# conversion functions
#######################
# convert Ls2 and SCET
def SCET_to_Ls2(scet, dfa=dfa):
"""Convert SCET to Ls2."""
if np.any(scet > dfa.data.SCET.max()) or np.any(scet < dfa.data.SCET.min()):
raise ValueError("SCET outside of range in data table.")
# interpolate Ls2
if _HAS_SCIPY:
return interp1d(dfa.data.SCET, dfa.data.Ls2, kind='linear')(scet)
return np.interp(scet, dfa.data.SCET, dfa.data.Ls2)
def Ls2_to_SCET(ls2, dfa=dfa):
"""Convert Ls2 to SCET."""
if np.any(ls2 > dfa.data.Ls2.max()) or np.any(ls2 < dfa.data.Ls2.min()):
raise ValueError("Ls2 outside of range in data table.")
# interpolate SCET
if _HAS_SCIPY:
return interp1d(dfa.data.Ls2, dfa.data.SCET, kind='linear')(ls2)
return np.interp(ls2, dfa.data.Ls2, dfa.data.SCET)
# convert Ls2 and JD
def JD_to_SCET(jd, dfa=dfa):
"""Convert JD to SCET."""
if np.any(jd > dfa.data.JDUT.max()) or np.any(jd < dfa.data.JDUT.min()):
raise ValueError("JD outside of range in data table.")
# interpolate SCET
if _HAS_SCIPY:
return interp1d(dfa.data.JDUT, dfa.data.SCET, kind='linear')(jd)
return np.interp(jd, dfa.data.JDUT, dfa.data.SCET)
def SCET_to_JD(scet, dfa=dfa):
"""Convert SCET to JD."""
if np.any(scet > dfa.data.SCET.max()) or np.any(scet < dfa.data.SCET.min()):
raise ValueError("SCET outside of range in data table.")
# interpolate JD
if _HAS_SCIPY:
return interp1d(dfa.data.SCET, dfa.data.JDUT, kind='linear')(scet)
return np.interp(scet, dfa.data.SCET, dfa.data.JDUT)
# convert datetime and Ls2
def datetime_to_Ls2(pdt):
"""Convert datetime to Ls2."""
return SCET_to_Ls2(datetime_to_SCET(pdt))
datetime_to_Ls2 = np.frompyfunc(datetime_to_Ls2, nin=1, nout=1)
def datetime_to_SYLs(pdt):
"""Convert datetime to (SY, Ls)."""
return Ls2_to_SYLs(SCET_to_Ls2(datetime_to_SCET(pdt)))
datetime_to_SYLs = np.frompyfunc(datetime_to_SYLs, nin=1, nout=2)
def datetime_to_Ls(pdt):
"""Convert datetime to Ls and SY."""
sy, ls = datetime_to_SYLs(pdt)
return ls
datetime_to_Ls = np.frompyfunc(datetime_to_Ls, nin=1, nout=1)
def Ls2_to_datetime(Ls2):
"""Convert Ls2 to datetime."""
return SCET_to_datetime(Ls2_to_SCET(Ls2))
def SYLs_to_datetime(ls, sy=3, dfa=dfa):
"""Convert SY and Ls to daetime."""
return SCET_to_datetime(Ls2_to_SCET(SYLs_to_Ls2(ls, sy)))
# convert date string to Ls
def datestr_to_Ls2(datestr):
"""Convert datestring to Ls2."""
return datetime_to_Ls2(pd.to_datetime(datestr))
def datestr_to_Ls(datestr, include_SY=False):
"""Convert datestring to Ls.
Returns Ls alone, or (SY, Ls) if include_SY is True."""
sy, ls = Ls2_to_SYLs(datestr_to_Ls2(datestr))
if include_SY:
return sy, ls
return ls
#
# horizons@ssd.jpl.nasa.gov
# the funcions in this file use output from NASA Horizons using the following
'''
!$$SOF
COMMAND= '699'
CENTER= '500@10'
MAKE_EPHEM= 'YES'
TABLE_TYPE= 'OBSERVER'
START_TIME= '1950-01-01'
STOP_TIME= '2070-01-01'
STEP_SIZE= '1 d'
CAL_FORMAT= 'BOTH'
TIME_DIGITS= 'FRACSEC'
ANG_FORMAT= 'DEG'
OUT_UNITS= 'KM-S'
RANGE_UNITS= 'KM'
APPARENT= 'AIRLESS'
SUPPRESS_RANGE_RATE= 'NO'
SKIP_DAYLT= 'NO'
EXTRA_PREC= 'NO'
R_T_S_ONLY= 'NO'
REF_SYSTEM= 'J2000'
CSV_FORMAT= 'YES'
OBJ_DATA= 'YES'
QUANTITIES= '21,44'
!$$EOF
'''
if __name__ == '__main__':
import argparse, re
from datetime import datetime
now_date = datetime.now()
parser = argparse.ArgumentParser()#exit_on_error=False)
parser.add_argument('-S', '--saturn-year', nargs=1, help="provide SY argument")
parser.add_argument('-j', '--julian', help="use Julian date", action='store_true')
parser.add_argument('-s', '--simple', help='simple output', action='store_true')
parser.add_argument('-x', '--extended', help='extended date range (1890 to 2160)', action='store_true')
parser.add_argument('date', nargs='*', help="earth date or Ls", default=[now_date], metavar='date_or_Ls')
try:
args = parser.parse_args()
except argparse.ArgumentError:
parser.print_help()
if args.extended:
dfa.data = load_csv_ephem(ephem_csv_long)
for date in args.date:
outpt = ''
# check if date is formatted like Ls
if re.match(r'\A[+-]*\d+[.]*\d*\Z', date):
if args.saturn_year is not None:
SY = int(args.saturn_year[0])
else:
SY, _ = datestr_to_Ls(now_date, include_SY=True)
Ls = date
date = SYLs_to_datetime(ls=float(Ls), sy=float(SY))
if args.julian:
date = f'{date.to_julian_date():f}'
outpt = f"{date}"
if not args.simple:
outpt = f" Date(SY {int(SY):d}, {float(Ls):g}°) = {date}"
else:
if args.saturn_year is not None:
print("Option -S only applies when calculating Earth dates.")
parser.exit()
if args.julian:
date = pd.to_datetime(date, unit='D', origin='julian')
SY, Ls = datestr_to_Ls(date, include_SY=True)
outpt = f"{int(SY):d} {Ls:g}"
if not args.simple:
outpt = f" Ls({date}) = SY {int(SY):d}, {Ls:g}°"
print(outpt)