tofil.py 4.99 KB
Newer Older
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
1
2
3
4
5
6
7
#!/usr/bin/env python2
from __future__ import print_function

import numpy as np
import filterbank
import sigproc
from astropy.time import Time
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
8
from astropy import units as u
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
9
10

import os.path
11
import os
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
12
13
14
15
16
import sys
import re

import argparse

17
18
19
20
21
22
23
telescope_found = False
try:
    from telescope import telescope
    telescope_found = True
except:
    print("Could not connect to console, no RA/DEC recorded")

24
def create_filterbank(outfile, source, infile=None, date=None, headeronly=False, centerfrequency=420., time=30):
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
25
26
27
28
    """
    Create a filterbank file from the CAMRAS-format, which is generated
    by pulsar_record.

29
    outfile:  name of the filterbank file to be created
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
30
31
32
33
34
35
    infile:   name of the file generated by pulsar_record
    date:     the date to be used; default is the file creation date of infile
              if given, the date must be in a format the astropy.time understands
    source:   name of the source to be used; default is guessed from filename
    """

36
    source = "PSR " + source
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
37
38

    if not date:
39
40
41
42
        if infile:
            date = Time(os.path.getctime(infile), format='unix')
        else:
            date = Time.now()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
43
44
45
    else:
        date = Time(date, format='isot')

46
47
48
49
50
    az_start = 0.
    za_start = 0.
    src_raj = 0.
    src_dej = 0.
    if telescope_found:
51
        dt = telescope()
52
53
        az_start = (180*u.deg + dt.getAzEl(waitForUpdate=True)[0]).to(u.deg).value
        za_start = (90*u.deg - dt.getAzEl()[1]).to(u.deg).value
54
55
56
        src_raj = dt.getRaDec().ra.to_string(unit=u.hour, sep='')
        src_dej = dt.getRaDec().dec.to_string(unit=u.degree, sep='')

57
58
    print("Input file: {}".format(infile))
    print("Output file: {}".format(outfile))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
59
    print("Source name: {}".format(source))
60
    print("Observation time: {}".format(date.isot))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
61
62

    date = date.mjd
63
64
65
66
67

    if infile:
        rawdatafile = infile
    else:
        rawdatafile = ""
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
68
69
70
71
72
73
 
    # Generate filterbank header
    fil_header = {}
    fil_header["telescope_id"]  = sigproc.telescope_ids["Effelsberg"]
    fil_header["machine_id"]    = sigproc.machine_ids["FAKE"]
    fil_header["data_type"]     = 1
74
    fil_header["rawdatafile"]   = rawdatafile
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
75
76
77
    fil_header["source_name"]   = source
    fil_header["barycentric"]   = 0
    fil_header["pulsarcentric"] = 0
78
79
80
81
    fil_header["az_start"]      = az_start
    fil_header["za_start"]      = za_start
    fil_header["src_raj"]       = src_raj
    fil_header["src_dej"]       = src_dej
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
82
83
    fil_header["tstart"]        = date
    fil_header["tsamp"]         = 0.00046811428571414528 # 1.0/35e6*256 chans*64 decimation
84
    fil_header["nbits"]         = 16
85
    fil_header["fch1"]          = centerfrequency+21.4
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
86
87
88
89
90
    fil_header["foff"]          = -35.0/256.0
    fil_header["nchans"]        = 256
    fil_header["nifs"]          = 1

    # Write header
91
    out = filterbank.create_filterbank_file(outfile, fil_header, nbits=16)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
92

93
94
95
    if headeronly:
        print("Wrote filterbank header only")
        out.close()
96
97
        return

98
99
    if not infile:
        out.close()
100
101
102
        cmd = "./pulsar_filterbank " + str(time) + " " + outfile
        print("Running "+cmd)
        os.system(cmd)
103
104
        return

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
105
    # Read file
106
    print("Converting "+infile)
107
    f = open(infile)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
108

109
110
    filesize = os.stat(infile).st_size
    offset = 0
111
    CHUNKSIZE = 30000
112
    while offset < filesize:
113
        # Loop through file in chunks of 256*10000 numbers
114
        # i.e. 256*CHUNKSIZE*4 bytes (10Mb) to avoid memory problems
115
        f.seek(offset, os.SEEK_SET)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
116

117
        data = np.fliplr(np.fromfile(f, dtype='>u4', count=256*CHUNKSIZE).reshape(-1,256))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
118

119
        # Check that no packets were dropped
120
121
122
        if np.sum(data[1:,-1] - data[:-1,-1] - 1) > 0.1:
            print("Packets were dropped:", np.count_nonzero(data[:-1,-1] - data[1:,-1] - 1),
                                           np.sum(data[1:,-1] - data[:-1,-1] - 1))
123
124

        # Write data
125
        out.append_spectra(data.astype('u2'))
126

127
        offset += 256*CHUNKSIZE*4
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
128
129
130
131

    out.close()

if __name__ == "__main__":
132
    parser = argparse.ArgumentParser(description="Create a filterbank file from the output of the CAMRAS backend in pulsar mode.")
133
    parser.add_argument("outfile", help="Name of output file, e.g. 'B0329.fil'")
134
    parser.add_argument("sourcename", help="Name of the source, e.g. 'B0329+54'")
135
136
137
138
139
    parser.add_argument("-i", "--infile", help="Input CAMRAS pulsar file. If empty, it will start a new recording.")
    parser.add_argument("-d", "--date", help="Start date/time of observation, in 'isot' format (defaults to now, or creation date of infile)", default=None)
    parser.add_argument("--header", help="Only record a header", default=False, action="store_true")
    parser.add_argument("-f", "--frequency", help="Center frequency in MHz, default: 420.", default=420., type=float)
    parser.add_argument("-t", "--time", help="Time for observing in seconds, default: 300s (5min)", default=300)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
140
141
142

    args = parser.parse_args()

143
    create_filterbank(args.outfile, args.sourcename, infile=args.infile, date=args.date, headeronly=args.header, centerfrequency=args.frequency, time=args.time)