Comparing Arecibo, VLA, and Effelsberg constraints on FRB 121102 burst spectra

Build broad band spectra of FRB for Law et al 2017

Note that this notebook assumes:

  • rtpipe <= v1.54
  • sdmpy<=1.35

This is due to an error in the DM scale that must be preserved to properly redetect the FRB. Starting at v1.55, the DM scale factor changed from 4.2e-3 to 4.1488e-3.

The DM scale error is 1.2%, so DM=560.5 (in v 1.54) is actually 567.4 (in v1.55 and in other codes and in reality).


In [1]:
%matplotlib inline

In [2]:
import numpy as np
import pylab as pl
import astropy.io.fits as fits

import rtpipe
import rtlib_cython as rtlib

import astropy.units as units
import astropy.coordinates as coord
from astropy.time import Time


2017-07-26 13:25:35,321 - py.warnings - WARNING - /Users/caseyjlaw/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:1405: UserWarning: 
This call to matplotlib.use() has no effect because the backend has already
been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)


In [3]:
# confirm version is is earlier than 1.54 if using old dm scale
print(rtpipe.__version__)


1.53

Useful functions


In [4]:
dmdelay_new = lambda dm, f0, f1: 4.1488e-3*dm*(1./f0**2 - 1./f1**2)  # inconsistent with rtpipe
dmdelay = lambda dm, f0, f1: 4.2e-3*dm*(1./f0**2 - 1./f1**2)  # consistent with rtpipe, but requires scaling by 1.2%

Read coherently dedispersed Arecibo dynamic spectrum


In [5]:
#name = 'puppi_57648_C0531+33_0048_3756.99.ar.paz'
name = 'puppi_57648_C0531+33_0048_3757.00.ar.paz.pT'
fits.info(name)


Filename: puppi_57648_C0531+33_0048_3757.00.ar.paz.pT
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      64   ()      
  1  HISTORY       1 BinTableHDU     71   3R x 28C   [24A, 256A, 8A, 8A, 1J, 1I, 1I, 1I, 1D, 1D, 1J, 1D, 1D, 1D, 1I, 1I, 1I, 1I, 1I, 32A, 32A, 32A, 256A, 32A, 32A, 1I, 32A, 1I]   
  2  PSRPARAM      1 BinTableHDU     12   7R x 1C   [128A]   
  3  POLYCO        1 BinTableHDU     38   1R x 13C   [24A, 16A, 1I, 1I, 1I, 8A, 1D, 1D, 1D, 1D, 1D, 1D, 15D]   
  4  SUBINT        1 BinTableHDU     72   1R x 10C   [1D, 1D, 1D, 1D, 1D, 512D, 512E, 512E, 512E, 2097152I]   

In [6]:
hdu = fits.open(name)
hdu0, hdu1, hdu2, hdu3, hdu4 = hdu[0], hdu[1], hdu[2], hdu[3], hdu[4]

In [7]:
hdu0.header


