Skip to content

Commit 394f235

Browse files
CharlesBellCharlesBell
authored
Update plotCOBSmagnitudes.py
Updated with more code comments. Added Nakano brightness model. Updated JPL Horizons, MPC, and Yoshida brightness parameters.
1 parent 41d0d5d commit 394f235

File tree

1 file changed

+209
-21
lines changed

1 file changed

+209
-21
lines changed

‎plotCOBSmagnitudes.py

+209-21
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,72 @@
11
# -*- coding: utf-8 -*-
22
"""
3-
Created on Tue Mar 17 08:17:08 2020
43
5-
@author: charb
4+
plotCometC2019Y4MagEstimates
5+
6+
- queries JPL horizons with astroquery
7+
- save the quey as a fits table
8+
- reads the query results from fits table
9+
- extracts columns pertaining to comets
10+
- converts date strings to Time objects
11+
for use with matplotlib.dates as mdates
12+
- extract and plots JPL Horizons Tmag and Nmag values
13+
- calculates and plots magnitudes for MPC, Yoshida, Nakano, and COBS brightness models
14+
- plots COBS observation data
15+
- plots Minor Planet Center observation data by magnitude band
16+
17+
Updated on 2020-Mar-23
18+
@author: Charles Bell
619
"""
720

8-
# plot COBS magnitudes
21+
#
922

10-
import numpy as np
23+
#horizons
24+
#https://astroquery.readthedocs.io/en/latest/
25+
from astroquery.jplhorizons import Horizons
26+
from astropy.table import Table
1127
import matplotlib.pyplot as plt
1228
from astropy.time import Time
13-
import datetime
29+
import math
1430
import matplotlib.dates as mdates
1531
from pandas.plotting import register_matplotlib_converters
1632
register_matplotlib_converters()
1733

34+
def getDateTimeUTHorizons(datestr):
35+
#2021-Jan-01 00:00
36+
yyyy = datestr[0:4]
37+
mm = datestr[5:8]
38+
if (mm == 'Jan'):
39+
mm = '01'
40+
if (mm == 'Feb'):
41+
mm = '02'
42+
if (mm == 'Mar'):
43+
mm = '03'
44+
if (mm == 'Apr'):
45+
mm = '04'
46+
if (mm == 'May'):
47+
mm = '05'
48+
if (mm == 'Jun'):
49+
mm = '06'
50+
if (mm == 'Jul'):
51+
mm = '07'
52+
if (mm == 'Aug'):
53+
mm = '08'
54+
if (mm == 'Sep'):
55+
mm = '09'
56+
if (mm == 'Oct'):
57+
mm = '10'
58+
if (mm == 'Nov'):
59+
mm = '11'
60+
if (mm == 'Dec'):
61+
mm = '12'
62+
dd = datestr[9:11]
63+
s = yyyy + '-' + mm + '-' + dd
64+
hrs = datestr[12:14]
65+
mins = datestr[15:17]
66+
s = s + ' ' + hrs + ':' + mins +":00"
67+
return s
1868

19-
def getDateTimeUT(dateOfObservation):
69+
def getDateTimeUTCOBS(dateOfObservation):
2070
# print(dateOfObservation)
2171
yyyy = dateOfObservation[0:4]
2272
mm = dateOfObservation[5:7]
@@ -39,47 +89,185 @@ def getDateTimeUT(dateOfObservation):
3989
s = s + str(secs)
4090
return s
4191

42-
datetimes = []
43-
magnitudes = []
92+
def getDateTimeUTMPC(dateOfObservation):
93+
yyyy = dateOfObservation[0:4]
94+
mm = dateOfObservation[5:7]
95+
dd = dateOfObservation[8:10]
96+
s = yyyy + '-' + mm + '-' + dd
97+
ff = float(dateOfObservation[10: ])
98+
seconds = int(ff * 86400)
99+
hrs = int(seconds/3600)
100+
mins = int((seconds - hrs * 3600)/60)
101+
secs = seconds - hrs * 3600 - mins * 60
102+
s = s + ' '
103+
if (hrs < 10):
104+
s = s + "0"
105+
s = s + str(hrs) + ":"
106+
if (mins < 10):
107+
s = s + "0"
108+
s = s + str(mins) + ":"
109+
if (secs < 10):
110+
s = s + "0"
111+
s = s + str(secs)
112+
return s
113+
114+
cobsdatetimes = []
115+
cobsobsmags = []
44116

