casa6

"""

Fit a constant tau component to the 4 ALMA WVR channels, 

then remove this from the 4 temperature measurements tsrc(0-3)

​

Based on rem_cloud.py by B. Dent version 12 Aug 2015, see ALMA CSV-3189

"""

​

from taskinit import *

​

import numpy as np

import pylab as pl

from matplotlib import pyplot

import scipy.optimize as opt

import os

import math

import time

​

def rmc_func1(x,a,b,c):    

    return a*(x-b)**2+c

​

def rmc_weighted_avg_and_std(values, weights):

    """

    Return the weighted average and standard deviation.

​

    values, weights -- Numpy ndarrays with the same shape.

    """

    average = np.average(values, weights=weights)

    variance = np.average((values-average)**2, weights=weights)

    return (average, math.sqrt(variance))

​

def rmc_approxCalc(tsrc0, tsrc1, tsrc2, tsrc3,

                   m_el, Tamb, verb=False):

    m_el=m_el/57.295   # convert to radians

    mean_pwt=0.0; raw_mean_pwt=0.0

​

    # correct for coupling to the sky, assuming T_loss=275K, eta_c=0.98

    eta_c=0.98

    T_loss=275.0

​

    tsrc0=(tsrc0-(1.0-eta_c)*T_loss)/eta_c

    tsrc1=(tsrc1-(1.0-eta_c)*T_loss)/eta_c

    tsrc2=(tsrc2-(1.0-eta_c)*T_loss)/eta_c

    tsrc3=(tsrc3-(1.0-eta_c)*T_loss)/eta_c

​

    #    m_az, m_el, m_ts= mc.actualAzEl()

​

    pw=[0.0,0.0,0.0,0.0];     pw_noc=[0.0,0.0,0.0,0.0]

    site="AOS"

    # set up fitting constants depending on the site:

    if site == "OSF":   # constants:

       tau0=[0.024,0.02,0.009,0.01]

       r=[1.205,0.402,0.177,0.116]

       # approximate physical temp of atmosphere:

       Tphys=267.0

​

    elif site == "AOS":  # constants:

       tau0=[0.027,0.02,0.010,0.01]

       r=[1.193,0.399,0.176,0.115]

       # approximate physical temp of atmosphere ,based on ambient temperature, Tamb, in centrigrade:

       Tphys=270.0 +Tamb

​

       if Tphys<tsrc0:

           Tphys=tsrc0+1.0

           if verb:

               casalog.post('  fixed physical temperature to be tsrc0', 'INFO')

​

    tsrcn=np.zeros(4)

    teln=np.zeros(4)

​

    tel=[0.0,0.0,0.0,0.0]

    tz=[0.0,0.0,0.0,0.0]

    wt=[0.0,0.0,0.0,0.0]

​

    # calculate transmissions:

    tel[3]=(1.0-tsrc3/Tphys)

    tel[2]=(1.0-tsrc2/Tphys)

    tel[1]=(1.0-tsrc1/Tphys)

    tel[0]=(1.0-tsrc0/Tphys)

​

    for i in range(4):

      if tel[i]<0.02:

         tel[i]=0.02

      wt[i]=1.0-(abs(tel[i]-0.5)/0.49)**2.0    # weights

​

    if verb:

        casalog.post( '  weights 0-3: '+str(wt), 'INFO')

        casalog.post( '  tsrc 0-3:    '+str(tsrc0)+', '+str(tsrc1)+', '+str(tsrc2)+', '+str(tsrc3), 'INFO')

​

    use=1