Out[7]:
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                    8 / number of bits per data pixel                  
NAXIS   =                    0 / number of data axes                            
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
COMMENT   FITS (Flexible Image Transport System) format defined in Astronomy and
COMMENT   Astrophysics Supplement Series v44/p363, v44/p371, v73/p359, v73/p365.
COMMENT   Contact the NASA Science Office of Standards and Technology for the   
COMMENT   FITS Definition document #100 and other FITS information.             
HDRVER  = '5.4             '   / Header version                                 
FITSTYPE= 'PSRFITS         '   / FITS definition for pulsar data files          
DATE    = '2017-03-10T14:23:15' / File creation date (YYYY-MM-DDThh:mm:ss UTC)  
OBSERVER= 'J.Hessels'          / Observer name(s)                               
PROJID  = 'p3054   '           / Project name                                   
TELESCOP= 'Arecibo '           / Telescope name                                 
ANT_X   =            882589.65 / [m] Antenna ITRF X-coordinate (D)              
ANT_Y   =          -4924872.32 / [m] Antenna ITRF Y-coordinate (D)              
ANT_Z   =          3943729.348 / [m] Antenna ITRF Z-coordinate (D)              
FRONTEND= 'lbw     '           / Receiver ID                                    
IBEAM   = '        '           / Beam ID for multibeam systems                  
NRCVR   =                    2 / Number of receiver polarisation channels       
FD_POLN = 'LIN     '           / LIN or CIRC                                    
FD_HAND =                   -1 / +/- 1. +1 is LIN:A=X,B=Y, CIRC:A=L,B=R (I)     
FD_SANG =                  45. / [deg] FA of E vect for equal sig in A&B (E)    
FD_XYPH =                   0. / [deg] Phase of A^* B for injected cal (E)      
BACKEND = 'PUPPI   '           / Backend ID                                     
BECONFIG= 'N/A     '           / Backend configuration file name                
BE_PHASE=                   -1 / 0/+1/-1 BE cross-phase:0 unknown,+/-1 std/rev  
BE_DCC  =                    0 / 0/1 BE downconversion conjugation corrected    
BE_DELAY=                   0. / [s] Backend propn delay from digitiser input   
TCYCLE  =                   0. / [s] On-line cycle time (D)                     
OBS_MODE= 'PSR     '           / (PSR, CAL, SEARCH)                             
DATE-OBS= '2016-09-17T09:26:33.000' / Date of observation (YYYY-MM-DDThh:mm:ss U
OBSFREQ =           1380.78125 / [MHz] Centre frequency for observation         
OBSBW   =                -800. / [MHz] Bandwidth for observation                
OBSNCHAN=                  512 / Number of frequency channels (original)        
CHAN_DM =                 557. / [cm-3 pc] DM used for on-line dedispersion     
PNT_ID  = '        '           / Name or ID for pointing ctr (multibeam feeds)  
SRC_NAME= 'J0531+33'           / Source or scan ID                              
COORD_MD= 'J2000   '           / Coordinate mode (J2000, GALACTIC, ECLIPTIC)    
EQUINOX =                2000. / Equinox of coords (e.g. 2000.0)                
RA      = '05:31:58.600'       / Right ascension (hh:mm:ss.ssss)                
DEC     = '+33:08:49.600'      / Declination (-dd:mm:ss.sss)                    
BMAJ    =                   0. / [deg] Beam major axis length                   
BMIN    =                   0. / [deg] Beam minor axis length                   
BPA     =                   0. / [deg] Beam position angle                      
STT_CRD1= '05:31:58.600'       / Start coord 1 (hh:mm:ss.sss or ddd.ddd)        
STT_CRD2= '+33:08:49.600'      / Start coord 2 (-dd:mm:ss.sss or -dd.ddd)       
TRK_MODE= 'TRACK   '           / Track mode (TRACK, SCANGC, SCANLAT)            
STP_CRD1= '05:31:58.6000'      / Stop coord 1 (hh:mm:ss.sss or ddd.ddd)         
STP_CRD2= '+33:08:49.6003'     / Stop coord 2 (-dd:mm:ss.sss or -dd.ddd)        
SCANLEN = '*       '           / [s] Requested scan length (E)                  
FD_MODE = 'FA      '           / Feed track mode - FA, CPA, SPA, TPA            
FA_REQ  =                   0. / [deg] Feed/Posn angle requested (E)            
CAL_MODE= '                '   / Cal mode (OFF, SYNC, EXT1, EXT2)               
CAL_FREQ= '*       '           / [Hz] Cal modulation frequency (E)              
CAL_DCYC= '*       '           / Cal duty cycle (E)                             
CAL_PHS = '*       '           / Cal phase (wrt start time) (E)                 
CAL_NPHS= '*       '           / Number of states in cal pulse (I)              
STT_IMJD=                57648 / Start MJD (UTC days) (J - long integer)        
STT_SMJD=                37750 / [s] Start time (sec past UTC 00h) (J)          
STT_OFFS=    0.336629511829426 / [s] Start time offset (D)                      
STT_LST =               14053. / [s] Start LST (D)                              

Define python names for Arecibo header info


In [8]:
nch0 = 512
obsfreq = 1780.
df = 1.5625
tint0 = 2.04780080448418e-05
#dm0 = 560.5
dm0 = 553.7  # equivalent to 560.5 in correct scaling

dmt_ao = dmdelay(dm0, obsfreq*1e-3, 1e4)/(24*3600)
bary_ao = 75.7555109 # seconds, from presto
mjd0 = int(57648) # subtract this for convenience. must be an integer mjd
#mjdfrac_ao = (37750 + 0.319370962786796)/(24*3600) - dmt_ao # start time from header
#mjdfrac_ao = (37750 + 0.336629511829426)/(24*3600) - dmt_ao # start time from header, second round
mjdfrac_ao = 57648.43692057008 - dmt_ao - mjd0 + bary_ao/(24*3600)

nint0 = 4096
dt_ao = tint0*nint0

Read dynamic spectrum and average down a bit

Define time/freq binning and barycentric time grid for Arecibo data


In [9]:
tbin = 8
nint1 = nint0/tbin
tint0 = dt_ao/nint0
print('Orig time resolution {0}'.format(tint0))
tint1 = tint0*tbin
print('New time resolution {0}'.format(tint1))

data = fits.getdata(name, ext=4)
spec = data[0][9]  # spectrum read as (npol, nchan, nint) array
specI = spec[:2].mean(axis=0)  # average over first two to get stokes I **confirm**
specI = specI.reshape(nch0, nint1, tbin).mean(axis=2)  # use array slicing to rebin

# flag bad channels
flagch = range(0, 16) + [99, 100, 101, 286, 287, 310, 311] + range(320, 356) + range(nch0-110, nch0)
specI[flagch] = 0
specI = np.ma.masked_equal(specI, 0)

fbin = 2
nch1 = nch0/fbin
print('Orig freq resolution {0}'.format(df))
print('New freq resolution {0}'.format(df*fbin))
specAO = specI.reshape(nch1, fbin, nint1).mean(axis=1)  # use array slicing to rebin
meanspec = specAO.mean(axis=1)
specAO = specAO - meanspec[:,None]

specAO = (specAO-np.ma.mean(specAO))/np.ma.std(specAO)  # convert to snr per pixel
#specAO = specAO/5  # scale by gain relative to VLA?
print(specAO.shape)

tbary_ao = np.linspace(mjdfrac_ao, mjdfrac_ao+(dt_ao/(24*3600)), nint1)
freqs_ao = np.linspace(obsfreq - fbin*df*nch1, obsfreq, nch1)

extent_ao = (tbary_ao[0], tbary_ao[-1], freqs_ao[0], freqs_ao[-1])


Orig time resolution 2.04780080448e-05
New time resolution 0.000163824064359
Orig freq resolution 1.5625
New freq resolution 3.125
(256, 512)

In [10]:
fig = pl.figure(figsize=(12,8))
pl.imshow(specAO, interpolation='nearest', aspect='auto', extent=extent_ao)
pl.colorbar()


Out[10]:
<matplotlib.colorbar.Colorbar at 0x126369f90>

In [11]:
pl.figure(figsize=(8,8))
pl.subplot(211)
pl.plot(tbary_ao, specAO.mean(axis=0))
pl.ylabel('Amp (arb)')
pl.xlabel('Time (s)')
pl.subplot(212)
pl.plot(freqs_ao, meanspec)
pl.xlabel('Freq (MHz)')
pl.ylabel('Amp (arb)')


Out[11]:
<matplotlib.text.Text at 0x124ac8250>

Arecibo burst SNR


In [12]:
sp = specAO.reshape(nch1, nint1/16, 16).mean(axis=2)
print(sp.mean(axis=0).max()/sp.mean(axis=0)[:20].std())


38.87800296

In [13]:
tpk_ao = tbary_ao[np.where(specAO.mean(axis=0) == specAO.mean(axis=0).max())][0]
print('%.10f' % tpk_ao)


0.4377895303

crude integrated flux measurement


In [14]:
sigma = 4
peakbins = np.where(specAO.mean(axis=0) >= sigma*specAO.mean(axis=0)[:256].std())
peaktimes = tbary_ao[peakbins]
window = (peaktimes.max()-peaktimes.min())*24*3600
print('{0} sigma limit selects {1} ms of pulse ({2} bins)'.format(sigma, 
                                                                  window*1e3, 
                                                                  len(peaktimes)))
print('Window width is a bit wider than properly-dedispersed pulse, due to extra sweep')


4 sigma limit selects 2.7904592093 ms of pulse (16 bins)
Window width is a bit wider than properly-dedispersed pulse, due to extra sweep

In [15]:
Sint_sys = specAO.mean(axis=0)[peakbins].mean()
noise = specAO.mean(axis=0)[:peakbins[0][0]].std()/np.sqrt(len(peaktimes))

Sint = (3./np.sqrt(600e6*window*2))*Sint_sys/noise
print('Integrated flux density over {0} ms pulse: {1} mJy'.format(window, Sint*1e3))


Integrated flux density over 0.0027904592093 ms pulse: 57.057635674 mJy

In [16]:
Sspec_sys = specAO[:, peakbins[0]].mean(axis=1)
noise_spec = specAO[:, :peakbins[0][0]].std()/np.sqrt(len(peaktimes))
Sspec = (3./np.sqrt(3.125e6*window*2))*Sspec_sys/noise_spec
print('A (less good) integrated flux density over 2 ms pulse: {0} mJy'.format(int(Sspec.mean()*1e3)))


A (less good) integrated flux density over 2 ms pulse: 62 mJy

Make nice Arecibo burst spectrum like the VLA ones


In [17]:
fig = pl.figure(figsize=(15,7))
ax = fig.add_subplot(111)
pl.plot(1e-3*freqs_ao[::-1], Sspec, 'k.')
pl.text(0.75, 0.88, '57648, Arecibo', horizontalalignment='left', fontsize=24,
        verticalalignment='center', transform=ax.transAxes)
pl.errorbar(1.600, 1.2*Sspec.max(), yerr=(3./np.sqrt(3.125e6*window*2)), fmt='k.', ecolor='k')
pl.ylim(-0.03, Sspec.max()*1.4)
pl.xlabel('Frequency (GHz)', fontsize=18)
pl.ylabel('Flux density (Jy)', fontsize=18)
xt = pl.setp(ax.get_xticklabels(), fontsize=18)
yt = pl.setp(ax.get_yticklabels(), fontsize=18)
ax.xaxis.set_tick_params(width=4, color='k')
ax.yaxis.set_tick_params(width=4, color='k')
fig.savefig('specAO_57648.pdf', format='pdf')


VLA

Set calibration table and set useful functions for rtpipe


In [18]:
calstring = """2.0520    2.89698    0.00279
2.1800    *******    *******
2.3080    *******    *******
2.4360    3.53585    0.00377
2.5640    3.69554    0.00376
2.6920    3.85507    0.00423
2.8200    4.00438    0.00486
2.9480    4.11069    0.00562
3.0520    4.20375    0.00631
3.1800    4.29385    0.00662
3.3080    4.36557    0.00715
3.4360    4.43684    0.00786
3.5640    4.46937    0.00850
3.6920    4.52488    0.00860
3.8200    4.53571    0.00969
3.9480    4.54625    0.00859""" 
# parse flux scale
freq = []
flux = []
eflux = []
for line in calstring.split('\n'): 
    if '*' not in line:
        result = line.split()
        freq.append(float(result[0]))
        flux.append(float(result[1]))
        eflux.append(float(result[2]))
calfreq = np.array(freq)
calflux = np.array(flux)
print(calfreq, calflux)


(array([ 2.052,  2.436,  2.564,  2.692,  2.82 ,  2.948,  3.052,  3.18 ,
        3.308,  3.436,  3.564,  3.692,  3.82 ,  3.948]), array([ 2.89698,  3.53585,  3.69554,  3.85507,  4.00438,  4.11069,
        4.20375,  4.29385,  4.36557,  4.43684,  4.46937,  4.52488,
        4.53571,  4.54625]))

In [19]:
def getscannum(sdmfile):
    sdm = rtpipe.parsesdm.getsdm(sdmfile)
    for scan in sdm.scans():
        try:
            print('Scan {0} binary data file: {1}'.format(scan.idx, scan.bdf.fname))
            bdfscan = int(scan.idx)
        except IOError:
            pass
    return bdfscan

def read_cut(sdmfile, scan, segment, dm=558., dt=1, gainfile=None, **kwargs):
    if not gainfile:
        gainfile = '.'.join(sdmfile.split('.')[:-1] + ['GN'])

    st = rtpipe.RT.set_pipeline(sdmfile, scan, dmarr=[dm], dtarr=[dt], flaglist=[('badap', 3., 0.2)], 
                                uvoversample=1.5, gainfile=gainfile, flagantsol=True, 
                                timesub='mean', logfile=False, savecands=False,
                                savenoise=False, **kwargs)
        
    data = rtpipe.RT.pipeline_reproduce(st, candloc=[segment,0,0,0,0], product='data')
    u, v, w = rtpipe.parsesdm.get_uvw_segment(st, segment)

    return st, data, u, v, w

In [20]:
def correctdata(st, data, u, v, w, corr='ph,dm', lm = (-3.835e-04,5.406e-04)):
    """ lm gives (ra, dec) = (5 31 58.703708986 33 8 52.5067634154)
    as quoted in Chatterjee et al (2017)
    """
    data2 = data.copy()
    
    if 'ph' in corr:
        l1, m1 = lm
        rtlib.phaseshift_threaded(data2, st, l1, m1, u, v)

    if 'dm' in corr:
        rtlib.dedisperse_par(data2, st['freq'], st['inttime'], st['dmarr'][0], [0, st['nbl']])
    
    return data2

# get array2 for bin in array near value
def find_nearest(array, array2, value):
    idx = (np.abs(array-value)).argmin()
    return array2[idx]

def getscale(st):
    # get flux scaling at nearest frequency
    scale = []
    for i in range(len(st['freq'])):
        freq = st['freq'][i]
        scale.append(find_nearest(calfreq, calflux, freq))
#         print(i, st['freq'][i], scale)
    scale = np.array(scale, dtype='complex64')[None,None,:,None]
    return scale

def correct_all(st, data, u, v, w):
    scale = getscale(st)
    dataph = correctdata(st, data*scale, u, v, w, corr='ph')
    dataphdm = correctdata(st, dataph, u, v, w, corr='dm')
    return dataphdm

Read data with rtpipe and phase it to FRB 121102


In [21]:
key = '57648'
read = {}
sdmfile = '16A-496_sb32698778_1_02h00m_000.57648.37452900463.cut/'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 7, npix_max=7400, chans=range(2,256))


