Note

This tutorial was generated from a Jupyter notebook that can be downloaded here.

Writing PSRFITS

import numpy as np
import pdat

Instantiate a draft file from a template

template = '../../../templates/search_template.fits'
new_psrfits = 'my_psrfits.fits'

psrfits1=pdat.psrfits(new_psrfits,from_template=template)

Making new SEARCH mode PSRFITS file using template from path:
'../../../templates/search_template.fits'.
Writing to path 'my_psrfits.fits'.
The Binary Table HDU headers will be written as they are added
     to the PSRFITS file.

The pdat package is very helpful for making PSRFITS file from drafts. Many of the parameters in the PRIMARY and SUBINT HDUs are codependent on one another, and the program keeps track of these dependencies for the user. When you instantiate a psrfits with a template a set of draft headers is made from the template so that you can edit them before writing to disk. This “template-to-write” scheme exists because many pieces of software that will analyze PSRFITS will baulk at files without a complete set of header information.

The template fits file is accessible as fits_template.

def set_primary_header(psrfits_object,prim_dict):
    """
    prim_dict = dictionary of primary header changes
    """
    PF_obj = psrfits_object
    for key in prim_dict.keys():
        PF_obj.replace_FITS_Record('PRIMARY',key,prim_dict[key])

def set_subint_header(psrfits_object,subint_dict):
    """
    prim_dict = dictionary of primary header changes
    """
    PF_obj = psrfits_object
    for key in subint_dict.keys():
        PF_obj.replace_FITS_Record('SUBINT',key,subint_dict[key])

def make_subint_BinTable(self):
    subint_draft = self.make_HDU_rec_array(self.nsubint, self.subint_dtype)
    return subint_draft

psrfits1.fits_template[0].read_header()['OBSERVER']

'GALILEOGALILEI'

The draft headers are editable, and can be changed until you write the file to disk. The ImageHDU that conatins the primary header is names ‘PRIMARY’. The others all go by the name of the EXTNAME. [‘PRIMARY’,‘SUBINT’,‘POLYCO’,‘HISTORY’,‘PARAM’]

psrfits1.draft_hdrs['SUBINT']

