-
Notifications
You must be signed in to change notification settings - Fork 1
/
find_sun.py
107 lines (84 loc) · 3.79 KB
/
find_sun.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
#!/usr/bin/env python
import logging
import numpy
import sys
from pyrap.tables import table
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.coordinates import solar_system_ephemeris, EarthLocation
from astropy.coordinates import get_body_barycentric, get_body, get_moon
def rad2deg(xx):
return 180.0*xx/numpy.pi
def get_fields(myms):
fieldtab = table(myms+'/FIELD',ack=False)
ids = fieldtab.getcol("SOURCE_ID")
names = fieldtab.getcol("NAME")
dirs = fieldtab.getcol("PHASE_DIR")
return ids,names,dirs
def match_field(ids,names,dirs,i):
idx = numpy.where(ids == i)[0]
radec = rad2deg(dirs[idx,0,:])[0]
return names[i],radec[0],radec[1]
def calcsep(ra0,dec0,ra1,dec1):
c0 = SkyCoord(ra0*u.deg,dec0*u.deg,frame='fk5')
c1 = SkyCoord(ra1*u.deg,dec1*u.deg,frame='fk5')
sep = round(c0.separation(c1).deg,4)
return sep
def format_coords(ra0,dec0):
c = SkyCoord(ra0*u.deg,dec0*u.deg,frame='fk5')
hms = str(c.ra.to_string(u.hour))
dms = str(c.dec)
return hms,dms
def main():
# MeerKAT
obs_lat = -30.71323598930457
obs_lon = 21.443001467965008
loc = EarthLocation.from_geodetic(obs_lat,obs_lon) #,obs_height,ellipsoid)
myms = sys.argv[1].rstrip('/')
maintab = table(myms,ack=False)
scans = list(numpy.unique(maintab.getcol('SCAN_NUMBER')))
ids,names,dirs = get_fields(myms)
logfile = 'sun_'+myms+'.log'
logging.basicConfig(filename=logfile, level=logging.DEBUG, format='%(asctime)s | %(message)s', datefmt='%d/%m/%Y %H:%M:%S ')
stream = logging.StreamHandler()
stream.setLevel(logging.DEBUG)
streamformat = logging.Formatter('%(asctime)s | %(message)s', datefmt='%d/%m/%Y %H:%M:%S ')
stream.setFormatter(streamformat)
mylogger = logging.getLogger(__name__)
mylogger.setLevel(logging.DEBUG)
mylogger.addHandler(stream)
mylogger.info(myms+' | '+str(len(ids))+' fields | '+str(len(scans))+' scans')
#header = 'Scan Field ID t[iso] t[s] t0[s] t1[s] int0 int1 Duration[m] N_int'
header = 't[iso] Scan Field Name SunRA[deg] SunDec[deg] SunRA[hms] SunDec[dms] SunSep[deg] MoonRA[deg] MoonDec[deg] MoonRA[hms] MoonDec[dms] MoonSep[deg]'
mylogger.info('-'*len(header))
mylogger.info(header)
mylogger.info('-'*len(header))
for scan in scans:
subtab = maintab.query(query='SCAN_NUMBER=='+str(scan))
field = numpy.unique(subtab.getcol('FIELD_ID'))[0]
name,field_ra,field_dec = match_field(ids,names,dirs,field)
field_hms,field_dms = format_coords(field_ra,field_dec)
t_scan = numpy.mean(subtab.getcol('TIME'))
t = Time(t_scan/86400.0,format='mjd')
with solar_system_ephemeris.set('builtin'):
sun = get_body('Sun', t, loc)
moon = get_body('Moon', t, loc)
sun_ra = sun.ra.value
sun_dec = sun.dec.value
moon_ra = moon.ra.value
moon_dec = moon.dec.value
sun_hms,sun_dms = format_coords(sun_ra,sun_dec)
moon_hms,moon_dms = format_coords(moon_ra,moon_dec)
delta_ra_sun = field_ra - sun_ra
delta_dec_sun = field_dec - sun_dec
delta_ra_moon = field_ra - moon_ra
delta_dec_moon = field_dec - moon_dec
sun_sep = calcsep(field_ra,field_dec,sun_ra,sun_dec)
moon_sep = calcsep(field_ra,field_dec,moon_ra,moon_dec)
mylogger.info('%-28s %-5i %-5i %-12s %-12f %-12f %-16s %-16s %-12f %-12f %-12f %-16s %-16s %-12f' %
(t.iso,scan,field,name,sun_ra,sun_dec,sun_hms,sun_dms,sun_sep,moon_ra,moon_dec,moon_hms,moon_dms,moon_sep))
mylogger.info('-'*len(header))
if __name__ == "__main__":
main()