Scan 25 binary data file: /Users/caseyjlaw/code/FRB121102_private/16A-496_sb32698778_1_02h00m_000.57648.37452900463.cut/ASDMBinary/uid____evla_bdf_1474107905881
2017-07-26 13:25:37,532 - rtpipe.parsesdm - INFO - Setting (standard) key logfile to False
2017-07-26 13:25:37,533 - rtpipe.parsesdm - INFO - Setting (standard) key timesub to mean
2017-07-26 13:25:37,534 - rtpipe.parsesdm - INFO - Setting (standard) key gainfile to 16A-496_sb32698778_1_02h00m_000.57648.37452900463.GN
2017-07-26 13:25:37,534 - rtpipe.parsesdm - INFO - Setting (standard) key dmarr to [558.0]
2017-07-26 13:25:37,535 - rtpipe.parsesdm - INFO - Setting (standard) key dtarr to [1]
2017-07-26 13:25:37,536 - rtpipe.parsesdm - INFO - Setting (standard) key chans to [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, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255]
2017-07-26 13:25:37,537 - rtpipe.parsesdm - INFO - Setting (standard) key savecands to False
2017-07-26 13:25:37,537 - rtpipe.parsesdm - INFO - Setting (standard) key flaglist to [('badap', 3.0, 0.2)]
2017-07-26 13:25:37,538 - rtpipe.parsesdm - INFO - Setting (standard) key savenoise to False
2017-07-26 13:25:37,539 - rtpipe.parsesdm - INFO - Setting (standard) key flagantsol to True
2017-07-26 13:25:37,539 - rtpipe.parsesdm - INFO - Setting (standard) key npix_max to 7400
2017-07-26 13:25:37,540 - rtpipe.parsesdm - INFO - Setting (standard) key uvoversample to 1.5
2017-07-26 13:25:37,687 - rtpipe.parsesdm - WARNING - No BDF found for scans [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33]
2017-07-26 13:25:38,090 - rtpipe.parsesdm - WARNING - No BDF found for scans [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33]
2017-07-26 13:25:38,229 - rtpipe.parsesdm - INFO - Calculating uvw for first integration of scan 25 of source FRB121102-off
2017-07-26 13:25:38,431 - rtpipe.parsesdm - INFO - 