XTENSION= 'BINTABLE'           / * Subintegration data  *
BITPIX  =                    8 / N/A
NAXIS   =                    2 / 2-dimensional binary table
NAXIS1  =             33636428 / width of table in bytes
NAXIS2  =                    4 / Number of rows in table (NSUBINT)
PCOUNT  =                    0 / size of special data area
GCOUNT  =                    1 / one data group (required keyword)
TFIELDS =                   17 / Number of fields per row
TTYPE1  = 'TSUBINT '           / Length of subintegration
TFORM1  = '1D      '           / Double
TTYPE2  = 'OFFS_SUB'           / Offset from Start of subint centre
TFORM2  = '1D      '           / Double
TTYPE3  = 'LST_SUB '           / LST at subint centre
TFORM3  = '1D      '           / Double
TTYPE4  = 'RA_SUB  '           / RA (J2000) at subint centre
TFORM4  = '1D      '           / Double
TTYPE5  = 'DEC_SUB '           / Dec (J2000) at subint centre
TFORM5  = '1D      '           / Double
TTYPE6  = 'GLON_SUB'           / [deg] Gal longitude at subint centre
TFORM6  = '1D      '           / Double
TTYPE7  = 'GLAT_SUB'           / [deg] Gal latitude at subint centre
TFORM7  = '1D      '           / Double
TTYPE8  = 'FD_ANG  '           / [deg] Feed angle at subint centre
TFORM8  = '1E      '           / Float
TTYPE9  = 'POS_ANG '           / [deg] Position angle of feed at subint centre
TFORM9  = '1E      '           / Float
TTYPE10 = 'PAR_ANG '           / [deg] Parallactic angle at subint centre
TFORM10 = '1E      '           / Float
TTYPE11 = 'TEL_AZ  '           / [deg] Telescope azimuth at subint centre
TFORM11 = '1E      '           / Float
TTYPE12 = 'TEL_ZEN '           / [deg] Telescope zenith angle at subint centre
TFORM12 = '1E      '           / Float
TTYPE13 = 'DAT_FREQ'           / [MHz] Centre frequency for each channel
TFORM13 = '2048E   '           / NCHAN floats
TTYPE14 = 'DAT_WTS '           / Weights for each channel
TFORM14 = '2048E   '           / NCHAN floats
TTYPE15 = 'DAT_OFFS'           / Data offset for each channel
TFORM15 = '8192E   '           / NCHAN*NPOL floats
TTYPE16 = 'DAT_SCL '           / Data scale factor for each channel
TFORM16 = '8192E   '           / NCHAN*NPOL floats
TTYPE17 = 'DATA    '           / Subint data table
TFORM17 = '33554432B'          / NBIN*NCHAN*NPOL*NSBLK int, byte(B) or bit(X)
INT_TYPE= 'TIME    '           / Time axis (TIME, BINPHSPERI, BINLNGASC, etc)
INT_UNIT= 'SEC     '           / Unit of time axis (SEC, PHS (0-1), DEG)
SCALE   = 'FluxDen '           / Intensity units (FluxDen/RefFlux/Jansky)
NPOL    =                    4 / Nr of polarisations
POL_TYPE= 'IQUV    '           / Polarisation identifier (e.g., AABBCRCI, AA+BB)
TBIN    = 0.000770559999999999 / [s] Time per bin or sample
NBIN    =                    1 / Nr of bins (PSR/CAL mode; else 1)
NBIN_PRD=                    0 / Nr of bins/pulse period (for gated data)
PHS_OFFS=                   0. / Phase offset of bin 0 for gated data
NBITS   =                    8 / Nr of bits/datum (SEARCH mode 'X' data, else 1)
NSUBOFFS=                    0 / Subint offset (Contiguous SEARCH-mode files)
NCHAN   =                 2048 / Number of channels/sub-bands in this file
CHAN_BW =            -0.390625 / [MHz] Channel/sub-band width
NCHNOFFS=                    0 / Channel/sub-band offset for split files
NSBLK   =                 4096 / Samples/row (SEARCH mode, else 1)
EXTNAME = 'SUBINT  '           / name of this binary table extension
TUNIT1  = 's       '           / Units of field
TUNIT2  = 's       '           / Units of field
TUNIT3  = 's       '           / Units of field
TUNIT4  = 'deg     '           / Units of field
TUNIT5  = 'deg     '           / Units of field
TUNIT6  = 'deg     '           / Units of field
TUNIT7  = 'deg     '           / Units of field
TUNIT8  = 'deg     '           / Units of field
TUNIT9  = 'deg     '           / Units of field
TUNIT10 = 'deg     '           / Units of field
TUNIT11 = 'deg     '           / Units of field
TUNIT12 = 'deg     '           / Units of field
TUNIT13 = 'MHz     '           / Units of field
TDIM17  = '(1,2048,4,4096)'    / Dimensions (NBITS or NBIN,NCHAN,NPOL,NSBLK)
TUNIT17 = 'Jy      '           / Units of subint data
EXTVER  =                    1 / auto assigned by template parser

In order to set the dimensions of the data arrays within the SUBINT HDU there is a convenience function called set_subint_dims. By setting the dimensions using this function the dependencies on these dimensions, inclduing memory allocation, will be propagated through the headers correctly.

First lets choose some dimensions for the data.

sample_size = 20.48e-3 # in milliseconds
ROWS = 30
N_Time_Bins = 2048*ROWS
Total_time = round(N_Time_Bins*sample_size)
dt = Total_time/N_Time_Bins
subband =1.5625
BW=200
N_freq = int(BW/subband)
Npols = 4
print('Total_time',Total_time/1e3)
print('N_freq',N_freq)

Total_time 1.258
N_freq 128

And then call the set_subint_dims method.

psrfits1.set_subint_dims(nsblk=2048,nchan=N_freq,nsubint=ROWS,npol=Npols)

Once we have set the SUBINT dimensions a subint_dtype list is made which we can then use to make a recarray to hold the data. Here nsubint is the same as above, and has been made an attribute.

subint_draft = psrfits1.make_HDU_rec_array(psrfits1.nsubint, psrfits1.subint_dtype)

All of the header cards can bet set by assigning them to the appropriate member of the draft header.

npol = psrfits1.draft_hdrs['SUBINT']['NPOL']

Here we set the time per subintegration (time length of an NSBLK) and the offsets, which are the times at the center of each subintegration from the beginning of the observation.

tsubint = data.shape[-1]*dt*1e-3 #in seconds
offs_sub_init = tsubint/2
offs_sub = np.zeros((ROWS))

for jj in range(ROWS):
    offs_sub[jj] = offs_sub_init + (jj * tsubint)

Here we just use the values from the template file.

