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==""):