2017-07-26 13:25:38,432 - rtpipe.parsesdm - INFO - Metadata summary:
2017-07-26 13:25:38,432 - rtpipe.parsesdm - INFO - 	 Working directory and data at /Users/caseyjlaw/code/FRB121102_private, 16A-496_sb32698778_1_02h00m_000.57648.37452900463.cut
2017-07-26 13:25:38,434 - rtpipe.parsesdm - INFO - 	 Using scan 25, source FRB121102-off
2017-07-26 13:25:38,435 - rtpipe.parsesdm - INFO - 	 nants, nbl: 27, 351
2017-07-26 13:25:38,435 - rtpipe.parsesdm - INFO - 	 Freq range (2.497 -- 3.509). 8 spw with 254 chans.
2017-07-26 13:25:38,436 - rtpipe.parsesdm - INFO - 	 Scan has 800 ints (4.0 s) and inttime 0.005 s
2017-07-26 13:25:38,437 - rtpipe.parsesdm - INFO - 	 2 polarizations: ['RR', 'LL']
2017-07-26 13:25:38,437 - rtpipe.parsesdm - INFO - 	 Ideal uvgrid npix=(7776,7776) and res=104 (oversample 1.5)
2017-07-26 13:25:38,466 - rtpipe - INFO - 
2017-07-26 13:25:38,467 - rtpipe - INFO - Pipeline summary:
2017-07-26 13:25:38,468 - rtpipe - INFO - 	 Products saved with 16A-496_sb32698778_1_02h00m_000.57648.37452900463.cut. telcal calibration with 16A-496_sb32698778_1_02h00m_000.57648.37452900463.GN
2017-07-26 13:25:38,469 - rtpipe - INFO - 	 Using 15 segments of 88 ints (0.4 s) with overlap of 0.2 s
2017-07-26 13:25:38,469 - rtpipe - INFO - 		 Lots of segments needed, since Max DM sweep (0.2 s) close to segment size (0.44 s)
2017-07-26 13:25:38,470 - rtpipe - INFO - 	 Downsampling in time/freq by 1/1 and skipping 0 ints from start of scan.
2017-07-26 13:25:38,471 - rtpipe - INFO - 	 Excluding ants []
2017-07-26 13:25:38,471 - rtpipe - INFO - 	 Using pols ['RR', 'LL']
2017-07-26 13:25:38,472 - rtpipe - INFO - 
2017-07-26 13:25:38,473 - rtpipe - INFO - 	 Search with image1 and threshold 7.0.
2017-07-26 13:25:38,474 - rtpipe - INFO - 	 Using 1 DMs from 558.0 to 558.0 and dts [1].
2017-07-26 13:25:38,474 - rtpipe - INFO - 	 Using uvgrid npix=(7400,7400) and res=104.
2017-07-26 13:25:38,475 - rtpipe - INFO - 	 Expect 0 thermal false positives per segment.
2017-07-26 13:25:38,476 - rtpipe - INFO - 
2017-07-26 13:25:38,476 - rtpipe - INFO - 	 Visibility memory usage is 0.5 GB/segment
2017-07-26 13:25:38,477 - rtpipe - INFO - 	 Imaging in 1 chunks using max of 35.9 GB/segment
2017-07-26 13:25:38,478 - rtpipe - INFO - 	 Grand total memory usage: 36.4 GB/segment
2017-07-26 13:25:38,653 - rtpipe.parsesdm - INFO - Reading scan 25, segment 7/14, times 10:29:09.225 to 10:29:09.665
2017-07-26 13:25:38,796 - rtpipe.parsesdm - INFO - Reading 88 ints starting at int 356
2017-07-26 13:25:39,385 - rtpipe.parsesdm - INFO - Found online flags for 387 antenna/time ranges.
2017-07-26 13:25:39,393 - rtpipe.parsesdm - INFO - Applied online flags to 0 ints.
2017-07-26 13:25:39,814 - rtpipe.parsesdm - WARNING - No BDF found for scans [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33]
2017-07-26 13:25:39,967 - rtpipe.parsesdm - INFO - Calculating uvw at 2016/09/17/10:29:09.445 for scan 25 of source FRB121102-off
2017-07-26 13:25:40,123 - rtpipe.parsecal - INFO - Read telcalfile 16A-496_sb32698778_1_02h00m_000.57648.37452900463.GN
2017-07-26 13:25:40,125 - py.warnings - WARNING - /Users/caseyjlaw/anaconda/lib/python2.7/site-packages/rtpipe/parsecal.py:418: RuntimeWarning: divide by zero encountered in divide
  badsols = n.where( (n.median(self.amp)/self.amp > threshold) & (self.flagged == False))[0]

