casa6

# opacal()

#

# Recipe to create a frequency dependent opacity calibration table.

# The gencal generated 'opac' table is constant in time and frequency

# (single value per spectral window), this function uses the WEATHER table

# to create an opacity table that can be interpolated in time and frequency.

# This recipe matches the Miriad calculation for the ATCA when the

# height parameter is left at the default and the asap parameter is set

# to True. It may give a reasonable approximation for other observatories

# if the height parameters is adjusted, but some of the atmospheric

# parameters in the code may need changing too.

#

# To access this function, type (at the CASA prompt):

#

# from recipes.opacity import opacal

# or 

# execfile(casadef.python_library_directory+'/recipes/opacity.py')

#

# MHW, Feb 2016

from taskinit import tbtool,cu,casac

from task_gencal import gencal as mygencal

import pylab as pl

​

def opacal(vis,calname,asap=True,interpolate=0,height=200):

  """Create an opacity calibration table for the MeasurementSet given.

  If calname is empty, one is created based on vis and returned.

  Use asap=True to use the ASAP/Miriad atmosphere model and

  asap=False to use the CASA ATM model.

  Use interpolate>1 to use piecewise linear interpolation across each spectrum

  or leave at 0 to use a constant value per spectral window.

  The height parameter should specify the observatory height above mean

  sea level (not height above WGS84 ellipsoid)"""

​

  mytb = tbtool()

  # get data we need from weather table

  try:

    mytb.open(vis+'/WEATHER')

    time = mytb.getcol('TIME')

    press = mytb.getcol('PRESSURE') #Pa

    temp = mytb.getcol('TEMPERATURE') # C (sometimes returns Kelvin)

    hum = mytb.getcol('REL_HUMIDITY') # %

  except Exception, e:

    print "Error: Cannot read required data from WEATHER table"

    return

  mytb.close()

  if (time.size==0):

    print "Error: no data in WEATHER table, cannot calculate opacities"

    return

  # check array name

  mytb.open(vis+'/OBSERVATION')

  tel=mytb.getcell('TELESCOPE_NAME',0)

  mytb.close()

  # ATCA weather station records pressure adjusted to sea level, ATM code wants

  # actual surface pressure - we undo the (fixed) correction here. ASAP code

  # does the correction itself when height>0.

  pfac=1.0

  if (tel=='ATCA'):

    pfac=0.975

  # create opac table and read the frequencies, extend opac table to fit an

  # entry for each timestamp

  if (calname==""):

    calname=vis.split(".")[0]+".opac.caltable"

  cu.removetable(calname)

  mygencal(vis,calname,'opac',infile='',spw='',antenna='',pol='',parameter=[0.0])

  mytb.open(calname+'/SPECTRAL_WINDOW')

  fref=mytb.getcol('REF_FREQUENCY')/1e9

  bw=mytb.getcol('TOTAL_BANDWIDTH')/1e9

  nchan=mytb.getcol('NUM_CHAN')

  mytb.close()

  nf=1

  if (interpolate): 

    nf=max(2,interpolate)

  f=pl.zeros(nf*fref.size)

  print "Using %i frequencies to calculate opacity for each spectral window" % nf

  if (interpolate):

    for i in range(nf):

      f[i::nf]=fref-bw/2+i*bw/(nf-1)

  else:

    f=fref

  mytb.open(calname,nomodify=False)

  spw=mytb.getcol('SPECTRAL_WINDOW_ID')

  ant1=mytb.getcol('ANTENNA1')

  nr=mytb.nrows()

  addnr=nr*(time.size-1)

  mytb.addrows(addnr)

  nrow2=mytb.nrows()

​

  # work out opacity for each time stamp in weather table and

  # write values to opacity table