-
Notifications
You must be signed in to change notification settings - Fork 0
/
eclipsefinder.py
255 lines (223 loc) · 10.9 KB
/
eclipsefinder.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
from __future__ import print_function
import numpy as np
from astropy.time import Time
from astropy.time import TimeDelta
import matplotlib.pyplot as plt
import pandas as pd
'''
Python program to calculate upcoming EB eclipse times and durations.
Useful for figuring out when to observe, or not.
Originally written by Meredith Rawls, March-April 2014
INPUT:
A text file with six columns:
ID - whatever you want to call your EB
P - orbital period
T0 - zeropoint (presumably in BJD)
pwidth - length of primary eclipse in phase units
swidth - length of secondary eclipse in phase units
sep - separation between two eclipses in phase units
RA - 10-character RA value 00:00:00.0
Dec - 10-character Dec value 00:00:00.0
OUTPUT:
A graph of when eclipses will occur.
Eclipse timing data in two text files: one that is human-friendly and one that is
for use as a control file by the 1m telescope at APO. (You will need to manually edit
the desired exposure times, etc. in this program if you care about the latter.)
Read the damn comments, Jean.
'''
# Put your desired start and end dates in UTC here!!!
startdate = '2015-10-08 00:00'
enddate = '2015-12-31 00:00'
startshort = startdate[0:10]
endshort = enddate[0:10]
infile = '../../1m_observations/RGEB_info_alpha.txt' # infile you must provide
outfile1 = '../../1m_observations/eclipses_2015_OctDec.txt' # human-friendly outfile
outfile2 = '../../1m_observations/RGEB_1m_2015_OctDec.inp' # 1m telescope friendly outfile
# Function to make a 2D list without stabbing your eye out
def make2dList(rows, cols):
a=[]
for row in xrange(rows): a += [[0]*cols]
return a
# Create a list of possible observation times
obs_tstart = Time(startdate, format='iso', scale='utc')
obs_tend = Time(enddate, format='iso', scale='utc')
obs_deltat = TimeDelta(0.5, format='jd')
# arbitrary fraction-of-a-day time resolution... 0.5 is probably sufficient for eclipse
# durations ~1 day (?)
# make sure to set obs_deltat < Porb for all systems or bad things might happen
# Read in data from file
kic, p, t0, pwid, swid, sep, RA, Dec = np.loadtxt(infile, comments='#', usecols=(0,1,2,3,4,5,6,7),
dtype={'names': ('kic', 'p', 't0', 'pwid', 'swid', 'sep', 'RA', 'Dec'),
'formats': (np.float64,np.float64,np.float64,np.float64,np.float64,np.float64,'|S11','|S11')}, unpack=True)
print ('Reading data from ' + infile)
# Create astropy Time objects for the zeropoints and orbital periods
bjd0 = Time(t0+2400000.0, format='jd', scale='utc') # time of primary eclipse
porb = TimeDelta(p, format='jd') # duration of one orbit
tobs = Time(np.arange(obs_tstart.jd, obs_tend.jd, obs_deltat.jd), format='jd', scale='utc')
print ('Forecasting eclipses for this time range:')
print (obs_tstart.iso + ' to ' + obs_tend.iso)
print (obs_tstart.jd, obs_tend.jd)
print (' ')
# Access one element of this observation time array in JD format:
#print tobs[0].utc.jd
# What we'd really like is a list of times when eclipses BEGIN and END that fall
# between obs_tstart and obs_tend. Want to work in real time, not phase-time.
# --> astropy Time values DO NOT play nice with numpy arrays!!!! Use lists instead.
# This gives WAY MORE eclipse times than you could ever want, assuming Porb > obs_deltat
newbjd0_float = []
newbjd0 = []
pri_eclipse_mid = make2dList(len(kic), len(tobs))
sec_eclipse_mid = make2dList(len(kic), len(tobs))
pri_eclipse_start = make2dList(len(kic), len(tobs))
pri_eclipse_end = make2dList(len(kic), len(tobs))
sec_eclipse_start = make2dList(len(kic), len(tobs))
sec_eclipse_end = make2dList(len(kic), len(tobs))
for j in range(0,len(kic)):
# Find the most recent bjd0 time right BEFORE the observation window of interest
newbjd0_float.append( np.floor((obs_tstart.jd - bjd0[j].jd)/porb[j].value) * porb[j].value + bjd0[j].jd )
newbjd0.append( Time(newbjd0_float[j], format='jd', scale='utc') )
i = 0
for i in range(0,len(tobs)):
# Save eclipse midpoints
pri_eclipse_mid[j][i] = newbjd0[j] + i*porb[j]
sec_eclipse_mid[j][i] = newbjd0[j] + i*porb[j] + sep[j]*porb[j]
# print j, i, pri_eclipse_mid[j][i].iso # IT WORKS!
# Save primary eclipse start & end times
pri_eclipse_start[j][i] = pri_eclipse_mid[j][i] - pwid[j]*porb[j]/2
pri_eclipse_end[j][i] = pri_eclipse_mid[j][i] + pwid[j]*porb[j]/2
# Save secondary eclipse start & end times
sec_eclipse_start[j][i] = sec_eclipse_mid[j][i] - swid[j]*porb[j]/2
sec_eclipse_end[j][i] = sec_eclipse_mid[j][i] + swid[j]*porb[j]/2
# The 1-m telescope 'inp' file needs lots of other info too. Here we go.
## FIRSTLINE INFO ##
Epoch = '\t 2000.0'
Xmin = '\t 1.01'
Xmax = '\t 2.00'
#Year = '\t 2007'
Blank = '\t 24.0' #yeah I don't know either
UTstart = '\t 0.00'
UTend = '\t 24.00'
HAmin = '\t -12.00'
HAmax = '\t 12.00'
MPhase = '\t 1.00'
MDist = '\t 5.00'
Nexp = '\t -1'
DT = '\t 3600' # for file1, make this ~ 1 hour = 3600 sec, and for file2, make this 0
Acquire = '\t -1 \t 0 \t 3 \t 0 \t 0'
# -1 is find guiding star & focus. -2 is just find a guiding star. -3 is just go.
# the trailing 0 3 0 0 are for something with guiding
## SECONDLINE INFO ##
#seq_repeat = '1' #this is calculated below instead
# NOTE exposure times manually correspond with the infile order (sorry)
nU = '0'
tU = '0'
nB = 1
tB = [240,120,120,120,75,75,240,90,240,210,75,10,10,90,5,240,75,75,10]
nV = 1
tV = [120,60,60,60,30,30,120,60,120,120,30,5,5,60,5,120,30,30,5]
nR = 1
tR = [45,30,30,30,15,15,45,30,45,60,15,5,5,30,5,60,15,15,5]
nI = 1
tI = [60,45,45,45,20,20,60,45,60,90,20,5,5,45,5,90,20,20,5]
end_line2 = '0 \t 0.0 \t 0 \t 0.0 \t 0 \t 0.0 \t 0 \t 0.0 \t 0 \t 0.0'
# Plot the eclipses as a function of time
# And write useful info to text files!
f1 = open(outfile1, 'w')
f2 = open(outfile2, 'w')
plt.yticks(range(1,len(kic)+1), ['%.0f' % a for a in kic])
tplotstart = obs_tstart
tplotend = obs_tend
plt.axis([tplotstart.plot_date, tplotend.plot_date, 0, len(kic)+1])
# Shade regions in plot corresponding to approximately A half and B half
# (this part written by Jean)
datelist = pd.date_range(start=startshort, end=endshort, freq='6H')
for k in range(0,len(datelist)):
if k % 4 == 0:
plt.bar(datelist[k], height=25, width=.25, color='silver', edgecolor='silver')
# Height is a function of how many stars are being plotted
# Width is 0.25 because it's a quarter of a day (6 hours)
if k % 4 == 1:
plt.bar(datelist[k], height=25, width=.25, color='grey', edgecolor='grey')
# epic loop for plotting AND writing info to files
for j in range(0,len(kic)):
# Calculate how many times to repeat the sequence for ~30 min (1800 sec)
seq_repeat = int(1800 / (tB[j]+20 + tV[j]+20 + tR[j]+20 + tI[j]+20) )
i = 0
while (pri_eclipse_start[j][i] < tplotend or sec_eclipse_start[j][i] < tplotend):
# the while loop below should work, but doesn't, so as it stands (above) there are
# some extra points that occur before tplotstart. sigh.
#while ( (pri_eclipse_start[j][i] < tplotend and pri_eclipse_start[j][i] > tplotstart) or (sec_eclipse_start[j][i] < tplotend and sec_eclipse_start[j][i] > tplotstart) ):
# Put stuff on the graph
# plt.plot([x1,x2],[y1,y2], 'colors-n-stuff') #syntax to draw some lines
plt.plot_date([pri_eclipse_start[j][i].plot_date, pri_eclipse_end[j][i].plot_date],[j+1, j+1], 'k-', lw=2)#, label='primary' if j==0 else '')
plt.plot_date([sec_eclipse_start[j][i].plot_date, sec_eclipse_end[j][i].plot_date],[j+1, j+1], 'r-', lw=2)#, label='secondary' if j==0 else '')
# Write just eclipse timing info to a file
print(int(kic[j]), '\t p \t', pri_eclipse_start[j][i].isot, '\t',
pri_eclipse_mid[j][i].isot, '\t', pri_eclipse_end[j][i].isot, file=f1)
print(int(kic[j]), '\t s \t', sec_eclipse_start[j][i].isot, '\t',
sec_eclipse_mid[j][i].isot, '\t', sec_eclipse_end[j][i].isot, file=f1)
# Define MS, DS, ME, DE (Month/Day Start/End) windows from eclipse times
Yearpri = pri_eclipse_start[j][i].iso[0:4]
MSpri = pri_eclipse_start[j][i].iso[5:7]
DSpri = pri_eclipse_start[j][i].iso[8:10]
MEpri = pri_eclipse_end[j][i].iso[5:7]
DEpri = pri_eclipse_end[j][i].iso[8:10]
Yearsec = sec_eclipse_start[j][i].iso[0:4]
MSsec = sec_eclipse_start[j][i].iso[5:7]
DSsec = sec_eclipse_start[j][i].iso[8:10]
MEsec = sec_eclipse_end[j][i].iso[5:7]
DEsec = sec_eclipse_end[j][i].iso[8:10]
# Create new RGEB.inp file for 1m observations
# LINE 1, primary eclipse
print(int(kic[j]), '\t \t', RA[j], '\t', Dec[j], Epoch, Xmin, Xmax, Yearpri, '\t',
MSpri, '\t', DSpri, '\t', MEpri, '\t', DEpri, '\t', Blank, UTstart, UTend,
HAmin, HAmax, MPhase, MDist, Nexp, DT, Acquire, file=f2)
# LINE 2, primary eclipse
print('\t', seq_repeat, '\t', nU, '\t', tU, '\t', nB, '\t', tB[j], '\t', nV, '\t',
tV[j], '\t', nR, '\t', tR[j], '\t', nI, '\t', tI[j], '\t', end_line2, file=f2)
# LINE 1, secondary eclipse
print(int(kic[j]), '\t \t', RA[j], '\t', Dec[j], Epoch, Xmin, Xmax, Yearsec, '\t',
MSsec, '\t', DSsec, '\t', MEsec, '\t', DEsec, '\t', Blank, UTstart, UTend,
HAmin, HAmax, MPhase, MDist, Nexp, DT, Acquire, file=f2)
# LINE 2, secondary eclipse
print('\t', seq_repeat, '\t', nU, '\t', tU, '\t', nB, '\t', tB[j], '\t', nV, '\t',
tV[j], '\t', nR, '\t', tR[j], '\t', nI, '\t', tI[j], '\t', end_line2, file=f2)
i = i + 1
# print each system one final time (with no date restrictions, and seq_repeat = 5)
print ('# The following is for file 3! Each system is printed once with no time restrictions', file=f2)
print ('# ', file=f2)
for j in range(0,len(kic)):
# LINE 1, no date restrictions
print(int(kic[j]), '\t \t', RA[j], '\t', Dec[j], Epoch, Xmin, Xmax, '2014', '\t',
'0', '\t', '0', '\t', '12', '\t', '31', '\t', Blank, UTstart, UTend,
HAmin, HAmax, MPhase, MDist, '\t 1', DT, Acquire, file=f2)
# LINE 2, no date restrictions
print('\t', '5', '\t', nU, '\t', tU, '\t', nB, '\t', tB[j], '\t', nV, '\t',
tV[j], '\t', nR, '\t', tR[j], '\t', nI, '\t', tI[j], '\t', end_line2, file=f2)
# finish the plot
#plt.xlabel("Time (BJD)")
#plt.ylabel("System")
plt.gcf().autofmt_xdate() #for slanty dates
f1.close()
#plt.legend()
## option to manually plot horizontal lines at dates of interest ##
# windowstart = Time('2015-04-07T09:16:00', format='isot', scale='utc')
# windowend = Time('2015-04-07T12:35:00', format='isot', scale='utc')
# windowstart2 = Time('2015-04-16T09:16:00', format='isot', scale='utc')
# windowend2 = Time('2015-04-16T12:24:00', format='isot', scale='utc')
# windowstart3 = Time('2015-04-27T07:16:00', format='isot', scale='utc')
# windowend3 = Time('2015-04-27T12:12:00', format='isot', scale='utc')
# windowstart4 = Time('2015-05-06T07:16:00', format='isot', scale='utc')
# windowend4 = Time('2015-05-06T12:03:00', format='isot', scale='utc')
# plt.axvline(windowstart.plot_date)
# plt.axvline(windowend.plot_date, ls=':')
# plt.axvline(windowstart2.plot_date)
# plt.axvline(windowend2.plot_date, ls=':')
# plt.axvline(windowstart3.plot_date)
# plt.axvline(windowend3.plot_date, ls=':')
# plt.axvline(windowstart4.plot_date)
# plt.axvline(windowend4.plot_date, ls=':')
plt.title('(black primary, red secondary)')
plt.show()
# primary eclipses are plotted in black
# secondary eclipses are plotted in red