2017-07-26 13:25:40,153 - rtpipe.parsecal - INFO - Frequency selection cut down to 2592 solutions
2017-07-26 13:25:40,156 - rtpipe.parsecal - INFO - Selection down to 432 solutions separated from given time by 7 minutes
2017-07-26 13:25:40,158 - rtpipe.parsecal - INFO - MJD: [ 57648.44188426]
2017-07-26 13:25:40,160 - rtpipe.parsecal - INFO - Source: ['J0555+3948']
2017-07-26 13:25:40,162 - rtpipe.parsecal - INFO - Solutions for 8 spw: ([ 2553.  2681.  2809.  2937.  3065.  3193.  3321.  3449.])
2017-07-26 13:25:40,163 - rtpipe.parsecal - INFO - Applying gain solution for chans from 0-30
2017-07-26 13:25:40,257 - rtpipe.parsecal - INFO - Applying gain solution for chans from 31-62
2017-07-26 13:25:40,341 - rtpipe.parsecal - INFO - Applying gain solution for chans from 63-94
2017-07-26 13:25:40,425 - rtpipe.parsecal - INFO - Applying gain solution for chans from 95-126
2017-07-26 13:25:40,512 - rtpipe.parsecal - INFO - Applying gain solution for chans from 127-157
2017-07-26 13:25:40,591 - rtpipe.parsecal - INFO - Applying gain solution for chans from 158-189
2017-07-26 13:25:40,676 - rtpipe.parsecal - INFO - Applying gain solution for chans from 190-221
2017-07-26 13:25:40,760 - rtpipe.parsecal - INFO - Applying gain solution for chans from 222-253
2017-07-26 13:25:40,847 - rtpipe - INFO - Flagging with flaglist: [('badap', 3.0, 0.2)]
2017-07-26 13:25:40,935 - rtpipe - INFO - Bad basepol flagging for chans 0-29 at 3.0 sigma: ants/pols [8]/[1], 0.44 % of total flagged
2017-07-26 13:25:40,986 - rtpipe - INFO - Bad basepol flagging for chans 0-29 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,082 - rtpipe - INFO - Bad basepol flagging for chans 30-61 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,144 - rtpipe - INFO - Bad basepol flagging for chans 30-61 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,236 - rtpipe - INFO - Bad basepol flagging for chans 62-93 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,293 - rtpipe - INFO - Bad basepol flagging for chans 62-93 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,385 - rtpipe - INFO - Bad basepol flagging for chans 94-125 at 3.0 sigma: ants/pols [8]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,445 - rtpipe - INFO - Bad basepol flagging for chans 94-125 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,544 - rtpipe - INFO - Bad basepol flagging for chans 126-157 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,600 - rtpipe - INFO - Bad basepol flagging for chans 126-157 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,692 - rtpipe - INFO - Bad basepol flagging for chans 158-189 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,752 - rtpipe - INFO - Bad basepol flagging for chans 158-189 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:41,851 - rtpipe - INFO - Bad basepol flagging for chans 190-221 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:41,909 - rtpipe - INFO - Bad basepol flagging for chans 190-221 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:42,004 - rtpipe - INFO - Bad basepol flagging for chans 222-253 at 3.0 sigma: ants/pols [14]/[1], 0.47 % of total flagged
2017-07-26 13:25:42,063 - rtpipe - INFO - Bad basepol flagging for chans 222-253 at 3.0 sigma: ants/pols []/[], 0.00 % of total flagged
2017-07-26 13:25:42,065 - rtpipe - INFO - Subtracting mean visibility in time...
2017-07-26 13:25:42,201 - rtpipe - INFO - Returning prepared data...
2017-07-26 13:25:42,352 - rtpipe.parsesdm - WARNING - No BDF found for scans [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33]
2017-07-26 13:25:42,496 - rtpipe.parsesdm - INFO - Calculating uvw at 2016/09/17/10:29:09.445 for scan 25 of source FRB121102-off

