Effelsberg 100m Teleskop

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
===== Read raw data (FFTS format) =====
<code python>
import numpy as np
import os

def getRawData(filename,metainfo=None):
    """read raw data from FFTS into numpy record arrays
    metainfo ... if provided (e.g. if known from a previous file 
    of the same type, this prevents parsing it again
    
    returns alldata,metainfo
    alldata contains the following fields (use alldata['FIELD'] to access)
    IEEE, DATFMT, DATALENGTH, BE_IDENT, TIMESTAMP, INTTIME, PHASENUM, BENUMSEC,
    BLOCKINGFACTOR, DATABLOCK
    DATABLOCK is itself a record, containing:
    SECTION, NUMCHANNELS, DATA
    
    to access the raw data use 'alldata['DATABLOCK']['DATA']' which is a 3D array,
    of shape (dump,baseband,specchannel), e.g., to access first dump, first baseband 
    spectrum 'alldata['DATABLOCK']['DATA'][0,0]', last dump, second baseband: 
    'alldata['DATABLOCK']['DATA'][-1,1]', etc.
    
    use: 
    
        import dateutil.parser
        alltimestamps=[dateutil.parser.parse(a) for a in alldata['TIMESTAMP']]
        allmjd=np.array([cal_mjd(dt.month,
                                 dt.day+(dt.hour+dt.minute/60.+dt.second/3600.+dt.microsecond/3.6e9)/24.,
                                 dt.year) for dt in alltimestamps])
    
    to compute MJD values from the timestamps
    
    """
    totalsize=os.path.getsize(fname)
    if metainfo is None:
        # getting endian
        ieee=np.fromfile(fname , dtype = np.dtype([('IEEE','|S4')]) , count = 1 )
        endian='>' if ieee=="IEEE" else '<'
        metainfo=np.fromfile(fname , dtype = np.dtype([ ('IEEE','|S4'), ('DATFMT','|S4'),
        ('DATALENGTH','%si4'%endian), ('BE_IDENT','|S8'), ('TIMESTAMP','|S28'),
        ('INTTIME','%si4'%endian), ('PHASENUM','%si4'%endian), ('BENUMSEC','%si4'%endian),
        ('BLOCKINGFACTOR','%si4'%endian), ('SECTION','%si4'%endian),
        ('NUMCHANNELS','%si4'%endian) ]) , count = 1 )
    else:
        endian='>' if metainfo['IEEE'][0]=="IEEE" else '<'
    #
    intorfloat=np.int32 if metainfo['DATFMT'][0]=="I" else np.float32
    datfmt=('%si4'%endian) if metainfo['DATFMT'][0]=="I" else ('%sf4'%endian)
    datalength=metainfo['DATALENGTH'][0]
    numchannels=metainfo['NUMCHANNELS'][0]
    
    datablock=np.dtype([('SECTION','%si4'%endian), ('NUMCHANNELS','%si4'%endian),
    ('DATA',datfmt,numchannels)])
    record=np.dtype([ ('IEEE','|S4'), ('DATFMT','|S4'), ('DATALENGTH','%si4'%endian), 
    ('BE_IDENT','|S8'), ('TIMESTAMP','|S28'), ('INTTIME','%si4'%endian), 
    ('PHASENUM','%si4'%endian), ('BENUMSEC','%si4'%endian), ('BLOCKINGFACTOR','%si4'%endian),
    ('DATABLOCK',datablock,metainfo['BLOCKINGFACTOR'][0]*metainfo['BENUMSEC'][0])  ])
    record_length=record.itemsize
    num_records=totalsize/record_length
    alldata=np.fromfile(fname , dtype =record, count = num_records )
    return alldata,metainfo
#    

fname="/be4-daten/AFFTS/2788.0001"
alldata,metainfo=getRawData(fname)

print alldata['DATABLOCK']['SECTION']
print alldata['DATABLOCK']['NUMCHANNELS']
print alldata['DATABLOCK']['DATA'].shape
</code>

If you want to convert the timestamps to MJD, you can do the following:

<code python>
import dateutil.parser

def cal_mjd(mn,dy,yr):
    """calculate MJD from a data
    mn ... month (int)
    yr ... year (int)
    dy ... day (float, i.e. including fractions of a day)
    
    returns MJD
    
    (code snippet adapted from pyephem, which is using libastro)"""
    m=int(mn)
    y=int(yr)
    if y<0: y+=1
    if mn<3: mn+=12
    if mn<3: y-=1
    b=0
    if not ( (yr < 1582) | ((yr == 1582) & ((mn < 10) | ((mn == 10) & (dy < 15))))):
        a=y//100
        b=2 - a + a//4
    
    c=0
    if y<0:
        c=int((365.25*y) - 0.75) - 694025L
    else:
        c=int(365.25*y) - 694025L
    d=int(30.6001*(m+1))
    return b + c + d + dy - 0.5  +15020L -0.5
#

alltimestamps=[dateutil.parser.parse(a) for a in alldata['TIMESTAMP']]
allmjd=np.array([cal_mjd(dt.month,
                         dt.day+(dt.hour+dt.minute/60.+dt.second/3600.+dt.microsecond/3.6e9)/24.,
                         dt.year) for dt in alltimestamps])
</code>
Effelsberg 100m Teleskop

User Tools

Site Tools

Page Tools
