generated from erosita/uds
220 lines
7.0 KiB
Python
220 lines
7.0 KiB
Python
#!/usr/bin/env python
|
|
|
|
__author__ = "Roman Krivonos"
|
|
__copyright__ = "Space Research Institute (IKI)"
|
|
|
|
""" Работает только в архиве """
|
|
|
|
import re
|
|
import sys
|
|
from glob import glob
|
|
from astropy.io import fits
|
|
|
|
from astropy.io.fits import getheader
|
|
|
|
|
|
from astropy.coordinates import SkyCoord # High-level coordinates
|
|
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
|
|
from astropy.coordinates import Angle, Latitude, Longitude # Angles
|
|
import astropy.units as u
|
|
|
|
|
|
def file_grep(pattern=None, filename=None, filt=None):
|
|
fl = open(filename, "r")
|
|
for line in fl:
|
|
if re.search(pattern, line):
|
|
if(filt):
|
|
if(filt in line):
|
|
return(line.replace("\n", ""))
|
|
else:
|
|
return(line.replace("\n", ""))
|
|
return(None)
|
|
|
|
bands={
|
|
'E01':'25.00 60.00',
|
|
'E02':'17.00 22.10',
|
|
'E03':'22.10 28.73',
|
|
'E04':'28.73 37.35',
|
|
'E05':'37.35 48.55',
|
|
'E06':'48.55 63.12',
|
|
'E07':'63.12 82.06',
|
|
'E08':'82.06 106.67',
|
|
'E09':'106.67 138.67',
|
|
'E10':'138.67 180.28',
|
|
'E11':'180.28 234.36',
|
|
'E12':'234.36 304.67',
|
|
}
|
|
|
|
crab_crd = SkyCoord(83.6287, 22.0147, frame=FK5(), unit="deg")
|
|
|
|
key = sys.argv[1]
|
|
|
|
do_mjd = False
|
|
|
|
revs = sorted(glob("????", recursive = False))
|
|
|
|
erange = bands[key]
|
|
|
|
fout="detcnts.{}.fits".format(key)
|
|
|
|
obs_id=[]
|
|
ra=[]
|
|
dec=[]
|
|
rota=[]
|
|
lon=[]
|
|
lat=[]
|
|
t_start=[]
|
|
t_stop=[]
|
|
exposure=[]
|
|
ph1=[]
|
|
ph2=[]
|
|
ph=[]
|
|
|
|
raw=[]
|
|
src=[]
|
|
clean=[]
|
|
npix=[]
|
|
rms=[]
|
|
orbit=[]
|
|
|
|
tstart=[]
|
|
tstop=[]
|
|
mjd=[]
|
|
crab_sep=[]
|
|
|
|
for rev in revs:
|
|
logs = sorted(glob(rev+"/????????????.log", recursive = False))
|
|
|
|
for log in logs:
|
|
res = file_grep(pattern="rev_phase1",filename=log)
|
|
phase1 = float(res.split()[1]) if(res) else 0.0
|
|
|
|
res = file_grep(pattern="rev_phase2",filename=log)
|
|
phase2 = float(res.split()[1]) if(res) else 0.0
|
|
|
|
|
|
res = file_grep(pattern="Final FSD RMS",filename=log)
|
|
rms0 = float(res.split()[-1]) if(res) else None
|
|
|
|
res = file_grep(pattern="OBTMISDC",filename=log)
|
|
obtmisdc = float(res.split()[-1]) if(res) else None
|
|
|
|
|
|
res = file_grep(pattern="DETCNTS2",filename=log,filt=erange)
|
|
d = res.split() if(res) else None
|
|
|
|
if not (d):
|
|
continue
|
|
# ['DETCNTS2', '011000130010', '00001', '266.570', '-26.994', '89.359', '2639.38', '28148860.0', '28152310.0', '3.55128E-03', '8.53612E-04', '2.69767E-03', '1.27125E+04', '28.73', '37.35']
|
|
|
|
|
|
|
|
obsid=d[1]
|
|
ra0=float(d[3])
|
|
dec0=float(d[4])
|
|
rota0=float(d[5])
|
|
texp=float(d[6])
|
|
t1=float(d[7])
|
|
t2=float(d[8])
|
|
raw0=float(d[9])
|
|
src0=float(d[10])
|
|
clean0=float(d[11])
|
|
npix0=float(d[12])
|
|
e1=float(d[13])
|
|
e2=float(d[14])
|
|
|
|
if not (raw0 > 0.0):
|
|
continue
|
|
|
|
obs_id.append(obsid)
|
|
orbit.append(int(obsid[:4]))
|
|
ra.append(ra0)
|
|
dec.append(dec0)
|
|
rota.append(rota0)
|
|
|
|
sc = SkyCoord(ra0, dec0, frame=FK5(), unit="deg")
|
|
lon.append(sc.galactic.l.value if(sc.galactic.l.value<180) else sc.galactic.l.value-360)
|
|
lat.append(sc.galactic.b.value)
|
|
|
|
crab_sep.append(crab_crd.separation(sc).deg)
|
|
|
|
t_start.append(t1)
|
|
t_stop.append(t2)
|
|
exposure.append(texp)
|
|
|
|
raw.append(raw0)
|
|
src.append(src0)
|
|
clean.append(clean0)
|
|
npix.append(npix0)
|
|
|
|
ph1.append(phase1)
|
|
ph2.append(phase2)
|
|
ph.append((phase1+phase2)/2.0)
|
|
rms.append(rms0)
|
|
|
|
if(do_mjd):
|
|
# ~/results1/osa11/scw/0101/010100040010.001/isgri_cor_events.fits.gz
|
|
corfile="/afs/mpa/project/integral/data1/scw/{}/{}.001/swg.fits".format(obsid[:4],obsid)
|
|
d = getheader(corfile, 1)
|
|
mjdref=float(d['MJDREF'])
|
|
# tjd2/86400+obtmisdc+mjdref-40000
|
|
tj1=float(d['TSTART'])
|
|
tj2=float(d['TSTOP'])
|
|
|
|
tstart.append(tj1)
|
|
tstop.append(tj2)
|
|
mjd.append((tj1+tj2)/2.0+mjdref)
|
|
#print(obtmisdc,mjdref)
|
|
#print("{} {} {} {}".format(log,phase1,phase2,obsid))
|
|
else:
|
|
mjd.append(0.0)
|
|
|
|
indices = sorted(
|
|
range(len(t_start)),
|
|
key=lambda index: t_start[index]
|
|
)
|
|
|
|
coldefs = fits.ColDefs([
|
|
fits.Column(name='REV', format='J', array=[orbit[index] for index in indices]),
|
|
fits.Column(name='OBSID', format='12A', array=[obs_id[index] for index in indices]),
|
|
fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]),
|
|
fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]),
|
|
fits.Column(name='ROTA', format='D', unit='deg', array=[rota[index] for index in indices]),
|
|
fits.Column(name='LON', format='D', unit='deg', array=[lon[index] for index in indices]),
|
|
fits.Column(name='LAT', format='D', unit='deg', array=[lat[index] for index in indices]),
|
|
#fits.Column(name='SC_TSTART', format='D', unit='s', array=[t_start[index] for index in indices]),
|
|
#fits.Column(name='SC_TSTOP', format='D', unit='s', array=[t_stop[index] for index in indices]),
|
|
#fits.Column(name='TSTART', format='D', unit='day', array=[tstart[index] for index in indices]),
|
|
#fits.Column(name='TSTOP', format='D', unit='day', array=[tstop[index] for index in indices]),
|
|
fits.Column(name='MJD', format='D', unit='day', array=[mjd[index] for index in indices]),
|
|
fits.Column(name='EXPOSURE', format='D', unit='s', array=[exposure[index] for index in indices]),
|
|
#fits.Column(name='PHASE1', format='D', unit='', array=[ph1[index] for index in indices]),
|
|
#fits.Column(name='PHASE2', format='D', unit='', array=[ph2[index] for index in indices]),
|
|
fits.Column(name='PHASE', format='D', unit='', array=[ph[index] for index in indices]),
|
|
#fits.Column(name='NROWS', format='J', unit='', array=[nrows[index] for index in indices]),
|
|
#fits.Column(name='LIVETIME', format='D', unit='s', array=[livetime[index] for index in indices]),
|
|
fits.Column(name='RAW', format='D', unit='cts/s', array=[raw[index] for index in indices]),
|
|
fits.Column(name='SRC', format='D', unit='', array=[src[index] for index in indices]),
|
|
fits.Column(name='CLEAN', format='D', unit='', array=[clean[index] for index in indices]),
|
|
fits.Column(name='RMS', format='D', unit='', array=[rms[index] for index in indices]),
|
|
fits.Column(name='CRAB_SEP', format='D', unit='deg', array=[crab_sep[index] for index in indices]),
|
|
|
|
])
|
|
|
|
hdu = fits.BinTableHDU.from_columns(coldefs, name='RATE')
|
|
hdu.header['MISSION'] = ('INTEGRAL', '')
|
|
hdu.header['ERANGE'] = (key, '')
|
|
hdu.header['E1'] = (erange.split()[0], '')
|
|
hdu.header['E2'] = (erange.split()[1], '')
|
|
hdu.header['TELESCOP'] = (key, '')
|
|
|
|
hdu.header['INSTITUT'] = ('IKI', 'Affiliation')
|
|
hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person')
|
|
hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail')
|
|
#hdu.add_checksum()
|
|
print(hdu.columns)
|
|
hdu.writeto(fout, overwrite=True)
|
|
|
|
with fits.open(fout, mode='update') as hdus:
|
|
hdus[1].add_checksum()
|