Select VLA data for comparison to Arecibo


In [22]:
st, data, u, v, w = read[key]
st['dmarr'] = [dm0]
scale = getscale(st)
dataph = correctdata(st, data*scale, u, v, w, corr='ph,dm')
intrange = (30, 60)
dint = intrange[1] - intrange[0]
specVLA = dataph[intrange[0]:intrange[1]].mean(axis=3).mean(axis=1).real
specVLA = (specVLA-specVLA.mean())/specVLA.std()  # in SNR units
print(specAO.shape, specVLA.shape)


((256, 512), (30, 254))

Calculate dm delay and topocentric correction for VLA


In [23]:
dmt_vla = dmdelay(dm0, st['freq'][-1], 1e4)/(24*3600)
topo_vla = 75.7533100 # seconds, from presto
mjdfrac_vla = st['segmenttimes'][7][0] + intrange[0]*st['inttime']/(24*3600) - mjd0 - dmt_vla + topo_vla/(24*3600)
print('VLA dMJD: {0}'.format(mjdfrac_vla))
print('AO dMJD: {0}'.format(mjdfrac_ao))
print('Diff: {0} s'.format((mjdfrac_ao-mjdfrac_vla)*24*3600))


VLA dMJD: 0.437788650962
AO dMJD: 0.437788874842
Diff: 0.0193432491544 s

