# # Author: Steven T. Myers (NRAO) 2001-12-14 # Set of functions to apply opacity and gain correction to vla data # # Contents: # vlagaincor(par) apply gain correction to DATA column # opac(tau,el) calculate opacity correction at vector el # gaincv(a1,a2,el) gain curve for vecs an1,an2,el # vlagainfill(band) fill gain coefficient array # # Usage: # # 1. define parameters in record par # 2. if needed, make ms from fits using makemsf(par) from automap.g # 3. run vlagaincor(par) # include 'measures.g'; include 'table.g'; include 'quanta.g'; # # Example set up of paramters # par.msf := '20010204AK.ms'; par.tau0 := 0.10; par.band := 'k'; # # Function to apply gain and opacity correction to ms # vlagaincor := function( par ) { # msf := par.msf; tau0 := par.tau0; band := par.band; # # Fill gain coefficients # c := vlagainfill( band ); # # Set the antenna number lookup # anam := spaste(msf, '/ANTENNA'); atab := table(anam); nant := atab.nrows(); anam := atab.getcol("NAME"); printf(' Found %d antennas in table\n',nant); atab.close(); # # Set the frame to be the appropriate telescope # onam := spaste(msf, '/OBSERVATION'); otab := table(onam); tel := otab.getcell("TELESCOPE_NAME",1); printf(' Data is from telescope %s\n',tel); telpos:=dm.observatory(tel); dm.measure(telpos,'WGS84'); # convert to geodetic (why????) dm.doframe(telpos); otab.close(); # # Get source directions # fnam := spaste(msf, '/FIELD'); ftab := table(fnam); nfld := ftab.nrows(); printf(' Found %d fields\n',nfld); fref := ftab.getcolkeyword("REFERENCE_DIR","MEASINFO"); printf(' Reference frame %s\n',fref.Ref); fdirs := ftab.getcol("REFERENCE_DIR"); fields := ftab.getcol("NAME"); ftab.close(); # # Process directions - store in record of directions # dirq := [=]; for ( i in 1:nfld ) { raq := dq.quantity(fdirs[1,,i],'rad'); deq := dq.quantity(fdirs[2,,i],'rad'); dirq[i] := dm.direction(fref.Ref,raq,deq); } # # Make ms table # mst := table(msf, readonly=F); # # Make ms table iterator and iterate # ti := tableiterator(mst, "TIME"); while ( ti.next() ) { st := ti.table(); t := st.getcell("TIME",1); tq := dq.quantity(t,'s'); tymd := dq.time(tq,form='ymd'); print ' Processing timestamp ',tymd; tim := dm.epoch('tai',tymd) ; dm.doframe(tim); # proxy source as first row source for now, get EL for this set ifld := st.getcell("FIELD_ID",1)+1; # (+1 corrects c++ indexing) azel := dm.measure(dirq[ifld],'AZEL'); el := azel.m1; # probably should use getvalue # # Get visibilties - this will be vector [shape=[4 1 nvis] ] # nvis := st.nrows(); vis := st.getcol("DATA"); # # Get AN number vectors (should do something about invalid antennas) # an1v := st.getcol("ANTENNA1")+1; an2v := st.getcol("ANTENNA2")+1; an1 := as_integer( anam[ an1v ]); an2 := as_integer( anam[ an2v ]); # # Construct el numeric vector # elv := array( el.value, nvis ); # # construct opacity corrections = ( 1 - tau/sin(EL) )^(-1) # oc := opac( tau0, elv ); # # construct gain corrections # gc := gaincv( an1, an2, elv, c ); vc := oc/gc; # # make compatible with vis vector # vcv := array(1,4,nvis); for ( i in 1:nvis ) { vcv[1:4,i] := array( vc[i], 4 ); } # # apply to data vector # vis *:= vcv; # # put back in ms # st.putcol("DATA",vis); printf(' Corrected %s with %d rows at EL = %.2f\n', fields[ifld], nvis, dq.convert(el,'deg').value ); } # # Done # ti.done(); mst.close(); print ' Completed gain and opacity correction\n'; } # # Define the opacity correction # opac := function( tau0, el ) { # # Do atmospheric absorption correction on vector of el # tau(el) = tau0/sin(el) -> oc = 1/[ 1 - tau(el) ] # maximum tau = 1e3 sel := sin( el ); tau := tau0/sel; if ( tau > 1.e4 ) { tau := 1.e3 }; oc := 1.0/( 1.0 - tau ); return oc; } # # Define a gain curve (would index these by antenna eventually) # gaincv := function( an1, an2, el, c) { # # Input vectors (same length) an1,an2,el # Assumes gain curve information is stored in array c[3,nant] # Elevation el in radians # Antenna number an1,an2 elements must correspond to columns in c # Example: Kband AN2 IF1 L (Oct00) # c[1:3,2] := [1.0041e+00, -6.0815e-04, 2.2647e-05 ]; # Need to pull out appropriate coefficients array pick g1 := c[1,an1] + c[2,an1]*el + c[3,an1]*el*el; g2 := c[1,an2] + c[2,an2]*el + c[3,an2]*el*el; g := sqrt( g1*g2 ); # Need to sort out abnormally low gains g[ g < 1.e-4 ] := 1.e-4; return g; } # # Define function to fill array of gain coefficients # vlagainfill := function( band ) { # fill array as vector if ( band == 'K' | band == 'k' ) { c := [ 1.0645E+00, -2.1800E-03, 1.9091E-06, # AN 1 9.9400E-01, 7.1267E-04, -2.1991E-05, # AN 2 9.0775E-01, 3.4519E-03, -3.2439E-05, # AN 3 8.7106E-01, 4.4819E-03, -3.9123E-05, # AN 4 9.5978E-01, 2.0537E-03, -2.6473E-05, # AN 5 9.6549E-01, 1.7126E-03, -2.1309E-05, # AN 6 9.0808E-01, 3.2659E-03, -2.9129E-05, # AN 7 9.2997E-01, 2.7050E-03, -2.6238E-05, # AN 8 9.4561E-01, 2.4109E-03, -2.7114E-05, # AN 9 9.3171E-01, 2.4068E-03, -2.1307E-05, # AN 10 9.9932E-01, 3.5149E-04, -2.0688E-05, # AN 11 8.4878E-01, 4.8514E-03, -3.9183E-05, # AN 12 9.1357E-01, 2.8711E-03, -2.4172E-05, # AN 13 8.9678E-01, 4.0089E-03, -3.9249E-05, # AN 14 9.8260E-01, 1.4266E-03, -3.0393E-05, # AN 15 8.7698E-01, 4.1427E-03, -3.5156E-05, # AN 16 8.5035E-01, 6.7215E-03, -7.5902E-05, # AN 17 9.9495E-01, 6.8603E-04, -2.4564E-05, # AN 18 9.7246E-01, 1.7074E-03, -2.6990E-05, # AN 19 9.3173E-01, 3.3235E-03, -4.0633E-05, # AN 20 1.0011E+00, 2.3808E-04, -2.2524E-05, # AN 21 9.1979E-01, 3.0976E-03, -3.0121E-05, # AN 22 9.8703E-01, 9.1300E-04, -1.6280E-05, # AN 23 9.9937E-01, 2.3530E-04, -1.4748E-05, # AN 24 8.5658E-01, 4.8396E-03, -4.1113E-05, # AN 25 9.7852E-01, 8.5349E-04, -8.7236E-06, # AN 26 9.2672E-01, 3.0128E-03, -3.1142E-05, # AN 27 9.4566E-01, 2.9426E-03, -4.0201E-05, # AN 28 1.0000E+00, 1.3900E-03, 0.0000E+00, # AN 29 9.4020E-01, 2.4481E-03, -2.9080E-05 ] # AN 30 = avg print ' Filled gain coefficients for VLA at Kband'; } else { if ( band == 'Q' | band == 'q' ) { c := [ 8.5390E-01, 6.9920E-03, -9.7635E-05, # AN 1 na = avg 8.5390E-01, 6.9920E-03, -9.7635E-05, # AN 2 na = avg 8.2699E-01, 8.6734E-03, -1.1053E-04, # AN 3 8.5282E-01, 6.9641E-03, -8.3056E-05, # AN 4 7.8219E-01, 6.6175E-03, -4.7822E-05, # AN 5 8.6940E-01, 7.1348E-03, -1.0007E-04, # AN 6 1.0000E+00, 0.0000E+00, 0.0000E+00, # AN 7 no rx 7.9055E-01, 1.0010E-02, -1.2145E-04, # AN 8 1.0000E+00, 0.0000E+00, 0.0000E+00, # AN 9 no rx 7.9209E-01, 7.7903E-03, -7.4262E-05, # AN 10 8.6475E-01, 3.7786E-03, -6.3711E-05, # AN 11 7.6261E-01, 1.1277E-02, -1.3603E-04, # AN 12 9.5367E-01, 3.4600E-03, -6.8573E-05, # AN 13 8.2851E-01, 9.2871E-03, -1.3003E-04, # AN 14 1.0000E+00, 0.0000E+00, 0.0000E+00, # AN 15 no rx 7.8560E-01, 1.0037E-02, -1.1899E-04, # AN 16 8.5390E-01, 6.9920E-03, -9.7635E-05, # AN 17 na = avg 8.7915E-01, 8.0377E-03, -1.7129E-04, # AN 18 8.7289E-01, 7.1793E-03, -1.0484E-04, # AN 19 7.8001E-01, 9.2160E-03, -9.7229E-05, # AN 20 1.2366E+00, -6.5267E-03, -2.5926E-05, # AN 21 7.5337E-01, 1.1281E-02, -1.3046E-04, # AN 22 8.5390E-01, 6.9920E-03, -9.7635E-05, # AN 23 na = avg 1.0165E+00, -7.3935E-04, -3.9395E-05, # AN 24 6.7155E-01, 1.3864E-02, -1.4794E-04, # AN 25 1.0983E+00, -3.1293E-03, -3.4942E-05, # AN 26 7.9194E-01, 9.6426E-03, -1.1291E-04, # AN 27 9.1300E-01, 5.5659E-03, -9.4216E-05, # AN 28 5.8500E-01, 1.1650E-02, 6.4210E-05, # AN 29 8.5390E-01, 6.9920E-03, -9.7635E-05 ] # AN 30 = avg print ' Filled gain coefficients for VLA at Qband'; } else { if ( band == 'U' | band == 'u' ) { c := [ 1.0213E+00, -1.6223E-04, -1.5142E-05, # AN 1 9.9380E-01, 5.2939E-04, -1.1844E-05, # AN 2 9.5077E-01, 1.7396E-03, -1.5404E-05, # AN 3 9.3945E-01, 1.6803E-03, -9.2158E-06, # AN 4 9.6259E-01, 1.5072E-03, -1.5226E-05, # AN 5 9.7797E-01, 9.2114E-04, -9.6542E-06, # AN 6 9.5072E-01, 1.6461E-03, -1.3873E-05, # AN 7 9.6315E-01, 1.3102E-03, -1.1699E-05, # AN 8 9.8908E-01, 7.4323E-04, -1.2785E-05, # AN 9 9.4594E-01, 1.9071E-03, -1.6902E-05, # AN 10 9.7378E-01, 1.2212E-03, -1.4299E-05, # AN 11 9.2886E-01, 2.0099E-03, -1.4326E-05, # AN 12 9.7188E-01, 1.2060E-03, -1.2958E-05, # AN 13 9.5259E-01, 1.5981E-03, -1.3701E-05, # AN 14 9.9007E-01, 7.6713E-04, -1.5014E-05, # AN 15 9.3265E-01, 2.0155E-03, -1.5228E-05, # AN 16 9.4203E-01, 2.3229E-03, -2.3322E-05, # AN 17 9.8450E-01, 1.3243E-03, -2.9131E-05, # AN 18 9.4425E-01, 2.1415E-03, -2.0617E-05, # AN 19 9.3657E-01, 2.9488E-03, -3.4398E-05, # AN 20 9.9317E-01, 5.8221E-04, -1.4247E-05, # AN 21 9.4996E-01, 1.7295E-03, -1.5152E-05, # AN 22 9.7844E-01, 9.0703E-04, -9.8916E-06, # AN 23 9.7704E-01, 8.9316E-04, -9.0036E-06, # AN 24 9.1715E-01, 2.1045E-03, -1.3433E-05, # AN 25 9.9063E-01, 4.9646E-04, -6.6641E-06, # AN 26 9.5552E-01, 1.6396E-03, -1.5122E-05, # AN 27 9.8287E-01, 8.9881E-04, -1.1829E-05, # AN 28 1.0000E+00, 6.9500E-04, 0.0000E+00, # AN 29 9.6350E-01, 1.3945E-03, -1.5014E-05 ] # AN 30 = avg print ' Filled gain coefficients for VLA at Uband'; } else { c := array( [1,0,0],3,30 ); print ' Filled gain coefficients for VLA to unity'; } } } # reshape as array c::shape := [3,30]; return c; }