45117
# COBS observation data can be found on line at https://cobs.si/
46-
f = open('C2019 Y4 (ATLAS)-COBS-2020Mar17.txt', 'r')
118+
f = open('C2019 Y4 (ATLAS)-COBS-2020Mar22.txt', 'r')
47119
for lines in f:
48120
if (lines.startswith('IIIYYYYMnL')):
49121
pass
50122
else:
51123
dateOfObservation = lines[11:26].strip()
52-
datetimeUT = getDateTimeUT(dateOfObservation)
124+
datetimeUT = getDateTimeUTCOBS(dateOfObservation)
53125
mag = lines[28:32].strip()
54126
# print(dateOfObservation, datetimeUT, mag)
55127
dateTime = Time(datetimeUT, format='iso').datetime
56-
datetimes.append(dateTime)
57-
magnitudes.append(float(mag))
58-
128+
cobsdatetimes.append(dateTime)
129+
cobsobsmags.append(float(mag))
130+
131+
#MPC observations
132+
f2 = open('D:\\Comets and Asteroids\\C2019 Y4 (ATLAS)\\C_2019_Y4.txt', "r")
133+
MPCm1dates = []
134+
MPCm1mags = []
135+
136+
MPCm2dates = []
137+
MPCm2mags = []
138+
139+
MPCRdates = []
140+
MPCRmags = []
141+
142+
MPCVdates = []
143+
MPCVmags = []
144+
145+
MPCGdates = []
146+
MPCGmags = []
147+
148+
for lines in f2:
149+
packed = lines[0:12].strip()
150+
orbitType = lines[4:5].strip()
151+
typeOfObservation = lines[14:15].strip()
152+
year = int(lines[15:19])
153+
# print(lines)
154+
155+
156+
if ((typeOfObservation != 's') and (typeOfObservation != 'S')):
157+
dateOfObservation = lines[15:32].strip()
158+
# print('dateOfObservation', dateOfObservation)
159+
datetimeUT = getDateTimeUTMPC(dateOfObservation)
160+
# print('datetimeUT', datetimeUT)
161+
dt = Time(datetimeUT, format='iso').datetime
162+
# print('dt', dt)
163+
#mpcdates.append(dt)
164+
observedRA = lines[32:44].strip()
165+
observedDecl = lines[44:56].strip()
166+
observedMagnitude = lines[65:70].strip(' ')
167+
#observedMagnitude = lines[65:70]
168+
#print('observedMagnitude', observedMagnitude)
169+
magnitudeBand = lines[70:71].strip()
170+
observatoryCode = lines[77:80].strip();
171+
reference = lines[72:77].strip();
172+
if ((len(observedMagnitude) > 0) and (magnitudeBand == "N")):
173+
MPCm2dates.append(dt)
174+
MPCm2mags.append(float(observedMagnitude))
175+
elif ((len(observedMagnitude) > 0) and (magnitudeBand == "T")):
176+
MPCm1dates.append(dt)
177+
MPCm1mags.append(float(observedMagnitude))
178+
elif ((len(observedMagnitude) > 0) and (magnitudeBand == "R")):
179+
MPCRdates.append(dt)
180+
MPCRmags.append(float(observedMagnitude))
181+
elif ((len(observedMagnitude) > 0) and (magnitudeBand == "V")):
182+
MPCVdates.append(dt)
183+
MPCVmags.append(float(observedMagnitude))
184+
elif ((len(observedMagnitude) > 0) and (magnitudeBand == "G")):
185+
MPCGdates.append(dt)
186+
MPCGmags.append(float(observedMagnitude))
187+
188+
# Geocentric [500]
189+
# obj = Horizons(id='C/2019 Y4', location='500',epochs={'start':'2019-12-28', 'stop':'2021-01-01','step':'1d'} )
190+
# eph = obj.ephemerides()
191+
# eph.write('ephem.fits', format='fits', overwrite=True )
192+
193+
eph = Table.read('ephem.fits')
194+
datetime_str = eph.columns[eph.index_column('datetime_str')]
195+
jd = eph.columns[eph.index_column('datetime_jd')]
196+
r = eph.columns[eph.index_column('r')]
197+
delta = eph.columns[eph.index_column('delta')]
198+
phase = eph.columns[eph.index_column('alpha')]
199+
200+
TP= 2459000.51510898
201+
datetimes = []
202+
jplTmags = []
203+
jplNmags = []
204+
mpcmags = []
205+
yoshidamags = []
206+
cobsdates = []
207+
cobsmags = []
208+
nakanodates = []
209+
nakanomags = []
210+
for datestr in datetime_str:
211+
datetimeUT = getDateTimeUTHorizons(datestr)
212+
dateTime = Time(datetimeUT, format='iso').datetime
213+
datetimes.append(dateTime)
214+
215+
i = 0
216+
while (i < len(datetimes)):
217+
mpcmags.append(6 + 5 * math.log10(delta[i]) + 2.5 * 4 * math.log10(r[i]))
218+
jplTmags.append(7.9 + 5 * math.log10(delta[i]) + 21 * math.log10(r[i]))
219+
jplNmags.append(13.3 + 5 * math.log10(delta[i]) + 5. * math.log10(r[i]) + .030*phase[i])
220+
yoshidamag = -0.5 + 5 * math.log10(delta[i]) + 36 * math.log10(r[i])
221+
if (jd[i] > TP-73):
222+
yoshidamag = 5.5 + 5 * math.log10(delta[i]) + 10 * math.log10(r[i])
223+
yoshidamags.append(yoshidamag)
224+
if ((jd[i] < (TP - 70))):
225+
cobsdates.append(datetimes[i])
226+
cobsmags.append(-1.62 + 5 * math.log10(delta[i]) + 2.5 * 16.94 * math.log10(r[i]))
227+
if ((jd[i] > (TP-30)) and (jd[i] < (TP+31))):
228+
nakanodates.append(datetimes[i])
229+
nakanomags.append(5.5 + 5 * math.log10(delta[i]) + 10 * math.log10(r[i]))
230+
231+
i = i + 1
232+
59233
font = {'family': 'serif',
60234
'color': 'black',
61235
'weight': 'normal',
62-
'size': 16
236+
'size': 16,
63237
}
64238