In [24]:
tbary_vla = np.linspace(mjdfrac_vla, mjdfrac_vla+(dint*st['inttime']/(24*3600)), dint)
freqs_vla = 1e3*st['freq']
extent_vla = (tbary_vla[0], tbary_vla[-1], freqs_vla[0], freqs_vla[-1])

In [25]:
# put them together
fig = pl.figure(figsize=(8,8))
fig.add_subplot(211)
pl.imshow(specVLA.transpose(), interpolation='nearest', origin='bottom', extent=extent_vla, aspect='auto')
fig.add_subplot(212)
pl.imshow(specAO, interpolation='nearest', aspect='auto', extent=extent_ao)


Out[25]:
<matplotlib.image.AxesImage at 0x124a379d0>

In [26]:
tpk_vla = tbary_vla[np.where(specVLA.mean(axis=1) == specVLA.mean(axis=1).max())][0]
print('%.10f' % tpk_vla)


0.4377894891

Regrid VLA and AO data to same fixed image grid


In [27]:
gap = 30
specregrid = np.zeros(shape=(nch1+st['nchan']+gap, len(tbary_ao)))

for idx in range(len(tbary_ao)):
    specregrid[254+gap:, idx] = specAO[:, idx]

idxs_vla = [np.argmin(np.abs(tbary_vla-tbary_ao[i])) for i in range(len(tbary_ao))]
for idx_ao in range(len(idxs_vla)):
    idx_vla = idxs_vla[idx_ao]
    specregrid[:254, idx_ao] += specVLA[idx_vla, ::-1]

Plot it


In [28]:
# assumes fixed relative gain between VLA and AO == 3
fig = pl.figure(figsize=(12,12))
ax = fig.add_subplot(211)
pl.imshow(specregrid, interpolation='nearest', aspect='equal', vmax=0.8*specregrid.max(), cmap='Greys')#, vmin=-0.8)
ax.fill_between(np.arange(0, len(tbary_ao)), 254*np.ones(len(tbary_ao)), 
                (254+gap)*np.ones(len(tbary_ao)), facecolor='k')

pl.xlabel('Time (ms)', fontsize=14)
pl.ylabel('Frequency (GHz)', fontsize=14)
xticks = np.arange(0, 600, 100)
pl.xticks(xticks, np.array(tbin*2e-2*xticks, dtype='int'))
ntot = nch1+len(st['freq'])+gap
yticks = np.arange(0, ntot, 80)
print(yticks)
pl.yticks(yticks)
pl.yticks(yticks, [st['freq'][-1], st['freq'][-80], st['freq'][-160], st['freq'][-240],
                   np.round(1e-3*freqs_ao[ntot-320], 3), 
                   np.round(1e-3*freqs_ao[ntot-400], 3), 
                   np.round(1e-3*freqs_ao[ntot-480], 3)])
