Date - 2001.12.12 Tester - S.T. Myers (AOC) Platform - Linux (sandrock) Note: bugs are denoted by lines with prefix: >>>BUG: and queries are denoted by lines with prefix: >>>QUERY: while comments are prefixed by: >>>COMMENT: Some enhancement requests are given after >>>REQUEST: Notes added later for clarification are indicated by >>>NOTE: ------------------------------------------------------------------------------- GOAL: Develop script for gain curve and opacity correction 1. Date - 2001.12.12 Version - Weekly (1.7 build 047) Directory - /home/sandrock/smyers/Testing/2001-12-12/ Data - 20010204AK.UVF Script - 2001-12-12-myers.g In light of the NAUG list discussions on the viewer, would try msplot() again and exercise the viewer more. In particular, Dave King suggested: I suggest using 'Power Cycles' (usually decreasing it -- lightens things up). It is the most useful control for making a large range of data visible, IMO, but is also probably the least used. It scales the data using power or log functions, before mapping to color. Try to find this... msplot() from GUI >>>BUG: Again this comes up in zoomed mode. I have the feeling that if I blow away the cache and everything else it would come up right. >>>NOTE: I could not make this repeat - bummer. Play with colormap Shift/Slope and Contrast - note that assigning the former to LeftMouse and latter to CenterMouse worked ok, but this could have been the assignment for a general Colormap control. Found the power cycles - it was just what I did before! Doh, thought there was something that I missed. Basically I cant seem to get these to work well. How does TVFLG do it? Those still seem better. Played around a bit and didn't do any better. Note that you might get better control by defining the transfer functions like is done in xv, but setting points on the wedge and controlling the transfer. There must be better ways out there... 2. Now back to the gain curve problem. We had previously from George the example recipe: # Example 7: measures: #--------------------- include 'measures.g' print vlapos:=dm.observatory('VLA'); print dm.measure(vlapos,'WGS84') # convert to geodetic dm.doframe(vlapos) print tim:=dm.epoch('utc','today') ; dm.doframe(tim) print srcdir:=dm.source('0234+285') print 'AZ =', dq.convert( dm.measure(srcdir,'AZEL').m0,'deg') print 'EL =', dq.convert( dm.measure(srcdir,'AZEL').m1,'deg') # should use getvalue (not .m0, etc.) above! What we will need to do is get some header info from the ms (like the observatory) and then go through each record, use the direction (from source) and time to determine the EL for that record, and apply opacity or gain (using ANT 1 and 2 for the baseline and some table of gain curves versus band). >>>REQUEST: Note that it would be a great help if EL was a field in the ms (or easily derived on-the-fly like HA is) as AZ-EL is commonly used (more so than HA!). (Submitted as AOCso03268) First, lets try and just get the table stuff and hardwire what we can. Select data based on FIELDID FLAG WEIGHT? FLAG_ROW? Will want to pull out the following columns for data that is not flagged: TIME ANTENNA1 ANTENNA2 CORRECTED_DATA Note - do we want to do this on the data = update the CORRECTED_DATA or work on a calibration table? For now let us do it on the data, I don't know what the interaction would be with the cal tables (would need to make a compatible 'g' cal table and apply it with calibrater - long term this is what should be done). One thing I am not sure about is that the CORRECTED_DATA column seems to be overwritten by calibrater, in which case I would need to make a table or to fudge the DATA column directly. Reading the description of calibration in Getting Results Vol2 - 1.2.2 suggests that indeed the CORRECTED_DATA column is computed by multiplying all the antenna-based Jones matrices, e.g. Vcor_ij = B_ij*G_ij*D_ij*P_ij*T_ij*V_ij where G_ij = G_i x G_j (outer product). Note that T_ij is the tropospheric calibration matrix we want for the opacity correction (actually, the "T" matrix in calibrater is used for gain variations that are common to polarizations, like in AIPS CALIB where you set APARM(3)=1 to average RR and LL, so never mind). Talked to George about this and he confirmed that I will need to fudge the DATA column itself, as doing a cal table might be complicated. Note that currently the calibration is not incremental (there is no equivalent of SN tables) but that they are working on being able to increment the G calibrations for example so an a-priori atmosphere and antenna gain G table could be worked in to this scheme. For now just muck with the data column. George supplied the following example he sent to Joe Lazio to copy DATA columns around (to get around the non-incremental thing): # To copy all or some of CORRECTED_DATA to DATA column # gmoellen (01Nov30) # (use after running calibrater.correct()) # (be careful if you dataset is very large!) include 'table.g'; # make table tool t:=table('yourms.ms',readonly=F); # extract existing data columns origdata:=t.getcol('DATA'); caldata:=t.getcol('CORRECTED_DATA'); # Make column that has same properties as DATA column: newcol.name:='ORIG_DATA'; newcol.desc:=t.getcoldesc('DATA'); t.addcols(newcol) # Put origdata into new column: t.putcol('ORIG_DATA',origdata); # Now put calibrated data into DATA column: t.putcol('DATA',caldata); # ALTERNATIVELY, if you only want to replace certain FIELDs' data: # make a mask that selects only the fields you want to treat: flds:=t.getcol('FIELD_ID')+1; # (+1 corrects c++ indexing) fmask:= (flds==3 || flds==4) # select fieldids 3 & 4 # put old data into caldata except for selected fields: caldata[,,!fmask]:=origdata[,,!fmask]; t.putcol('DATA',caldata); # that should do it Note that tableiterator should be used in the column operations to cut down the amount of data that needs to be kept in glish. Thus, the best approach is to extract the data for a single timestamp, work out the source (should be one source if single subarray, but might eventually do this in general) and thus the elevation, then for each row get ANTENN1 and ANTENN2 derive EL gain correction and opacity correction, multiply data by this and replace in table. The drawback is you would have to revert the dataset if you did this wrong. Good enough for now. Start coding: include 'measures.g'; include 'table.g'; include 'quanta.g'; # # Define the opacity correction (would use a par.tau0) # 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 ) { # Example: # ANTEN 2 0; STOKES = 'L'; BIF=1; EIF=1; # CLCOR 1.0041e+00, -6.0815e-04, 2.2647e-05, 0; # Could use polyfitter.eval function (in mathematics module) but # is overkill for this c := [1.0041e+00, -6.0815e-04, 2.2647e-05 ]; g := c[1] + c[2]*el + c[3]*el*el; g[ g < 1.e-4 ] := 1.e-4; return g; } >>>NOTE: The statement g[ g < 1.e-4 ] := 1.e-4; is fixed from the original (ncorrect) version. Now try to test path Start from UVF par.msf := '20010204AK.ms'; par.tau0 := 0.10; - makemsf(par) msf := par.msf; tau0 := par.tau0; # # 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); fdirs := ftab.getcol("REFERENCE_DIR"); fields := ftab.getcol("NAME"); ftab.close(); >>>BUG: Note - on debugging, I did a ftab.getcolkeyword on a keyword which did not exist, and this seemed to kill the ftab table object. This shouldn't happen. - fdref := ftab.getcolkeyword("REFERENCE_DIR","direction"); : get_keyword: Keyword direction in column REFERENCE_DIR does not exist File: tableserver.g, Line 71 Stack: .(), tableserver.g line 953 .(), table.g line 844 .() - fdirs := ftab.getcol("REFERENCE_DIR"); : get_column: Illegal table-id sent (no such table) File: tableserver.g, Line 71 Stack: .(), tableserver.g line 739 .(), table.g line 762 .() - ftab.close(); : close: Illegal table-id sent (no such table) File: tableserver.g, Line 71 Stack: .(), tableserver.g line 297 .(), table.g line 478 is_fail(), table.g line 478 .() >>>NOTE: not sure I want to try to make this repeat. If I get this again I will bug() it. # # 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('j2000',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 72] ] # nvis := st.nrows(); vis := st.getcol("DATA"); # # Get AN number vectors # an1 := st.getcol("ANTENNA1")+1; # need to convert these to real ANs an2 := st.getcol("ANTENNA2")+1; # # Construct el numeric vector # elv := array( el.value, nvis ); # # construct opacity corrections # oc := opac( tau0, elv ); # # construct gain corrections # gc := gaincv( an1, an2, elv ); vc := oc/gc; # # make compatible with vis vector (I dont know how to do this better) # 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); } ti.done(); mst.close(); Put this in as a function gaincor(par) in 2001-12-12-myers.g par.msf := '20010204AK.ms'; par.tau0 := 0.10; - makemsf(par) Data loaded into 20010204AK.ms F - gaincor(par) Data is from telescope VLA Found 8 fields Processing timestamp 2001/02/04/08:27:30.000 Corrected 1159+292 with 72 rows at EL = 67.26 Processing timestamp 2001/02/04/08:27:33.334 Corrected 1159+292 with 156 rows at EL = 67.27 Processing timestamp 2001/02/04/08:27:36.666 Corrected 1159+292 with 380 rows at EL = 67.28 Processing timestamp 2001/02/04/08:27:40.001 Corrected 1159+292 with 380 rows at EL = 67.29 ... Processing timestamp 2001/02/04/10:41:06.668 Corrected 1159+292 with 650 rows at EL = 82.24 Processing timestamp 2001/02/04/10:41:10.000 Corrected 1159+292 with 650 rows at EL = 82.23 Completed gain and opacity correction OK - finally seems to work! Took a bit less than a minute (its not very efficient). Next need to check that the corrections being applied are sensible. Pundits could probably point out where I can be more efficient. One thing I need to know is how to extract the epoch from the FIELD column keywords or in the header somewhere (Im assuming J2000). In the painful course of trying to figure out how to do the measures and quanta and table gymnastics, it became clear that examples like this will need to be included in the documentation as well as in recipes to illustrate how to do stuff with real data, not just toy examples like manipulating the 'today' epoch and inputting positions by hand. Sigh. I'll call this a victory for now... now need to get the real antenna based gain curves from Bryan, put those in somewhere, and also figure out how to get the real antenna numbers (which the gain curves will be indexed by). See the attached script 2001-12-12-myers.g for the final function. Also will eventually go to the scriptorum. ------------------------------------------------------------------------------ Testing logs archived at http://www.aoc.nrao.edu/~smyers/aips++/testlogs/ This file can be found as 2001-12-12-myers.txt The script can be found as 2001-12-12-myers.g