tofil.py 6.46 KB
Newer Older
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
1
#!/usr/bin/env python3
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
2
3
4
5
6

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

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
9
10
11
12
13
14
15
try:
    from camrasdevices import Receiver
    have_hpib = True
except ImportError:
    print("Could not connect to HPIB, skipping that")
    have_hpib = False

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
16
import os.path
17
import os
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
18
19
import sys
import re
20
import subprocess
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
21
22
23

import argparse

24
25
26
27
28
29
30
telescope_found = False
try:
    from telescope import telescope
    telescope_found = True
except:
    print("Could not connect to console, no RA/DEC recorded")

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
31
def create_filterbank(source, infile=None, outfile=None, date=None, headeronly=False, centerfrequency=420., duration=30, showprogressbar=False):
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
32
33
34
35
36
    """
    Create a filterbank file from the CAMRAS-format, which is generated
    by pulsar_record.

    infile:   name of the file generated by pulsar_record
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
37
    outfile:  name of the filterbank file to be created
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
38
39
40
41
42
    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
    """

43
    source = "PSR " + source
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
44

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

56
57
    if not date:
        if infile:
58
            raise("When specifying an infile, also specify a date, e.g. '2018-01-30T11:06:30.651'")
59
60
        else:
            date = Time.now()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
61
            print("Waiting for a time multiple of 10 seconds...")
62
63
            while date.datetime.second % 10 != 0:
                date = Time.now()
64
65
66
    else:
        date = Time(date, format='isot')

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
67
68
69
    if not outfile:
        outfile = date.iso.replace(' ','-')[0:19]+'_'+source[4:]+'.fil'

70
71
    print("Input file: {}".format(infile))
    print("Output file: {}".format(outfile))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
72
    print("Source name: {}".format(source))
73
    print("Observation time: {}".format(date.isot))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
74
    print("Center frequency: {}".format(centerfrequency))
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
75
76

    date = date.mjd
77
78
79
80
81

    if infile:
        rawdatafile = infile
    else:
        rawdatafile = ""
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
82
83
84
 
    # Generate filterbank header
    fil_header = {}
85
    fil_header["telescope_id"]  = sigproc.telescope_ids["Dwingeloo"]
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
86
87
    fil_header["machine_id"]    = sigproc.machine_ids["FAKE"]
    fil_header["data_type"]     = 1
88
    fil_header["rawdatafile"]   = rawdatafile
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
89
90
91
    fil_header["source_name"]   = source
    fil_header["barycentric"]   = 0
    fil_header["pulsarcentric"] = 0
92
93
94
95
    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
96
97
    fil_header["tstart"]        = date
    fil_header["tsamp"]         = 0.00046811428571414528 # 1.0/35e6*256 chans*64 decimation
98
    fil_header["nbits"]         = 16
99
    fil_header["fch1"]          = centerfrequency+21.4
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
100
101
102
103
104
    fil_header["foff"]          = -35.0/256.0
    fil_header["nchans"]        = 256
    fil_header["nifs"]          = 1

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

107
108
109
    if headeronly:
        print("Wrote filterbank header only")
        out.close()
110
111
        return

112
113
    if not infile:
        out.close()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
114
        cmd = "/home/dijkema/pulsar/pulsar_filterbank " + str(duration) + " " + outfile
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
115
        #cmd = "sleep "+str(duration)
116
        print("Running "+cmd)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
117
118
        measure_process = subprocess.Popen(cmd.split())

119
120
121
122
123
124
125
126
127
128
129
130
131
132
        try:
            if showprogressbar:
                from tqdm import tqdm
                import time
                time.sleep(1) # Give measure process time to print something on the screen

                with tqdm(total=duration) as pbar:
                    for time_chunk in range(duration):
                        pbar.update(1)
                        if measure_process.poll() is not None:
                            break
                        time.sleep(1)
        finally:
            measure_process.wait()
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
133
134
135
        if measure_process.returncode != 0:
            raise RuntimeError("Measure process exited abnormally")

136
        print("\a")
137
138
        return

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
139
    # Read file
140
    print("Converting "+infile)
141
    f = open(infile)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
142

143
144
    filesize = os.stat(infile).st_size
    offset = 0
145
    CHUNKSIZE = 30000
146
    while offset < filesize:
147
        # Loop through file in chunks of 256*10000 numbers
148
        # i.e. 256*CHUNKSIZE*4 bytes (10Mb) to avoid memory problems
149
        f.seek(offset, os.SEEK_SET)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
150

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

153
        # Check that no packets were dropped
154
155
156
        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))
157
158

        # Write data
159
        out.append_spectra(data.astype('u2'))
160

161
        offset += 256*CHUNKSIZE*4
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
162
163
164
165

    out.close()

if __name__ == "__main__":
166
    parser = argparse.ArgumentParser(description="Create a filterbank file from the output of the CAMRAS backend in pulsar mode.")
167
    parser.add_argument("sourcename", help="Name of the source, e.g. 'B0329+54'")
168
    parser.add_argument("-o", "--outfile", help="Name of output file, e.g. 'B0329.fil'. Default: '2018-01-17-20:38:50_B0329+54.fil'")
169
170
171
172
    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)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
173
    parser.add_argument("-t", "--time", help="Duration for observing in seconds, default: 300s (5min)", default=300, type=int)
Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
174
175
176

    args = parser.parse_args()

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
177
178
179
180
    if have_hpib:
        rec = Receiver()
        rec.frequency = args.frequency*1.e6

Tammo Jan Dijkema's avatar
Tammo Jan Dijkema committed
181
    create_filterbank(args.sourcename, infile=args.infile, outfile=args.outfile, date=args.date, headeronly=args.header, centerfrequency=args.frequency, duration=args.time, showprogressbar=True)