casa6

#!/usr/bin/env python

#

# linfeedpolhelpers.py

#

# History:

#  v1.0 (gmoellen, 2012Oct24) == initial version

#

# This script defines several methods currently needed for

# instrumental polarization calibration for the linear

# feed basis.

# To import these functions, type (at the CASA prompt):

#

# from recipes.linfeedpolhelpers import *

#

#

import os

​

def qufromgain(gt,badspw=[]):

​

    me.doframe(me.observatory('alma'))

​

    tb.open(gt+'/ANTENNA')

    pos=tb.getcol('POSITION')

    meanpos=pl.mean(pos,1)

    frame=tb.getcolkeyword('POSITION','MEASINFO')['Ref']

    units=tb.getcolkeyword('POSITION','QuantumUnits')

    mpos=me.position(frame,

                     str(meanpos[0])+units[0],

                     str(meanpos[1])+units[1],

                     str(meanpos[2])+units[2])

    me.doframe(mpos)

​

    # _geodetic_ latitude

    latr=me.measure(mpos,'WGS84')['m1']['value']

​

    tb.open(gt+'/FIELD')

    nfld=tb.nrows()

    dirs=tb.getcol('DELAY_DIR')[:,0,:]

    tb.close()

    print 'Found as many as '+str(nfld)+' fields.'

​

    tb.open(gt+'/SPECTRAL_WINDOW')

    nspw=tb.nrows()

    tb.close()

    print 'Found as many as '+str(nspw)+' spws.'

​

    R=pl.zeros((nspw,nfld))

    Q=pl.zeros((nspw,nfld))

    U=pl.zeros((nspw,nfld))

    mask=pl.ones((nspw,nfld),dtype=bool)

​

    if (len(badspw)>0):

        mask[badspw,:]=False

​

    tb.open(gt)

    for ifld in range(nfld):

        for ispw in range(nspw):

            if not mask[ispw,ifld]:

                continue

            st=tb.query('FIELD_ID=='+str(ifld)+' && SPECTRAL_WINDOW_ID=='+str(ispw))

            nrows=st.nrows()

            if nrows > 0:

                rah=dirs[0,ifld]*12.0/pi

                decr=dirs[1,ifld]

                times=st.getcol('TIME')

                gains=st.getcol('CPARAM')

                ants=st.getcol('ANTENNA1')

​

                ntimes=len(pl.unique(times))

                nants=ants.max()+1

​

                print nrows, ntimes, nants

​

                # times

                time0=86400.0*floor(times[0]/86400.0)

                rtimes=times-time0

​

                # amplitude ratio

                amps=pl.absolute(gains)

                amps[amps==0.0]=1.0

                ratio=amps[0,0,:]/amps[1,0,:]

                # parang

                parang=pl.zeros(len(times))

                for itim in range(len(times)):

                    tm=me.epoch('UTC',str(times[itim])+'s')

                    last=me.measure(tm,'LAST')['m0']['value']

                    last-=floor(last)  # days

                    last*=24.0  # hours