-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathPlotJPLHorizonsEphemeris.py
157 lines (135 loc) · 5.36 KB
/
PlotJPLHorizonsEphemeris.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
# -*- coding: utf-8 -*-
"""
PlotJPLHorizonsEphemeris
- queries JPL horizons with astroquery
- save the quey as a fits table
- reads the query results from fits table
- extracts columns pertaining to comets
- converts date strings to Time objects
for use with matplotlib.dates as mdates
- plots each selected column in a subplot
- mdates share x-axis on each sub plot
Created on Sat Mar 21 07:41:24 2020
@author: Charles Bell
"""
import matplotlib.pyplot as plt
#https://astroquery.readthedocs.io/en/latest/
# from astroquery.jplhorizons import Horizons
from astropy.table import Table
from astropy.time import Time
import matplotlib.dates as mdates
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
"""
2021-Jan-01 00:00
returns 2021-01-01 00:00:00.000
which can be input to astropy.time.Time iso
to plot dates with matplotlib
"""
def getDateTimeUTHorizons(datestr):
yyyy = datestr[0:4]
mm = datestr[5:8]
if (mm == 'Jan'):
mm = '01'
if (mm == 'Feb'):
mm = '02'
if (mm == 'Mar'):
mm = '03'
if (mm == 'Apr'):
mm = '04'
if (mm == 'May'):
mm = '05'
if (mm == 'Jun'):
mm = '06'
if (mm == 'Jul'):
mm = '07'
if (mm == 'Aug'):
mm = '08'
if (mm == 'Sep'):
mm = '09'
if (mm == 'Oct'):
mm = '10'
if (mm == 'Nov'):
mm = '11'
if (mm == 'Dec'):
mm = '12'
dd = datestr[9:11]
s = yyyy + '-' + mm + '-' + dd
hrs = datestr[12:14]
mins = datestr[15:17]
s = s + ' ' + hrs + ':' + mins +":00.000"
return s
# uncomment the next three lines to query horizons and save the ephemeris to a fits table
# obj = Horizons(id='C/2019 Y4', location='500',epochs={'start':'2020-01-01', 'stop':'2021-01-01','step':'1d'} )
# eph = obj.ephemerides()
# eph.write('ephem.fits', format='fits', overwrite=True )
# assumes the ephemeris fits table has been saved
# if you are modifying this demo code, you should only query the horisons server once
# and then use the file copy to modify to your needs.
eph = Table.read('ephem.fits')
objectname = eph.columns[eph.index_column('targetname')]
datetimestr = eph.columns[eph.index_column('datetime_str')]
jd = eph.columns[eph.index_column('datetime_jd')]
r = eph.columns[eph.index_column('r')]
delta = eph.columns[eph.index_column('delta')]
tmag = eph.columns[eph.index_column('Tmag')]
nmag = eph.columns[eph.index_column('Nmag')]
phase = eph.columns[eph.index_column('alpha')]
elong = eph.columns[eph.index_column('elong')]
psang = eph.columns[eph.index_column('sunTargetPA')]
psamv = eph.columns[eph.index_column('velocityPA')]
plang = eph.columns[eph.index_column('OrbPlaneAng')]
trueanomaly = eph.columns[eph.index_column('true_anom')]
ra = eph.columns[eph.index_column('RA')]
dec = eph.columns[eph.index_column('DEC')]
constellation = eph.columns[eph.index_column('constellation')]
datetimes = []
for datestr in datetimestr:
datetimeUT = getDateTimeUTHorizons(datestr)
dateTime = Time(datetimeUT, format='iso').datetime
datetimes.append(dateTime)
fig, axs = plt.subplots(9, sharex=True, figsize=(20, 18))
fig.suptitle('JPL Horizons Ephemeris C/2019 Y4 (ATLAS)', x=0.50, y=0.91, fontsize=24)
locator = mdates.AutoDateLocator()
formatter = mdates.AutoDateFormatter(locator)
axs[0].xaxis.set_major_locator(locator)
axs[0].xaxis.set_major_formatter(formatter)
axs[1].xaxis.set_major_locator(locator)
axs[1].xaxis.set_major_formatter(formatter)
axs[2].xaxis.set_major_locator(locator)
axs[2].xaxis.set_major_formatter(formatter)
axs[3].xaxis.set_major_locator(locator)
axs[3].xaxis.set_major_formatter(formatter)
axs[4].xaxis.set_major_locator(locator)
axs[4].xaxis.set_major_formatter(formatter)
axs[5].xaxis.set_major_locator(locator)
axs[5].xaxis.set_major_formatter(formatter)
axs[6].xaxis.set_major_locator(locator)
axs[6].xaxis.set_major_formatter(formatter)
axs[7].xaxis.set_major_locator(locator)
axs[7].xaxis.set_major_formatter(formatter)
axs[8].xaxis.set_major_locator(locator)
axs[8].xaxis.set_major_formatter(formatter)
# plot magnitude plots with inverted y-axis
axs[0].invert_yaxis() # Total magnitude
axs[1].invert_yaxis() # Nuclear magnitude
# The x-axis formatted dates are shared in all plots
axs[8].set_xlabel('Date', fontsize=20)
axs[0].set_ylabel('Tmag', fontsize=20)
axs[1].set_ylabel('Nmag', fontsize=20)
axs[2].set_ylabel('r', fontsize=20)
axs[3].set_ylabel('delta',fontsize=20)
axs[4].set_ylabel('Phase',fontsize=20)
axs[5].set_ylabel('Elong',fontsize=20)
axs[6].set_ylabel('PsAng',fontsize=20)
axs[7].set_ylabel('PsAMV',fontsize=20)
axs[8].set_ylabel('PlAng',fontsize=20)
axs[0].plot_date(datetimes, tmag, '.', xdate=True, ydate = False)
axs[1].plot_date(datetimes, nmag, '.', xdate=True, ydate = False)
axs[2].plot_date(datetimes, r, '.', xdate=True, ydate = False)
axs[3].plot_date(datetimes, delta, '.', xdate=True, ydate = False)
axs[4].plot_date(datetimes, phase, '.', xdate=True, ydate = False)
axs[5].plot_date(datetimes, elong, '.', xdate=True, ydate = False)
axs[6].plot_date(datetimes, psang, '.', xdate=True, ydate = False)
axs[7].plot_date(datetimes, psamv, '.', xdate=True, ydate = False)
axs[8].plot_date(datetimes, plang, '.', xdate=True, ydate = False)