;+ ; NAME: ; NEWSCAN ; PURPOSE: ; Procedure to be called on starting a new scan. Reads header and other ; segments and sets up variables as needed. ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm[,geometry,pRecent,pOld,cycle]$ ; [,pshift=pshift][,xpol=xpol]) ; INPUTS: ; a A file associated variable as returned from OPENARC. ; nrec The number of records in the file, as returned from OPENARC ; hrec The record number of the header record for the desired scan ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; pshift A two-element array giving the amount to shift the phase ; center, in arcsec, as [x,y], with x positive W and ; y positive N. If omitted, [0,0] is assumed. ; xpol A switch that states whether cross-polarization measurements ; (i.e. RL or LR) are analyzed. If omitted or set to zero, ; cross-polarization measurements (if any) are ignored. If ; set to 1, parallel-polarization measurements are ignored ; and the crossed-polarization measurements RL are placed ; where RR would go, and LR are placed where LL would go. This ; is implemented by changing the OBSEQ structure. ; ROUTINES CALLED: ; getdata, tl_decode, decode, findseg, getparms, poly_fit, get_index, ; ovsa_lun_assoc ; OUTPUTS: ; rec The record number of the next data record, or -1 if an error. ; header A structure containing the information in the scan HEADER segment ; obseq A structure containing the information in the OBSEQ segment ; cfg A structure containing the information in the CONFIG segment ; traj A structure containing the information in the TRAJECTORY segment ; refcal A structure containing the baseline and total power calibration data ; gparm A structure containing the GAINCAL parameters ; geometry A structure containing the baseline geometry, source params, etc. ; pRecent A pointer to the full-resolution data ; pOld A pointer to time-averaged data (not used--not defined) ; cycle A structure containing pointers to key tables used in the interpretation ; of one cycle of data ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 01-Sep-1998 by Dale E. Gary ; 03-Nov-1998 DG ; Fixed problem with SOLAID hanging in real-time mode? Inserted ; suitable wait for CONFIG and OBSEQ segments--up to 10 seconds ; 09-Nov-1998 DG ; Add wait for trajectory as well. ; 20-Dec-1998 DG ; Added TRAJ segment return argument and removed redundant SRCNAM ; argument. Philosophy is that NEWSCAN can be called in a quite ; general case of non-solar and non-pcal data, with arguments ; excluding GEOMETRY, pRECENT, pOLD and CYCLE. Converted NEWSCAN ; to a function, which returns the record number of the first ; data record, or -1 if an error. ; 07-Feb-1999 DG ; Implement decoding/accessing of EPHEMERIS segment and related ; issues. ; 14-Feb-1999 DG ; Fixed a couple of problems that surfaced when no TRAJ segment was ; found in data file. ; 15-Mar-1999 DG ; Added code for DLASCAN scan type, and also made TEST scan type a ; little more clever. ; 26-Jul-1999 DG ; Added handling for TPCAL, AMPCAL, and PHZCAL segments, which includes ; a new argument REFCAL made up of data from these three segments. ; 27-Jul-1999 DG ; Finished apportionment of ref calibration data into REFCAL structure. ; 09-Sep-1999 DG ; Fixed small bug that caused a crash when more than one calibrator ; entry had a given source name. Now issues appropriate error. ; 13-Sep-1999 DG ; Changed to reflect changes to TPCAL segment, which has been enlarged ; to 2 records. The code can handle either 1- or 2-record versions for ; now. ; 13-Nov-1999 DG ; REFCAL factor calculation made independent of number of antennas/baselines ; in the data. Also other changes for independence on active antennas. ; 24-Nov-1999 DG ; Add reading of GAINPARM segment ; 02-Dec-1999 DG ; Changes to REFCAL structure to allow inclusion of update parameters ; separately from primary parameters. This allows optional application of ; update parameters. ; 11-Jan-2000 DG ; Eliminated hardwired directory locations. ; 03-Feb-2000 DG ; Simplified code by moving INDEX finding to a new routine, GET_INDEX ; Also updated reading of the TPUPDATE segment, to read from the file using ; the INDEX record if available. ; 04-Feb-2000 DG ; Skip dlcal file handling on ioerror ; 13-Mar-2000 DG ; Remove reading of PLD file for "offsun" measurement, use GAINPARM entries ; instead. This removes a dependence on number of antennas. ; 14-Mar-2000 DG ; Slight change to make RA and DEC of celestial sources slightly more ; correct by bypassing POLY_FIT routine in that case. This probably has ; no real consequence. Added DUT (from EPHEM structure) to GEOMETRY. ; 23-Mar-2000 DG ; Fix bugs related to use of pre-GAINPARM-segment data (replace rec eq -1 ; with rec le 0 in several places, and set gcalflg when old-style gcal ; files are being used). ; 17-Apr-2000 DG ; Fixed bug with GPARM segment (T/L field was not removed from records ; 2 and 3 before decoding). ; 23-May-2000 DG ; Change GEOMETRY structure to anonymous so that change in number of ; antennas does not cause a problem. ; 08-Jun-2000 DG ; Implement handling of DAILYPHZ segment and associated REFCAL PhaseUpd ; NB: The REFCAL structure has changed due to these changes. ; 16-Jul-2000 DG ; Tracked down and fixed long-standing bug in BC phases--a rather ; complicated issue relating to phase closure to -135 degrees. Also ; changed the PhaseUpd definition and handling. ; 23-Jul-2000 DG ; Fixed a similar problem with the PhaseUpd, and verified that it is ; correct. ; 28-Jul-2000 DG ; Slight change to make DAILYPHZ reading work with older-style files ; (those with DAILYPHZ segments written at the end of the file). ; 23-Sep-2000 DG ; Remove restriction of 5 MB for pRecent size (since this gives ; less than 1 hour of data for our normal mode, now, at 10 ms sampint) ; 07-Oct-2000 DG ; Fixed phase closure of BC baseline to reflect new sign convention ; implemented in GAINCOR on 26-Sep-2000. This was just a reversal of ; sign when cfg.mux = -1, in the REFCAL section. ; 18-Mar-2001 DG ; I find that reversal of the BC baseline sign in the phase update ; parameters of the REFCAL section was NOT correct, so I removed it. ; 20-Mar-2001 DG ; Added solar parameters to GEOMETRY structure, from EPHEM.SOLAR. ; 29-Mar-2001 DG ; Corrected calculation of amplitude RMS factors for REFCAL structure. ; Changed from ds = sqrt(da^2 + db^2) to ds = sqrt[(da*b)^2 + (a*db)^2] ; 03-Apr-2001 DG ; Added call to GET_TIMING and added MSOFF to cycle structure ; 18-Jul-2001 GN ; Replaced get_index by ovsa_get_index to avoid conflict with Yohkoh ; 15-Jan-2002 DG ; Slight changes to work with new SATELLITE more (determined from ; TRAJ.SRCNAME). Main change is to add SRCNAME to GEOMETRY structure ; so that GAINCOR can skip phase correction if SRCNAME='SATELLITE' ; 06-Apr-2002 DG ; Fixed a bug that never became important until new 426 mux sequence ; was used. The POL indexing used ANT2 instead of IANTS[ANT2]. ; 08-Jun-2002 DG ; Added PSHIFT keyword to allow shifting the phase center. ; 08-Jul-2002 DG ; Put in a "kludge" to force a scancode of TEST if the source is a ; satellite. ; 14-Dec-2002 DG ; Applied PSHIFT to TRAJ structure itself, so that shift will be ; propagated to mapping stage. ; 23-Dec-2002 DG ; Add code for handling special MUX mode 17 (GAINPARM swap). ; 03-Jan-2003 DG ; Made (final?) change to REFCAL section to convert antenna-based ; to baseline-based phase. ; 04-Jan-2003 DG ; Fixed the above change due to discovery of a conceptual error. ; 11-Jan-2003 DG ; Went to a slightly different scheme in which AB channel offset is ; given, instead of making BC and CA relative to AB. ; 26-Jul-2003 DG ; Added more diagnostic information to error message when CONFIG ; segment is not found in time. ; 04-Aug-2003 DG ; I learned that the problem is that NREC is not changed as the ; file grows past 10000 records. I added code to determine the ; new NREC, at least when the CONFIG segment is being sought. The ; same error should just occur later, but we will see. ; 17-Aug-2003 DG ; New version of IDL (6.0) required iant and jant to be scalar, so ; fixed that. ; 07-Nov-2003 DG ; Added XPOL keyword and associated code to determine whether to ; analyze cross-polarization measurements. ; 14-Nov-2003 DG ; Change GPARM part of code to read 6 records after 2003.310, to ; use data from new dual-polarization front ends. ; 04-Dec-2003 DG ; Fix bug in dual-polarization code, to check nrr, nrl, nlr, nll ; before using corresponding indexes. ; 10-Jan-2004 DG ; Dual polarization scheme changed on 2003 Nov 18, so that obseq ; polarization code is now 3-bit corresponding to antennas n21. ; A set bit means LCP, a 0 bit is RCP. Changed code to deal with ; cross-polarization in the new scheme. ; 04-Apr-2004 DG ; Change setting of APOL to reflect new dual polarization feeds ; on ants 5 and 6. ; 07-Aug-2004 DG ; Coded reading of GPS satellite ephemeris into GEOMETRY structure. ; 11-Sep-2004 DG ; Change setting of APOL to reflect new dual polarization feeds ; on ants 7 and 8. ; 26-Sep-2004 DG ; Added code to regularize handling of polarization for data taken ; after all antennas are outfitted with dual polarization feeds. The ; handling of polarization after this time is much simpler both for ; XPOL = 0 (no cross polarization) and XPOL = 1 (cross polarization ; only) cases. New code is activated when header.tls.yrday > 2004.231 ; 07-Oct-2004 DG ; Standardize on 2 polarizations after 2004.231, even for data that ; were taken only in R or L. Missing polarization will be treated ; as NaN. ; 15-Mar-2005 DG ; Added code to deal with new satellite scheme, in which SATPOS lines ; are read from the dialog history and fit coefficients are entered ; into the GEOMETRY structure. This supersedes changes of 07-Aug-2004. ; 23-Aug-2005 DG ; I had the application of the ND_eff upside down on the OFFSUN adjustment. ; 27-Aug-2005 DG ; Updated to handle 4-record TPCAL and TPUPD segments. ; 07-Nov-2005 DG ; Fixed error in conversion of REFCAL phases to current MUX setting (and ; corresponding error in APCAL_ADD_2M_RL). Also apply CFG.BLDIR to phase slope ; in PHASEUPD. ; 20-Feb-2006 DG ; Updated to handle 4-record AMPCAL and 2-record PHZCAL segments. Also ; fixed code to work with R and L on all antennas. ; 29-Jun-2006 DG ; Remove all calls to WIDGET_MESSAGE and WSTR, to avoid software hanging up ; waiting for user input. Benign defaults assumed, or else error return. ; 18-Jul-2006 DG ; The geometry structure for GPS satellites was too short in some cases due to ; reliance on DIALOG HISTORY for SATPOS information. Changed to still read ; DIALOG HISTORY, but now checks for on-disk file. If it exists and has ; the same first line information, then the disk file is used instead. This ; approach allows the scheme to work both in real-time at OVRO and for later ; analysis. ; 19-Jun-2007 DG ; Fixed minor bug in reading the new dual polarization PHZCAL data. ;- function newscan,a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cycle,$ pshift=pshift,xpol=xpol common dlabit,afac,bfac ; Find and return the INDEX segment as a structure. If no INDEX segment ; exists, the user is prompted for the option to create a new one from the ; information in the file. IDX is an empty INDEX structure if the user ; does not choose to create a new one. idx = ovsa_get_index(a,nrec,/create) ; -------------------------------- HEADER -------------------------------------- rec = hrec - 1 data = getdata(rec,a) ; Read the header tls = tl_decode(data) ; Decode the header segment header = decode(data,!SEGM.HEADER) nant = header.nant ; -------------------------------- CONFIG -------------------------------------- ; Find the configuration segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.CONFIG,idx=idx) try = 0 while (rec eq -1) do begin wait,3 if (nrec gt 10000) then begin ; There was a fatal flaw in a growing file with more than 10000 ; records, since nrec was always equal to the header record number. ; This code determines the current nrec, at the time of the attempt. alun = ovsa_lun_assoc(a) if (alun eq -1) then begin print,'NEWSCAN: Error reading ARC file!' return,-1 endif nrec = (fstat(alun)).size/2048L endif rec = findseg(a,hrec,nrec,!SEGM.CONFIG,idx=idx) try = try + 1 if (try gt 10) then begin stime = systime() print,'NEWSCAN: Error reading CONFIG segment after 10 tries (header rec='+string(hrec)+') at '+stime return,-1 endif endwhile ; Decode the config segment data = getdata(rec-1,a) cfg = decode(data,!SEGM.CONFIG,header=header) ; ----------------------------- EMPHEMERIS ------------------------------------- ; Since an EPHEM segment is written immediately after the CONFIG segment ; (if written at all) we do not have to wait for it ; Find the emphemeris segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.EPHEM,idx=idx) if (rec le 0) then begin ephem = get_eph_struct() ; Get an empty structure endif else begin data = getdata(rec-1,a) ephem = decode(data,!SEGM.EPHEM) endelse ; -------------------------------- TPCAL --------------------------------------- ; Since a TPCAL segment is written immediately after the EPHEM segment ; (if written at all) we do not have to wait for it ; Find the tpcal segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.TPCAL,idx=idx) if (rec le 0) then begin tpcal = get_tpcal_struct() ; Get an empty structure endif else begin data = getdata(rec-1,a) data2 = getdata(rec,a) ; TPCAL segment is two records long if (tls.yrday ge 2004.231) then begin ; TPCAL segment has changed to 4 records data3 = getdata(rec+1,a) data4 = getdata(rec+2,a) endif ; Interim code to handle old files with TPCAL segments of only 1 record tptls = tl_decode(data2) if (tptls.segmentcode eq !SEGM.TPCAL) then begin ; This one does have two records ; Simply concatenate the two records, after removing the TLS data ; from DATA2. data = [data[0:1015],data2[12:*]] endif ; Interim code to handle case of two records, when four are expected tptls = tl_decode(data3) if (tptls.segmentcode eq !SEGM.TPCAL) then begin ; This one does have three, so assume four, and make two dimensional data = reform([data,data3[0:1015],data4[12:*]],2028,2) endif tpcal = decode(data,!SEGM.TPCAL) endelse ; -------------------------------- TPUPDATE --------------------------------------- ; Find the tpupdate segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.TPUPD,idx=idx) if (rec le 0) then begin tpupd = get_tpcal_struct() ; Get an empty structure--same as TPCAL structure endif else begin data = getdata(rec-1,a) data2 = getdata(rec,a) ; TPCAL segment is two records long if (tls.yrday ge 2004.231) then begin ; TPUPD segment has changed to 4 records data3 = getdata(rec+1,a) data4 = getdata(rec+2,a) endif data = [data[0:1015],data2[12:*]] ; Interim code to handle case of two records, when four are expected tptls = tl_decode(data3) if (tptls.segmentcode eq !SEGM.TPUPD) then begin ; This one does have three, so assume four, and make two dimensional data = reform([data,data3[0:1015],data4[12:*]],2028,2) endif tpupd = decode(data,!SEGM.TPUPD) endelse ; ------------------------------- RCVRCAL --------------------------------------- ; Read the GCSP.PLD file from c:\working (this is temporary) ; gcal = get_tpcal_struct() ; Get an empty structure--same as TPCAL structure ; ON_IOERROR,skipgcal ; openr,gcalun,/get_lun,!defaults.workdir+'gcsp.pld' ; hed = strarr(16) ; readf,gcalun,hed ; gcaldat = fltarr(16,86) ; readf,gcalun,gcaldat ; icol = [1,4,7,10,13] ; ndeff = fltarr(5) ; for i = 0,4 do begin ; ndlist = reform(gcaldat(icol[i]+1,[7,9,14,25,41,53]-5)) ; ndeff[i] = 10^(total(alog10(ndlist))/6.0) ; endfor ; for i = 0, 85 do begin ; gcal.tpfac[[0,3,4,5,6],i] = reform(gcaldat[icol,i])*500./ndeff ; gcal.tprms[[0,3,4,5,6],i] = reform(gcaldat[icol+2,i])*500./ndeff ; gcal.tpfac[[1,2],i] = 0.0 ; gcal.tprms[[1,2],i] = !values.f_nan ; endfor ; free_lun,gcalun ;skipgcal: ; ------------------------------- AMPCAL --------------------------------------- ; Since an AMPCAL segment is written immediately after the TPCAL segment ; (if written at all) we do not have to wait for it ; Find the ampcal segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.AMPCAL,idx=idx) if (rec le 0) then begin ampcal = get_ampcal_struct() ; Get an empty structure endif else begin data = getdata(rec-1,a) data2 = getdata(rec,a) ; AMPCAL segment is two records long if (tls.yrday ge 2004.231) then begin ; AMPCAL segment has changed to 4 records data3 = getdata(rec+1,a) data4 = getdata(rec+2,a) endif ; Simply concatenate the two records, after removing the TLS data ; from DATA2. data = [data[0:1015],data2[12:*]] ; Interim code to handle case of two records, when four are expected amptls = tl_decode(data3) if (amptls.segmentcode eq !SEGM.AMPCAL) then begin ; This one does have three, so assume four, and make two dimensional data = reform([data,data3[0:1015],data4[12:*]],2028,2) endif ampcal = decode(data,!SEGM.AMPCAL) endelse ; ------------------------------- PHZCAL --------------------------------------- ; Since an PHZCAL segment is written immediately after the AMPCAL segment ; (if written at all) we do not have to wait for it ; Find the phzcal segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.PHZCAL,idx=idx) if (rec le 0) then begin phzcal = get_phzcal_struct() ; Get an empty structure endif else begin data = getdata(rec-1,a) if (tls.yrday ge 2004.231) then begin ; PHZCAL segment has changed to 2 records data2 = getdata(rec,a) endif phztls = tl_decode(data2) if (phztls.segmentcode eq !SEGM.PHZCAL) then begin ; This one does have two, so make two dimensional data = reform([data,data2],1024,2) endif phzcal = decode(data,!SEGM.PHZCAL) endelse ; ------------------------------- DAILYPHZ --------------------------------------- ; If this is a PCAL scan, read its DAILYPHZ segment if any. If it is ; not a PCAL scan, read the DAILYPHZ segments from the nearest two PCAL ; scans before and after this scan, and do a linear interpolation between ; them. If either the earlier or later DAILYPHZ segment does not exist, ; simply use the nearest DAILYPHZ segment. npcal = 0 if (idx.nscans ne 0) then begin ; Case of existing INDEX record ipcal = where((*idx.pscan).scancode eq !SCAN.PCAL,npcal) if (npcal ne 0) then srec = (*idx.pscan)[ipcal].srec endif else begin ; The hard way, when no INDEX record is available seglist = get_segments(a,nrec) ; Get indexes of PCAL header lines (lines containing string 'Header:') ipcal = where(strpos(seglist,'Header: PHASE CAL') ne -1,npcal) if (npcal ne 0) then srec = uint(strmid(seglist[ipcal],13,5)) endelse ; At this point, srec is an array containing the start record numbers ; for each PCAL scan, or is undefined if there are no PCAL scans. ptref = [0L,0L] ; Reference time from the DAILYPHZ segments. ; Unless NPCAL = 0, we have a list of PCAL scan header record numbers, so ; now we can read the DAILYPHZ segments from them (should immediately follow ; the headers if they exist). if (npcal ne 0) then begin ; DP will contain the two DAILYPHZ segments on either side of current scan dailyphz = get_daily_struct() dp = replicate(dailyphz,2) ; Loop over PCALs, which are in ascending time order for i = 0, npcal-1 do begin data = getdata(srec[i],a) ; Gets the record after the header record dailyphz = decode(data,!SEGM.DAILYPHZ) if (size(dailyphz,/type) ne 8) then begin ; The scan after the PCAL header was not a DAILYPHZ, so see if ; the DAILYPHZ segments are at the end of the file... This is ; a bit of a kluge. If files of this type, the last NPCAL segments ; should be the DAILYPHZ segments corresponding to the NPCAL scans. data = getdata(nrec-npcal+i,a) dailyphz = decode(data,!SEGM.DAILYPHZ) endif if (size(dailyphz,/type) eq 8) then begin ; If DAILYPHZ is a structure, this was a good DAILYPHZ segment. ; Check the start record to verify that it is consistent with the ; PCAL scan we are looking for: if (dailyphz.tls.nrc eq srec[i]) then begin ; This is the right one, so compare its reference time with the ; current scan time. If before, enter it into the first DP entry. ; If after, enter it into the second DP entry and exit. if ((dailyphz.tls.msec+dailyphz.tls.day*86400000D) lt $ (header.tls.msec +header.tls.day*86400000D)) then begin ; Got a good DAILYPHZ with reference time BEFORE this scan dp[0] = dailyphz ; The calculation below is going to fail in the highly unusual ; case that the PCAL scan and the current scan are in different ; years! Note, if on different days, ptref can be negative. ptref[0] = dailyphz.tls.msec + (header.tls.day-dailyphz.tls.day)*86400000L endif else begin ; Got a good DAILYPHZ with reference time AFTER this scan dp[1] = dailyphz ; The calculation below is going to fail in the highly unusual ; case that the PCAL scan and the current scan are in different ; years! ptref[1] = dailyphz.tls.msec + (header.tls.day-dailyphz.tls.day)*86400000L if (ptref[0] eq 0L) then begin ; There must not have been a good DAILYPHZ before the scan, ; so use this one in both entries of DP. dp[0] = dp[1] ptref[0] = ptref[1] endif goto, done endelse endif endif endfor ; We get here only if the loop ran out without finding a good DAILYPHZ ; after the scan, so use this one (if any) in both entries of DP dp[1] = dp[0] ptref[1] = ptref[0] endif done: ; -------------------------------- GAINPARM ---------------------------------- ; Find the gain parameter segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.GPARM,idx=idx) if (rec le 0) then begin gainparm = get_gparm_struct() ; Get an empty structure endif else begin ; Decode the gain parameter segment (3 records) if (header.tls.yrday gt 2003.310) then n_rec = 6 else n_rec = 3 data = getdata(rec-1,a) for i = 1, n_rec-1 do data = [data,(getdata(rec-1+i,a))[12:*]] gainparm = decode(data,!SEGM.GPARM) if (cfg.muxseqid eq 17) then begin ; In special MUX mode 17, for chanphz calibration, ; ants 1, 2 and 4 masquerage as 6, 7 and 8 (*gainparm.pgnp)[4:6] = (*gainparm.pgnp)[0:2] o = (*gainparm.poffsets) aboff1 = o[2,1] aboff2 = o[1,2] bcoff1 = o[1,0] bcoff2 = o[0,1] caoff1 = o[2,0] caoff2 = o[0,2] o[1,0] = (o[6,0] = bcoff1) o[0,1] = (o[0,6] = bcoff2) o[3,0] = (o[2,1] = (o[4,1] = (o[6,5] = aboff1))) o[0,3] = (o[1,2] = (o[1,4] = (o[5,6] = aboff2))) o[2,0] = (o[5,0] = (o[4,3] = caoff1)) o[0,2] = (o[0,5] = (o[3,4] = caoff2)) (*gainparm.poffsets) = o endif endelse ; If there was no GAIN PARAMETER segment, use the "old-style" routines to ; read in the gain parameters (offsets, nonlinearity, attenuation factors) gcalflg = 0 if (gainparm.epoch.yr eq 0) then begin gcalflg = 1 ; Set flag to indicate old-style routines were used gainparm = getparms(tls) endif ; ptr_free,hsk.gainparm ; hsk.gainparm = ptr_new(gainparm) ; -------------------------------- OBSEQ -------------------------------------- ; Find the observing sequence segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.OBSEQ,idx=idx) try = 0 while (rec eq -1) do begin wait,1 rec = findseg(a,hrec,nrec,!SEGM.OBSEQ,idx=idx) try = try + 1 if (try gt 10) then begin print,'NEWSCAN: Error reading OBSEQ segment after 10 tries' ; ans = widget_message('NEWSCAN: Error reading OBSEQ segment',/ERROR) return,-1 endif endwhile ; Decode the observing sequence data = getdata(rec-1,a) obseq = decode(data,!SEGM.OBSEQ) if (header.tls.yrday ge 2003.310 and header.tls.scancode eq !SCAN.PEAKUP) then begin ; If this is a PEAKUP scan, set the polarization to 0 (RR). *obseq.ppol = 0 endif if (header.tls.yrday ge 2003.262 and header.tls.yrday lt 2003.322) then begin ; Code to handle special circumstance for cross-polarization rridx = where(*obseq.ppol eq 0,nrr) ; Both Ant 1 and Ant 2 are R lridx = where(*obseq.ppol eq 3,nlr) ; Ant 1 is L and Ant 2 is R rlidx = where(*obseq.ppol eq 12,nrl) ; Ant 1 is R and Ant 2 is L llidx = where(*obseq.ppol eq 15,nll) ; Both Ant 1 and Ant 2 are L endif else if (header.tls.yrday ge 2003.322) then begin ; Code to handle special circumstance for cross-polarization rridx = where(*obseq.ppol eq 0,nrr) ; Both Ant 1 and Ant 2 are R lridx = where(*obseq.ppol eq 1,nlr) ; Ant 1 is L and Ant 2 is R rlidx = where(*obseq.ppol eq 2,nrl) ; Ant 1 is R and Ant 2 is L llidx = where(*obseq.ppol eq 3,nll) ; Both Ant 1 and Ant 2 are L endif ; If not set, ignore the cross-polarized measurements if (not keyword_set(xpol)) then xpol = 0 if (header.tls.yrday ge 2003.262 and header.tls.yrday lt 2004.231) then begin ; For the period when some antennas have dual polarization front ends, while ; others do not, we actually change the OBSEQ depending on the value of XPOL ; to cause the data to be placed in the correct polarization channel. if (nrl gt 0 or nlr gt 0) then begin ; The XPOL keyword specifies what to do with cross-polarization ; measurements. If XPOL=0, then the cross-polarization measurements ; are flagged and not used. If XPOL=1, then any non-crossed ; measurements are flagged and the crossed measurements are ; buggered to look like non-crossed measurements. if (xpol eq 0) then begin ; Flag the cross-polarization measurements by replacing their HARM ; entry with -1 if(nrl gt 0) then (*obseq.pharm)[rlidx] = -1 if(nlr gt 0) then (*obseq.pharm)[lridx] = -1 if(nrl gt 0) then (*obseq.pnd)[rlidx] = 0 if(nlr gt 0) then (*obseq.pnd)[lridx] = 0 endif else if (xpol eq 1) then begin ; Flag the non-crossed polarization measurements by replacing their HARM ; entry with -1 if(nrr gt 0) then (*obseq.pharm)[rridx] = -1 if(nll gt 0) then (*obseq.pharm)[llidx] = -1 if(nrr gt 0) then (*obseq.pnd)[rridx] = 0 if(nll gt 0) then (*obseq.pnd)[llidx] = 0 ; Change the POL entries of the RL measurements to masquerade as RR if(nrl gt 0) then (*obseq.ppol)[rlidx] = 0 ; Change the POL entries of the LR measurements to masquerade as LL if(nrl gt 0) then (*obseq.ppol)[lridx] = 15 endif endif endif oldrec = rec ; ------------------------------ TRAJECTORY ------------------------------------- ; Find the trajectory segment record that corresponds to this scan rec = findseg(a,hrec,nrec,!SEGM.TRAJ,idx=idx) try = 10 while (rec eq -1) do begin wait,1 rec = findseg(a,hrec,nrec,!SEGM.TRAJ,idx=idx) try = try + 1 if (try gt 10) then goto,continue endwhile continue: if (try le 10) then begin ; Decode the trajectory segment, if any data = getdata(rec-1,a) traj = decode(data,!SEGM.TRAJ) srcname = traj.srcname ; PSHIFT is the phase center shift in arcsec, as a 2-element ; array, [x,y], with x positive W, y positive N. if (keyword_set(pshift)) then begin (*traj.pha)[0] = (*traj.pha)[0] + pshift[0]/3.6 (*traj.pdec)[0] = (*traj.pdec)[0] + pshift[1]/3.6 endif hao = ((*traj.pha)(0))/1000. deco = ((*traj.pdec)(0))/1000. endif else begin ; Return an empty TRAJECTORY structure (SRCNAME entry may be ; filled in below). traj = get_traj_struct() rec = oldrec endelse ; --------------------------- REFCAL structure -------------------------------- ; All of the contents of the reference calibrations are now in hand (even if ; just empty structures), so build the REFCAL structure from the contents ; of the calibration, but pared down to the minimum needed for this scan ; First adjust the 2m amplitude calibration factors to the 2m total power ; calibration adj2mcal,ampcal,tpcal,tpupd refcal = get_refcal_struct() refcal.nant = nant refcal.nbl = nant*(nant-1)/2 refcal.npol = 3 refcal.tls27m = tpcal.tls27m refcal.tls02m = tpcal.tls02m refcal.tlsi = ampcal.tlsi refcal.tlsrl = ampcal.tlsrl ; Free current pointers in the REFCAL structure ptr_free,refcal.pfrq ptr_free,refcal.pfactors ptr_free,refcal.prms ptr_free,refcal.pphase ptr_free,refcal.pfacUpd ptr_free,refcal.prmsUpd ptr_free,refcal.phaseUpd.ppoff ptr_free,refcal.phaseUpd.pdpdf ptr_free,refcal.poffsun ptr_free,refcal.pbldir ptr_free,refcal.pchanphz refcal.pbldir = ptr_new(intarr(refcal.nbl)) refcal.pchanphz = ptr_new(fltarr(refcal.nbl)) ; Fill in the frequency list refcal.pfrq = ptr_new((*obseq.phord)*0.2) ih = ((*obseq.phord)-5)>0 ; Indexes into harmonic array corresponding to h in obseq nf = n_elements(*obseq.phord) ; Fill in the factor list with NaN, then enter real values where appropriate factors = fltarr(refcal.npol,refcal.nant+refcal.nbl,nf) + !values.f_nan facupd = factors ; Fill in the rms list with NaN, then enter real values where appropriate rms = fltarr(refcal.npol,refcal.nant+refcal.nbl,nf) + !values.f_nan rmsupd = rms ; Fill in the phase list with NaN, then enter real values where appropriate phase = fltarr(refcal.npol,refcal.nbl,nf) + !values.f_nan ; Fill in the PhaseUpd entries with zeroes, then enter real values as approp. phaseupd = replicate({pupd},2) phaseupd[0].pdpdf = ptr_new(fltarr(refcal.nbl)) phaseupd[0].ppoff = ptr_new(fltarr(refcal.nbl)) phaseupd[0].tref = ptref[0] ; Value set during reading of DAILYPHZ segment (may be 0) phaseupd[1].pdpdf = ptr_new(fltarr(refcal.nbl)) phaseupd[1].ppoff = ptr_new(fltarr(refcal.nbl)) phaseupd[1].tref = ptref[1] ; Value set during reading of DAILYPHZ segment (may be 0) ; Fill in the offsun list with NaN, then enter real values where appropriate offsun = fltarr(refcal.npol,refcal.nant,nf) + !values.f_nan ; Get the effective noise diode values for each antenna from the GAINPARM ; structure. If any of the ND_eff values are zero (generally means there is no ; measurement) then set them to 1.0 as a nominal ND_eff value. nd_eff = (*gainparm.pgnp).ndfac bad = where(nd_eff eq 0,nbad) if (nbad gt 0) then nd_eff[bad] = 1.0 aa = header.aatab(where(header.aatab ne 255)) ; Active antennas if (header.tls.yrday ge 2004.231) then begin ; Translation from ipol,iant to feed number (just repeat L in third column) feedtrans = [[ 0, 1, 1],$ ; Antenna 1 (Lin, LCP) [ 2, 3, 3],$ ; Antenna 2 (RCP, Lin) [ 4, 4, 4],$ ; Antenna 4 (Lin) [ 5, 5, 5],$ ; Antenna 5 (Lin) [ 6, 6, 6],$ ; Antenna 6 (Lin) [ 7, 7, 7],$ ; Antenna 7 (Lin) [ 8, 8, 8]] ; Antenna 8 (Lin) endif else begin ; Translation from ipol,iant to feed number feedtrans = [[ 0, 1, 0],$ ; Antenna 1 (Lin, LCP) [ 2, 3, 3],$ ; Antenna 2 (RCP, Lin) [ 4, 4, 4],$ ; Antenna 4 (Lin) [ 5, 5, 5],$ ; Antenna 5 (Lin) [ 6, 6, 6],$ ; Antenna 6 (Lin) [ 7, 7, 7],$ ; Antenna 7 (Lin) [ 8, 8, 8]] ; Antenna 8 (Lin) endelse antlist = [ 1, 2, 4, 5, 6, 7, 8] ; Translation from antenna number to index ; Fill in baseline direction and backend channel phase offset from CFG segment ibl = -1 for m = 0, nant-1 do begin for n = m, nant-1 do begin ; iant and jant are antenna indexes corresponding to indexes m, n iant = (where(antlist eq aa[m]))[0] jant = (where(antlist eq aa[n]))[0] if (m ne n) then begin ibl = ibl+1 (*refcal.pBLDir)[ibl] = cfg.BLDir[m,n] (*refcal.pChanPhz)[ibl] = cfg.ChanPhz[m,n] endif endfor endfor ; Read the channel phases from cfg.chanphz for i = 0, 7 do cfg.muxarr[i,i] = -2 ; Flag non-baseline (total power) AB = where(cfg.muxarr eq 0,nab) if (nab eq 0) then begin phi_ab = 0. print,'NEWSCAN: Warning, file contains no AB baseline?' endif else phi_ab = cfg.chanphz[ab[0]]/!dtor BC = where(cfg.muxarr eq -1,nbc) if (nbc eq 0) then begin phi_bc = 0. print,'NEWSCAN: Warning, file contains no BC baseline?' endif else phi_bc = cfg.chanphz[bc[0]]/!dtor CA = where(cfg.muxarr eq 1,nca) if (nca eq 0) then begin phi_ca = 0. print,'NEWSCAN: Warning, file contains no CA baseline?' endif else phi_ca = cfg.chanphz[ca[0]]/!dtor pol = intarr(3) ; Loop over frequencies for i = 0, nf-1 do begin j = ih[i] ; Loop over active antennas ibl = -1 for m = 0, nant-1 do begin for n = m, nant-1 do begin ; iant and jant are antenna indexes corresponding to indexes m, n iant = (where(antlist eq aa[m]))[0] jant = (where(antlist eq aa[n]))[0] if (m ne n) then ibl = ibl+1 ; Determine polarization states for this pair of antennas for kpol = 0, 1 do begin ipol = cfg.antpol[kpol,iant] jpol = cfg.antpol[kpol,jant] pol[kpol] = (ipol0 ne 0,nksc) ; indexes of lines preceding current if (nsc gt 0) then begin psc = isc[ksc[nksc-1]] ; this is the index of the nearest preceding SATCLR line psc1 = isc[(ksc[nksc-1]+1)<(nsc-1)] ; index of next following SATCLR, if any endif else begin print,'NEWSCAN: Dialog History contains no lines for current day!' print,'NEWSCAN: Cannot unwind GPS satellite phase.' goto, skip_sat endelse ; Restrict lines to those between the preceding SATCLR and the next, if any if (psc1 eq psc) then dialog = dialog[psc+1:*] else dialog = dialog[psc+1:psc1-1] isp = where(strpos(dialog,'SATPOS') ne -1,nsp) ; Lines with SATPOS if (nsp ne 0) then begin ; Keep only lines that are contiguous (this should be unnecessary, but safe) bad = where((isp-shift(isp,1))[1:*] ne 1,nbad) if (nbad ne 0) then isp = isp[0:bad[0]] ; Whew. We finally have the good, contiguous lines. ; Eliminate leading part of each line (part with time and date) nsp = n_elements(isp) dd = 0 & hhmm = 0 & ha = 0. & dec = 0. for k = 0, nsp-1 do dialog[k] = strmid(dialog[k],17) ; See if a disk file exists, and compare dialog history with it. THis is to ; solve a problem with dialog history being incomplete. If the disk file ; contains the same information, read and keep the disk file... satposfile = !defaults.workdir+'satpos.fth' newdialog = '' if (file_test(satposfile)) then begin openr,satlun,/get,satposfile filedat = '' readf,satlun,newdialog while (not eof(satlun)) do begin readf,satlun,filedat newdialog = [newdialog,filedat] endwhile free_lun,satlun ; Eliminate first line (comment) if (n_elements(newdialog) gt 1) then newdialog = newdialog[1:*] ndd = 0 & nhhmm = 0 & nha = 0. & ndec = 0. reads,newdialog[0],ndd,nhhmm,nha,ndec reads,dialog[0],dd,hhmm,ha,dec ; If the first line time and coords agree, replace DIALOG with NEWDIALOG if (dd eq ndd and hhmm eq nhhmm and ha eq nha and ndec eq dec) then begin dialog = newdialog nsp = n_elements(dialog) endif endif tsat = (hasat = (decsat = dblarr(nsp))) sattab = dblarr(3,nsp,3) for k = 0, nsp-1 do begin reads,dialog[k],dd,hhmm,ha,dec sattab[0,k,0] = hms2sec(hhmm*100.)*1000d ; Time in ms ; Shift by phase shift amount if (keyword_set(pshift)) then ha = ha + pshift[0]/3600. if (keyword_set(pshift)) then dec = dec + pshift[1]/3600. hasat[k] = ha*3600000d/15. ; Convert degrees to ms decsat[k] = dec*3600000d ; Convert degrees to masec endfor traj.srcname = 'GPS ' for k = 0, nsp-3 do begin t = [0d,300000d,600000d] ha = hasat[k:k+2] dec = decsat[k:k+2] sattab[1,k,*] = poly_fit(t,ha,2) ; HA fit coefficients sattab[2,k,*] = poly_fit(t,dec,2) ; Dec fit coefficients endfor sattab[1,nsp-2:nsp-1,*] = sattab[1,[nsp-3,nsp-3],*] ; Repeat last two times in case of extrapolation sattab[2,nsp-2:nsp-1,*] = sattab[2,[nsp-3,nsp-3],*] ; Repeat last two times in case of extrapolation endif else begin print,'NEWSCAN: Dialog History contains no SATPOS line for scan!' print,'NEWSCAN: Cannot unwind GPS satellite phase.' goto, skip_sat endelse endif else begin print,'NEWSCAN: Dialog History contains no SATPOS lines!' print,'NEWSCAN: Cannot unwind GPS satellite phase.' goto, skip_sat endelse endif endelse skip_sat: endif ; Encapsulate it into a GEOMETRY structure for passing to BDOTS and DELAY routines. geometry = {bx:bx, by:by, bz:bz, htdiff:htdiff, a_ra:a_ra, a_dec:a_dec, psattab:ptr_new(sattab),$ year:tls.year, doy:doy, dut:ephem.dut, solar:ephem.solar, srcname:traj.srcname} *geometry.psattab = sattab ; Determine the size of the data arrays needed to hold RECENT and OLD data nchan = nant^2 nfrq = n_elements(*obseq.phord) npol = n_elements(uniq(*obseq.ppol,sort(*obseq.ppol))) ; NPOL = 2 means R and L data, but such data also includes I (linear) ; measurements, so there will really be 3 polarization states. if (npol eq 2 and header.tls.yrday lt 2004.231) then npol = 3 if (header.tls.yrday ge 2004.231) then npol = 2 ; After dual-poln front-ends installed ; ; Use the smaller of 1 hour's worth, or 5 Mbytes total data ; ntimes = long(5.e6/float(nchan*nfrq*npol*4)) dt = obseq.nx*cfg.sampintms ; ntimes = ntimes<(3600000L/dt) ntimes = 3600000L/dt navg = 120000L/dt ; Number of points in a 120 s average ; Declare storage via pointers, for the RECENT and OLD data array. pRecent = ptr_new(fltarr(nchan,nfrq,npol,ntimes)) pOld = ptr_new(fltarr(nchan,nfrq,npol,ntimes/navg)) ; Create pointer to Bigdata and Pol arrays. nndw = header.ndbw + header.nfgw + 1 ; Number of non-data words nws = header.nws ; Number of words per sample nsamp = obseq.nx ; Number of samples in one cycle pbigdata = ptr_new(fltarr(nws,nsamp)) ; Big enough to hold one cycle ; Obtain the polarization distribution apol = bytarr(nsamp,8) ; Create a poln list for each antenna pol = bytarr(nws,nsamp) ; Polarization list for each sample if (header.tls.yrday lt 2004.231) then begin ; Prior to full dual polarization on all antennas, the polarization situation was ; extremely complicated. This code tries to deal with all of the possibilities for i = 0, nsamp-1 do apol(i,*) = cfg.antpol(0,*) ; Set all polarizations to R apol(*,0) = cfg.antpol(*obseq.ppol and 1,0) ; Replace ant 1 pol with obseq value apol(*,1) = cfg.antpol(ishft(*obseq.ppol and 4,-2),1) ; Replace ant 2 pol with obseq value if (header.tls.yrday ge 2003.322) then begin apol[*,1] = cfg.antpol(ishft(*obseq.ppol and 2,-1),1) ; Ant 2 apol[*,2] = cfg.antpol(ishft(*obseq.ppol and 4,-2),1) ; Ant 4 endif ; On Mar 18, 2004, Ant 5 feed was changed to R and L, so change apol accordingly if (header.tls.yrday ge 2004.078) then apol[*,3] = cfg.antpol(ishft(*obseq.ppol and 4,-2),1) ; Ant 5 ; On Mar 22, 2004, Ant 6 feed was changed to R and L, so change apol accordingly if (header.tls.yrday ge 2004.082) then apol[*,4] = cfg.antpol(ishft(*obseq.ppol and 4,-2),1) ; Ant 6 ; On Aug 11, 2004, Ant 7 feed was changed to R and L, so change apol accordingly if (header.tls.yrday ge 2004.224) then apol[*,5] = cfg.antpol(ishft(*obseq.ppol and 4,-2),1) ; Ant 7 ; In special chanphz calibration MUX sequence 17, insert info for ; ants 1, 2 and 4 into ant locations 6, 7 and 8. if (cfg.muxseqid eq 17) then apol[*,4:6] = apol[*,0:2] for j = 0, nant-1 do begin for i = 0, j do begin k = header.lchoff(iants[i],iants[j]) if (i eq j) then begin pol(k,*) = apol(*,iants[i]) if (header.aatab[i] ge 4) then begin ; This is a 2m TP channel ant1 = where(header.aatab eq 1,n1) ant2 = where(header.aatab eq 2,n2) if (n1 gt 0) then begin lcp1 = where(pol[header.lchoff[iants[ant1],iants[ant1]],*] eq 1,nlcp) if (nlcp ne 0) then pol[k,lcp1] = 4 ; Change this 2m TP channel to poln 2 ; which signifies save as LCP and LIN endif if (n2 gt 0) then begin rcp2 = where(pol[header.lchoff[iants[ant2],iants[ant2]],*] eq 0,nrcp) if (nrcp ne 0) then pol[k,rcp2] = 2 ; Change this 2m TP channel to poln 4 ; which signifies save as RCP and LIN endif endif endif else begin ; Channel poln will always be lesser of two (only possibilities ; are 0-0, 0-3, 1-1, 1-3, or 3-3, corresp. to polns 0 0 1 1 3) pol(k,*) = apol(*,iants[i]) < apol(*,iants[j]) pol(k+1,*) = pol(k,*) endelse endfor endfor endif else begin ; Case of all antennas having dual polarization. This handles both XPOL=0 and XPOL=1 cases apol[*,0] = cfg.antpol[*obseq.ppol and 1,0] ; Ant 1 polarization apol[*,1] = cfg.antpol[ishft(*obseq.ppol and 2,-1),1] ; Ant 2 polarization for i = 2, 7 do apol[*,i] = cfg.antpol[ishft(*obseq.ppol and 4,-2),i] ; Small ant poln for j = 0, nant-1 do begin for i = 0, j do begin k = header.lchoff[iants[i],iants[j]] if (i eq j) then begin pol[k,*] = apol[*,iants[i]] endif else begin if (xpol) then begin pol[k,*] = (apol[*,iants[i]] xor apol[*,iants[j]]) - 1 + apol[*,iants[i]]*(1-apol[*,iants[j]]) pol[k+1,*] = pol[k,*] endif else begin pol[k,*] = (apol[*,iants[i]] and apol[*,iants[j]]) - (apol[*,iants[i]] xor apol[*,iants[j]]) pol[k+1,*] = pol[k,*] endelse endelse endfor endfor endelse msoff = get_timing(pol,header,cfg,obseq) ; Replace current polarization array with new one, first ; freeing the memory associated with the pointer ppol = ptr_new(pol) ; Parallel array to bigdata, containing poln code ; Gather information about indexes of a cycle at each frequency h = *obseq.pharm hord = *obseq.phord nf = n_elements(hord) ; Number of frequencies, regardless of polarization nrpt = intarr(nf) ; Array of repeated indexes in h, for each value in hord for i = 0, nf-1 do begin nrpt(i) = n_elements(where(h eq hord(i))) endfor nmax = max(nrpt) ; Greatest number of "repeats" at a single frequency indx = intarr(nmax,nf) ; Array to hold indexes for repeated frequencies indx = indx*0 - 1 ; Set all entries to -1 for i = 0, nf-1 do begin indx(0:nrpt(i)-1,i) = where(h eq hord(i)) endfor pnrpt = ptr_new(nrpt) pindx = ptr_new(indx) pmsoff = ptr_new(msoff) ; Return the pointers in a CYCLE structure cycle = {pbigdata:pbigdata, ppol:ppol, pnrpt:pnrpt, pindx:pindx, pmsoff:pmsoff} return,rec end