65239
boxprops = dict(boxstyle='round', facecolor='#ffffff', ec='#CCCCCC', alpha=0.9)
66-
credittext = 'Credit: COBS Comet Observation Database \n CC BY-NA-SA 4.0 https://cobs.si/'
240+
241+
observers = 'Observers: Christian Harder, Timo Karhula, Maciej Kwinta, \nPiotr Guzik, Maciej Reszelski, Nick James, Denis Buczynski, \nPiotr Nowak, Maik Meyer, Sandor Szabo, Kevin Hills,\n Teerasak Thaluang, Carl Hergenrother, Carlos Labordena, \nMariusz Swietnicki, Juan Jose Gonzalez Suarez, Jerzy Bohusz, \nArtyom Novichonok, Steffen Fritsche, Uwe Pilz,Alex Scholten, \nThomas Lehmann, Jose Pablo Navarro Pina, David Swan, \nJohan Warell, Nirmal Paul, Pieter-Jan Dekelver,\nMartin Masek, Gideon van Buitenen'
242+
credittext = 'Credit: COBS Comet Observation Database – CC BY-NA-SA 4.0\n'
67243

68244
fig, ax = plt.subplots(figsize=(20, 13))
245+
69246
locator = mdates.AutoDateLocator()
70247
formatter = mdates.AutoDateFormatter(locator)
71248
ax.xaxis.set_major_locator(locator)
72249
ax.xaxis.set_major_formatter(formatter)
73250

74251
ax.grid(b=True, which='major', axis='both')
75-
ax.set_title('COBS Magnitude Data Comet C/2019 Y4 (ATLAS)', fontdict=font)
252+
ax.set_title('Comet C/2019 Y4 (ATLAS) Magnitude Estimates (2020-Mar-22) ', fontsize=24)
76253

77-
if (len(magnitudes) > 0):
78-
ax.plot_date(datetimes, magnitudes, label='C/2019 Y4 magnitude', xdate=True, ydate=False)
79-
plt.text('2020-02-15', 17, credittext, fontdict=font, bbox=boxprops)
254+
# plt.text('2020-07-15', -5, observers, fontdict=font, bbox=boxprops)
255+
# plt.text('2020-04-01', 24, credittext, fontdict=font, bbox=boxprops)
80256

81-
plt.xlabel('Date')
257+
ax.plot_date(datetimes, mpcmags, '.', label='MPC M1 mag H0=6 n=4', xdate=True, ydate = False)
258+
ax.plot_date(datetimes, jplTmags,'.', label='JPL Horizons T-mag H0=7.7 n=5.4', xdate=True, ydate = False)
259+
ax.plot_date(datetimes, jplNmags, '.',label='JPL Horizons N-mag H0=13.3 n=1.25 phcoeff=0.030', xdate=True, ydate = False)
260+
ax.plot_date(datetimes, yoshidamags,'.', label='S. Yoshida Total mag H0=-0.5 n=14.4; H0=5.5 n=4', xdate=True, ydate = False)
261+
ax.plot_date(cobsdatetimes, cobsobsmags, '*', label='COBS observed mag', xdate=True, ydate=False)
262+
ax.plot_date(cobsdates, cobsmags, '.', label='COBS mag H0=-1.62 n=16.94', xdate=True, ydate = False)
263+
ax.plot_date(nakanodates, nakanomags, '.', label='NAKANO mag H0=-5.5 n=4 ', xdate=True, ydate = False)
264+
ax.plot_date(MPCm1dates, MPCm1mags, '^', label='MPC observed m1 mag', xdate=True, ydate = False)
265+
ax.plot_date(MPCm2dates, MPCm2mags, 'v', label='MPC observed m2 mag', xdate=True, ydate = False)
266+
ax.plot_date(MPCRdates, MPCRmags, '.', label='MPC observed R mag', xdate=True, ydate = False)
267+
ax.plot_date(MPCVdates, MPCVmags, '.', label='MPC observed V mag', xdate=True, ydate = False)
268+
ax.plot_date(MPCGdates, MPCGmags, '.', label='MPC observed G mag', xdate=True, ydate = False)
269+
plt.xlabel('Date', fontsize=20)
82270
plt.gca().invert_yaxis()
83-
plt.ylabel('Magnitude')
271+
plt.ylabel('Magnitude', fontsize=20)
84272
ax.legend(loc=1)
85273
plt.show()

0 commit comments

Comments
 (0)