lst_sub = psrfits1.fits_template[1]['LST_SUB'].read()[0]
ra_sub = psrfits1.fits_template[1]['RA_SUB'].read()[0]
dec_sub = psrfits1.fits_template[1]['DEC_SUB'].read()[0]
glon_sub = psrfits1.fits_template[1]['GLON_SUB'].read()[0]
glat_sub = psrfits1.fits_template[1]['GLAT_SUB'].read()[0]
fd_ang = psrfits1.fits_template[1]['FD_ANG'].read()[0]
pos_ang = psrfits1.fits_template[1]['POS_ANG'].read()[0]
par_ang = psrfits1.fits_template[1]['PAR_ANG'].read()[0]
tel_az = psrfits1.fits_template[1]['TEL_AZ'].read()[0]
tel_zen = psrfits1.fits_template[1]['TEL_ZEN'].read()[0]

ones = np.ones((ROWS))
#And assign them using arrays of the appropriate sizes
subint_draft['TSUBINT'] = tsubint * ones
subint_draft['OFFS_SUB'] = offs_sub
subint_draft['LST_SUB'] = lst_sub * ones
subint_draft['RA_SUB'] = ra_sub * ones
subint_draft['DEC_SUB'] = dec_sub * ones
subint_draft['GLON_SUB'] = glon_sub * ones
subint_draft['GLAT_SUB'] = glat_sub * ones
subint_draft['FD_ANG'] = fd_ang * ones
subint_draft['POS_ANG'] = pos_ang * ones
subint_draft['PAR_ANG'] = par_ang * ones
subint_draft['TEL_AZ'] = tel_az * ones
subint_draft['TEL_ZEN'] = tel_zen * ones

Here we’ll just make some data of the correct shape.

data = np.random.randn(ROWS,1,N_freq,Npols,2048)

And now we can assign the data arrays

for ii in range(subint_draft.size):
    subint_draft[ii]['DATA'] = data[ii,:,:,:,:]
    subint_draft[ii]['DAT_SCL'] = np.ones(N_freq*npol)
    subint_draft[ii]['DAT_OFFS'] = np.zeros(N_freq*npol)
    subint_draft[ii]['DAT_FREQ'] = np.linspace(1300,1500,N_freq)
    subint_draft[ii]['DAT_WTS'] = np.ones(N_freq)

subint_hdr=psrfits1.draft_hdrs['SUBINT']

from decimal import *
getcontext().prec=12
a=Decimal(S1.TimeBinSize*1e-3)
a.to_eng_string()

'0.00002047526041666666474943582498813299253015429712831974029541015625'

b='{0:1.18f}'.format(Decimal(a.to_eng_string()))
b

'0.000020475260416667'

pri_dic= {'OBSERVER':'GALILEOGALILEI','OBSFREQ':S1.f0,'OBSBW':S1.bw,'OBSNCHAN':S1.Nf}
subint_dic = {'TBIN':b,'CHAN_BW':S1.freqBinSize}

subint_dic['TBIN']

'0.000020475260416667'

psrfits1.make_FITS_card(subint_hdr,'TBIN',subint_dic['TBIN'])

{'card_string': 'TBIN    = 0.000020475260416667 / [s] Time per bin or sample',
 'class': 150,
 'comment': '[s] Time per bin or sample',
 'dtype': 'F',
 'name': 'TBIN',
 'value': 2.0475260416667e-05,
 'value_orig': 2.0475260416667e-05}

psrfits1.draft_hdrs['SUBINT'].records()[65]

{'card_string': "TUNIT8  = 'deg     '           / Units of field",
 'class': 70,
 'comment': 'Units of field',
 'dtype': 'C',
 'name': 'TUNIT8',
 'value': 'deg     ',
 'value_orig': 'deg     '}

set_primary_header(psrfits1,pri_dic)

set_subint_header(psrfits1,subint_dic)

psrfits1.draft_hdrs['SUBINT'].records()[47]['value'] = '0.000020483398437500'
psrfits1.draft_hdrs['SUBINT'].records()[47]['value_orig'] = '0.000020483398437500'
psrfits1.draft_hdrs['SUBINT'].records()[47]

{'card_string': 'TBIN    = 0.000020475260416667 / [s] Time per bin or sample',
 'class': 150,
 'comment': '[s] Time per bin or sample',
 'dtype': 'F',
 'name': 'TBIN',
 'value': '0.000020483398437500',
 'value_orig': '0.000020483398437500'}

psrfits1.HDU_drafts['SUBINT'] = subint_draft

psrfits1.write_psrfits()

psrfits1.close()