Initial commit

from pyclass import comm,gdict

from math import *

import numpy as np

from astropy.time import Time,TimeDelta

from scipy.interpolate import interp1d

import sys

import os

gain_curve_poly=[0.7929185,0.005900533,-4.203179E-5,0]

polarization_index={'LCP':1,'RCP':2}   

def gain_curve(el,c):

    return c[0]+c[1]*el+c[2]*el**2

def readOpacityFile(taufile):

    my_data = np.genfromtxt(taufile,dtype=float)

    print my_data

    return my_data

def findMinDist(time,time_tau,tau):

    indice=0

    minimo=fabs(time.mjd-time_tau[0])

    print minimo

    for idx,val in enumerate(time_tau):

         deltatime=fabs((time.mjd-val))

         if (deltatime <minimo):

             minimo=deltatime

             indice=idx

    tau_min=tau[indice]

    time_tau_min=time_tau[indice]

    return time_tau_min,tau_min

def readGains(gainfile):

  '''

  Read gains.txt

  '''

  my_data = np.genfromtxt(gainfile,dtype=float)

  return my_data

def findPolarizationIndex(polarization='RCP'):

    return polarization_index[polarization]

def gildasFind(polarization='RCP'):

    gildas_findpolarization_string='find /line *'+polarization

    return gildas_findpolarization_string

def calibrate(file_in,taufile,gainfile):

   #tau=0.09

   gildas_findpolarization_string='find'

   print "gildas:>>>>>>>>>>>",gildas_findpolarization_string 

   tau_list=readOpacityFile(taufile)

   gain_jy=readGains(gainfile)

   time_int= tau_list[:,0]

   time_gain =gain_jy[:,0]

   filename=file_in.split('.')[0]

   #################################

   outfile=filename+'_tstar.med'

   print filename

   comm('file in '+file_in)

   comm('file ou '+outfile + ' single /overwrite')

   comm(gildas_findpolarization_string)

   for i in range(0,gdict.found): 

            print i

            comm('get n')

            elevation=gdict.elevation  

            ut=degrees(gdict.r.head.gen.ut)/360. # ut is in radians. Converted to degrees and then in fraction of the day  

            spectrum_polarization=str(gdict.line)

            if 'RCP' in spectrum_polarization:

                tau_int=  tau_list[:,polarization_index['RCP']]

            else:

                tau_int=  tau_list[:,polarization_index['LCP']]

            time_scan=Time(gdict.r.head.gen.dobs+2460549.5+ut,format='jd') # .5 of difference because jd starts from midday. 

            time_tau,tau_min=findMinDist(time_scan,time_int,tau_int)            

            coeff=exp(tau_min/sin(elevation))

            print tau_min

            comm('mul '+str(coeff))

            print 'Ta*  mul '+str(coeff),'tau',tau_min,degrees(elevation)

            comm('write')

   # 

   file_in=outfile # the tstar output file is the

   filename=file_in.split('.')[0]

   outfile=filename+'_gain.med'

   print "gildas:>>>>>>>>>>>",gildas_findpolarization_string 

   comm('file in ' + file_in)

   comm('file ou '+outfile + ' single /overwrite')

   comm(gildas_findpolarization_string)

   for i in range(1,gdict.found+1): 

            print i

            comm('get '+str(i))

            elevation=degrees(gdict.elevation)

            coeff=1./ gain_curve(elevation,gain_curve_poly)  

            comm('mul '+str(coeff))

            print 'Gain Curve: mul '+str(coeff),1/coeff,elevation

            comm('write')

   file_in=outfile # the tstar output file is the

   filename=file_in.split('.')[0]

   outfile=filename+'_jy.med'

   print file_in, outfile

   comm('file in ' + file_in)

   comm('file ou '+outfile + ' single /overwrite')

   comm(gildas_findpolarization_string)

   for i in range(1,gdict.found+1): 

            print i

            comm('get '+str(i))

            elevation=gdict.elevation  

            ut=degrees(gdict.r.head.gen.ut)/360. # ut is in radians. Converted to degrees and then in fraction of the day  

            spectrum_polarization=str(gdict.line)

            if 'RCP' in spectrum_polarization:

                gain_int = gain_jy[:,polarization_index['RCP']]

            else:

                gain_int=  gain_jy[:,polarization_index['LCP']]

            time_scan=Time(gdict.r.head.gen.dobs+2460549.5+ut,format='jd') # .5 of difference because jd starts from midday. 

            time_gain_min,gain_min=findMinDist(time_scan,time_gain,gain_int)

            coeff=1./ gain_min 

            comm('mul '+str(coeff))

            print 'Gain: mul '+str(coeff),gain_min

            comm('write')

def printUsage(argv):

    print "Usage:"

    print "python "+argv[0]+" classfile taufile gainfile"

    print "The classfile is the file generated by discos2class"

def main(argv):

   if len(argv)==4:

      calibrate(argv[1],argv[2],argv[3])

   else:

      printUsage(argv)

if __name__ == "__main__":

     main(sys.argv)

 No newline at end of file

58008.73098  0.093  0.093

58008.79219  0.097  0.097

58009.68114  0.094  0.098

58009.68824  0.092  0.091

58013.78494  0.095  0.098

 No newline at end of file

58008.716983   0.125   0.128

58008.733996   0.125   0.128

58008.778534   0.118   0.120

58008.794852   0.116   0.119

58009.695779   0.154   0.156

58013.796575   0.129   0.140