pl.xlim(150, len(tbary_ao))
pl.ylim(ntot-50, 0)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='k')
ax.yaxis.set_tick_params(width=3, color='k')
ax.text(170, 330, "Arecibo", rotation=90, fontsize=20)
ax.text(170, 130, "VLA", rotation=90, fontsize=20)
fig.savefig('aovla_spec.pdf', format='pdf')


[  0  80 160 240 320 400 480]

Calculate residual sweep and compare to DM model


In [29]:
# some residual dm sweep?
print(24*3600*(tbary_ao[np.where(specAO[128:].mean(axis=0) == specAO[128:].mean(axis=0).max())][0] -
               tbary_ao[np.where(specAO[:128].mean(axis=0) == specAO[:128].mean(axis=0).max())][0]))

dmt_ao_delta = dmdelay(dm0, obsfreq*1e-3, st['freq'][-1]) - dmdelay(560, obsfreq*1e-3, st['freq'][-1])
dmt_ao_inter = dmdelay(dm0, obsfreq*1e-3, (obsfreq-500)*1e-3) - dmdelay(560, obsfreq*1e-3, (obsfreq-500)*1e-3)
print(dmt_ao_delta, dmt_ao_inter)


0.00147730193518
(-0.0062022839361107751, 0.007798684063229699)

Burst spectra with and without detections

Compare Effelsberg, VLA, Arecibo measured burst fluxes and limits


In [30]:
# 4 bursts with ao, vla, effelsberg coverage
# (ami-la covers all 4 vla bursts)
# 57643 (AO-C, *VLA-S)
# 57645 (AO-L, *VLA-S)
# 57648 (*AO-L, *VLA-S, Eff-C)
# 57649 (AO-L, *VLA-S, Eff-C)
# * shows detections
# limits assume 2 ms pulse width

#s43d = np.array([[3.0, 4.0], [0.508, 0.0036*5]])
#s43l = np.array([[4.0], [0.0036*5]])  # useless data!
s45d = np.array([[1.38, 3.0], [0.002*5, (5/2.)*0.064]])
s45l = np.array([[1.38], [0.002*5]])  # too high?
s48d = np.array([[1.4, 3.0, 4.85], [0.057, (5/2.)*0.111, 0.028*5]])  # fixed for 2 ms width
s48l = np.array([[4.85], [0.028*5]])
s49d = np.array([[1.42, 3.0, 4.9], [0.002*5, (5/2.)*0.167, 0.028*5]])
s49l = np.array([[1.42, 4.9], [0.002*5, 0.028*5]])

In [31]:
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(111)
# overplot upper limit symbols
#ax.plot(s43l[0], s43l[1], 'cv', ms=10)
ax.plot(s45l[0], s45l[1], 'kv', ms=10)
ax.plot(s48l[0], s48l[1], 'kv', ms=10)
ax.plot(s49l[0], s49l[1], 'kv', ms=10)
# plot lines with points first
#ax.plot(s43d[0], s43d[1], 'c.-', ms=10, label='12 Sep 2016')
ax.plot(s45d[0], s45d[1], 'k.--', ms=10, label='57645')
ax.plot(s48d[0], s48d[1], 'k.-.', ms=10, label='57648')
ax.plot(s49d[0], s49d[1], 'k.-', ms=10, label='57649')
ax.set_xlabel('Frequency (GHz)', fontsize=14)
ax.set_ylabel('Integrated Flux density (Jy; 2 ms)', fontsize=14)
ax.legend(fontsize=14)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='k')
ax.yaxis.set_tick_params(width=3, color='k')
fig.savefig('multispec.pdf', format='pdf')


Calculate implied spectral index and limits


In [32]:
def speclim(s):
    freqs = s[0]
    fluxes = s[1]
    for i in range(len(freqs)-1):
        freqi = freqs[i]
        freqi1 = freqs[i+1]

        # correct for plotting offset
        if freqi <= 1.5:
            freqi = 1.4
        elif freqi1 > 4.2:
            freqi1 = 4.85

        print(freqi, freqi1)
        print(fluxes[i], fluxes[i+1])
        print(np.log10(fluxes[i]/fluxes[i+1])/np.log10(freqi/freqi1))

In [33]:
speclim(s45d)


(1.4, 3.0)
(0.01, 0.16)
3.63789924804

In [34]:
speclim(s48d)


(1.4, 3.0)
(0.057000000000000002, 0.27750000000000002)
2.07674384925
(3.0, 4.85)
(0.27750000000000002, 0.14000000000000001)
-1.42428464478

In [35]:
speclim(s49d)


(1.4, 3.0)
(0.01, 0.41750000000000004)
4.89634344931
(3.0, 4.85)
(0.41750000000000004, 0.14000000000000001)
-2.27460139728