;+ ; NAME: ; TPCALCHEK2 ; PURPOSE: ; Analyze total power calibration data (called by TPCALCHEK main ; routine, which is just a date-sensitive shell). Responds ; to scans of type REFTPCAL (scan code 5). ; CATEGORY: ; OVRO APC DATA CALIBRATION ; CALLING SEQUENCE: ; tpcalchek2[,filename][,hrec | ,after=after][,debug][,/cmdfile] ; INPUTS: ; filename the name of the file containing TPCAL data. If ; omitted, the file DAILY.ARC in directory ; !DEFAULTS.WORKDIR is assumed. ; hrec the record number of the TPCAL scan header. Data ; will be processed up to the next EOS segment. ; If this argument is given, then AFTER keyword ; is ignored. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; after an optional time string of the standard form ; [yyyy.ddd ]hh:mm[:ss]m after which to start ; looking for a valid TPCAL scan. ; cmdfile a switch that indicates that the routine was called ; from a command file, so issue no modal messages, and ; print any error messages to \TPCAL.MSG. ; debug a switch that activates some debugging statements ; for getting an indication of intermediate steps. ; This switch is not compatible with /CMDFILE switch. ; ROUTINES CALLED: ; openarc, lasthrec, getdata, tl_decode, decode, tpnorm, gettpnd, ; tpsolve2, make_tpcal2, plot_tpcal, write_tpseg ; OUTPUTS: ; COMMENTS: ; This version works with 2-element data, where the TRAJECTORY is ; defined to take on-source data for two full cycles, then off-source ; data for two full cycles. ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 28-Aug-2000 by Dale E. Gary ; 25-Nov-2000 DG ; Changed call to PLOT_TPCAL, to reflect changes to that routine. ; 30-Nov-2000 DG ; Routine was not tolerant of extra non-data records in a scan. ; Changed loop reading BIGDATA to ignore non-data records. ; 18-Jul-2001 GN ; Replaced get_index by ovsa_get_index to avoid conflict with Yohkoh ;- pro tpcalchek2,filename,hrec,after=after,debug=debug,cmdfile=cmdfile ; If the /cmdfile switch is set, open a "new" TPCAL.MSG before doing ; anything. If it exists it will be cleared. msglun = 0 ; Indicates TPCAL.MSG does not exist IF (keyword_set(cmdfile)) THEN $ openw,msglun,/get_lun,!defaults.workdir+'TPCAL.MSG' ; If no filename, assume standard DAILY.ARC file IF (n_elements(filename) EQ 0) THEN $ filename=!defaults.workdir+'DAILY.ARC' ; Open the file lun = openarc(filename,a,nrec) IF (n_elements(hrec) EQ 0) THEN BEGIN IF (keyword_set(after)) THEN hrec = lasthrec(a,nrec,tstr=after) $ ELSE hrec = lasthrec(a,nrec) ENDIF ; Read header info, etc. rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm) tls = header.tls ; We need the actual data records of the current TPCAL segment idx = ovsa_get_index(a,nrec,/create) rectp = findseg(a,hrec,nrec,!SEGM.TPCAL,idx=idx) oldtp = intarr(1024,2) oldtp[*,0] = getdata(rectp-1,a) oldtp[*,1] = getdata(rectp,a) ; TPCAL segment is two records long ; Verify that the specified record has a REFTPCAL scan code IF (header.tls.scancode NE !SCAN.REFTPCAL) THEN BEGIN errstr = 'TPCALCHEK: Header found is of wrong type.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) return ENDIF nant = header.nant ; Frequency list f = (*obseq.phord)*200. ; The TPCAL should have obseq.nx*obseq.nt (700*10) entries, each with ; nchan (=header.nws-nndw=25) channels, so declare an appropriately ; sized array. nndw = header.ndbw + header.nfgw + 1 ; Number of non-data words nchan = header.nws-nndw ; Number of channels (25) nbl = (nchan-nant)/2 ; Number of baselines nsamp = obseq.nx*obseq.nt ; Number of samples (178*28) bigdata = intarr(nchan,nsamp) ; Big enough to hold all data ; The number of records to expect should be nrecs=nsamp/header.nsr ; (700*10/33 = 212). Note that there should be samples in the 213th ; record, but CPC doesn't send it. nrecs = nsamp/header.nsr ; rec = rec + 1 datoff = header.rdoff+nndw ndat = header.nws-nndw ; Read data from the file into a single large array, BIGDATA, of ; size (25,7000), which is just a long list of all 25 data channels ; Loop over the 212 records for i = 0,nrecs-1 do begin data = getdata(rec,a) ; Read a record tls = tl_decode(data) if (tls.segmentcode eq !SEGM.DATA) then begin ; Loop over the 33 samples in each record for j = 0, header.nsr-1 do begin soff = j*header.nws ; Additional offset for each sample bigdata(*,i*header.nsr+j) = data(datoff+soff:datoff+soff+ndat-1) endfor endif else i = i - 1 rec = rec + 1 ; Increment record number endfor free_lun,lun nrpt = obseq.nt ; Number of times entire cycle is repeated npolstate = 4 ; Number of polarization states ntp = 2 ; Number of useful total power channels ; Get index of active antennas IANTS (e.g. [1,2] => [0,1], [4,5,6] => [2,3,4] antlist = [1,2,4,5,6,7,8] ; Translation from antenna number to index aa = header.aatab(where(header.aatab ne 255)) iants = intarr(nant) for i = 0, nant-1 do begin iants[i] = where(antlist eq aa[i]) endfor offsets = (*gparm.poffsets) ; Array of offsets from GCAL, one for each data channel nlparm = (*gparm.pgnp).nlparm ; Array of nonlinearity parameters from GCAL, one for each ant ; Nant x nant matrix of channel offsets. choff = intarr(nant,nant) for i = 0, nant-1 do begin for j = 0, nant-1 do begin choff[i,j] = header.lchoff(iants[i],iants[j]) endfor endfor ; Form array representing the indexes of upper part of matrix upper = intarr(nbl) k = 0 for i = 0, nant-2 do begin for j = i+1, nant-1 do begin upper[k] = i*nant+j k = k + 1 endfor endfor choff(upper) = choff(upper)+1 ; Trick to put cos channel numbers into choff ; Subtract number of non-data words (nndw) from all channel offsets choff = choff-nndw ; Convert offsets from nant x nant array to linear array arranged in channel order off = fltarr(nant*nant) off(reform(choff,nant*nant)) = offsets ; Find the most current offsets and nonlinearity ; off = getoff(tls.yrday) ; gnp = getgnp(tls.yrday) ; Set defaults in case no offsets or nonlinearity could be found ; IF (off(0) eq -1) THEN BEGIN ; off = [10.,7.] ; errstr = 'TPCALCHEK: No channel offsets available--using defaults.' ; IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ; ELSE ans = widget_message(errstr,/ERROR) ; ENDIF ELSE off = off(0:ntp-1) ; IF ((size(gnp))(0) eq 0) THEN BEGIN ; nlparm = [0.,0.] ; errstr = 'TPCALCHEK: No nonlinearity parm available--using ZERO.' ; IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ; ELSE ans = widget_message(errstr,/ERROR) ; ENDIF ELSE nlparm = gnp.nlparm(0:ntp-1) ; List of indexes where noise diode is off (OFFIDX) and on (ONIDX) ; for various polarization combinations as obtained from observing ; sequence (OBSEQ structure) offrr = where(*obseq.pnd eq 2 and *obseq.ppol eq 0) offrl = where(*obseq.pnd eq 2 and *obseq.ppol eq 5) onrr = where(*obseq.pnd eq 1 and *obseq.ppol eq 0) onrl = where(*obseq.pnd eq 1 and *obseq.ppol eq 5) offlr = where(*obseq.pnd eq 2 and *obseq.ppol eq 10) offll = where(*obseq.pnd eq 2 and *obseq.ppol eq 15) onlr = where(*obseq.pnd eq 1 and *obseq.ppol eq 10) onll = where(*obseq.pnd eq 1 and *obseq.ppol eq 15) ; Declare storage for main TPCAL array, TPTEMP, a 5-dimensional array of ; size (2,4,2,85,10) where: ; first index 0..1 ND off/on ; second index 0..3 RR,RL,LL,LR (center feed/feed select) ; third index 0..1 TP Channel number ; fourth index 0..84 Frequency (all except H89) ; fifth index 0..9 Repeat of entire cycle tptemp = fltarr(2,npolstate,ntp,n_elements(offrr),nrpt) koff = 2 ; Number of samples before DOSEQ starts ; Apportion data in BIGDATA array into TPTEMP. Note that data that ; are expected to be bad are flagged with zeroes in the *OBSEQ.PND array. ; Thus, these entries do not appear in either the OFFIDX or ONIDX list ; of indexes, and therefore are excluded from the TPTEMP array. for k = 0, nrpt-1 do begin ; Loop over the 10 measurements tptemp(0,0,*,*,k) = bigdata(0:ntp-1,koff+offrr+k*obseq.nx) ;ND off RR state tptemp(0,1,*,*,k) = bigdata(0:ntp-1,koff+offrl+k*obseq.nx) ;ND off RL state tptemp(0,2,*,*,k) = bigdata(0:ntp-1,koff+offll+k*obseq.nx) ;ND off LL state tptemp(0,3,*,*,k) = bigdata(0:ntp-1,koff+offlr+k*obseq.nx) ;ND off LR state tptemp(1,0,*,*,k) = bigdata(0:ntp-1,koff+onrr +k*obseq.nx) ;ND on RR state tptemp(1,1,*,*,k) = bigdata(0:ntp-1,koff+onrl +k*obseq.nx) ;ND on RL state tptemp(1,2,*,*,k) = bigdata(0:ntp-1,koff+onll +k*obseq.nx) ;ND on LL state tptemp(1,3,*,*,k) = bigdata(0:ntp-1,koff+onlr +k*obseq.nx) ;ND on LR state endfor ; Save a copy of TPTEMP and set all saturated data ('7FFF'X) to zero ; so that it does not affect total. Note that the WHERE() function ; returns -1 if the condition is not true for any indices, so the ; IF statement is needed--otherwise it is "crash-city." This general ; technique is used often below. tpgood = tptemp saturated = where(tptemp eq '7FFF'X) ; Subtract channel offsets and apply nonlinearity parameters for ; each channel right up front for i = 0, ntp-1 do begin tpgood(*,*,i,*,*) = tpgood(*,*,i,*,*)-off(i) tpgood(*,*,i,*,*) = tpgood(*,*,i,*,*)*(1 + nlparm(i)*tpgood(*,*,i,*,*)/2047.) endfor ; Set the saturated values to zero IF (saturated(0) NE -1) THEN tpgood(saturated) = 0 ; Copy the array with good values into TPFLAG, then set all good values in ; TPFLAG to 1. TPFLAG then has only two values--0 where data were saturated ; and 1 where they are good. This is handy for indicating good vs bad ; data, but also for summing total numbers of good data. This general ; technque is used often below. tpflag = tpgood good = where(tptemp NE '7FFF'X) IF (good(0) EQ -1) THEN tpflag(*) = 1 $ ELSE tpflag(good) = 1 ; Divide the data into two parts--the first half is the on-source ; data, and the second half is the off-source data tpon = tpgood[*,*,*,*,0:nrpt/2-1] tponflg = tpflag[*,*,*,*,0:nrpt/2-1] tpoff = tpgood[*,*,*,*,nrpt/2:*] tpofflg = tpflag[*,*,*,*,nrpt/2:*] ; Work on ON-SOURCE data ;--------------------------------------------------------------------- ; Obtain information and statistics on the noise diode increment, from ; the data where the noise diode is on. IF (gettpnd(tpon,tponflg,ndfac,ndavg,ndrms) NE 0) THEN BEGIN errstr = 'TPCALCHEK: Error returned from routine GETTPND.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) ENDIF ; Now that we have the noise diode information, we can eliminate the ; data where the noise diode is on, and consider only the data where ; the noise diode is off tpgood = reform(tpon(0,*,*,*,*)) tpflag = reform(tponflg(0,*,*,*,*)) ; Take TPGOOD array, of size (4,2,85,NRPT), determine ON-OFF difference ; (RR-RL and LL-LR), and average over the NRPT independent measurements, ; returning the result of size (2,2,85) in TPAVG. The RMS deviation ; of the NRPT independent measurements is returned in TPRMS. IF (tpnorm(tpgood,tpflag,ndfac,ndavg,tpon,onrms) NE 0) THEN BEGIN errstr = 'TPCALCHEK: Error returned from routine TPNORM.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) ENDIF ; The noise diode increments for RR and LR (i.e. feed select R) are ; essentially identical. Likewise, those for LL and RL (feed select L) ; are also identical except for digitization noise of 1 or 2 DN. We ; can therefore average these. We do this in such a way that if one or ; the other of the two measurements is missing, the single good value is ; used. ; Form flag array ndflg = ndavg*0 good = where(finite(ndavg),ngood) IF (ngood gt 0) THEN ndflg(good) = 1 ;Zero the NaNs in NDAVG bad = where(finite(ndavg) eq 0,nbad) IF (nbad gt 0) THEN ndavg(bad) = 0 ; Now get sums ndavgr = reform(ndavg(0,*,*)+ndavg(3,*,*)) nptsr = reform(ndflg(0,*,*)+ndflg(3,*,*)) ndavgr = ndavgr/(nptsr>1) ndavgl = reform(ndavg(1,*,*)+ndavg(2,*,*)) nptsl = reform(ndflg(1,*,*)+ndflg(2,*,*)) ndavgl = ndavgl/(nptsl>1) ; Combine into a single array ndavg = ndavg(0:1,*,*) ndavg(0,*,*) = ndavgr ndavg(1,*,*) = ndavgl ; Set 0's back to NaN bad = where(ndavg eq 0,nbad) IF (nbad gt 0) THEN ndavg(bad) = !values.f_nan ; Do all of this for RMS also ; Form flag array ndflg = ndrms*0 good = where(finite(ndrms),ngood) IF (ngood gt 0) THEN ndflg(good) = 1 ;Zero the NaNs in NDRMS bad = where(finite(ndrms) eq 0,nbad) IF (nbad gt 0) THEN ndrms(bad) = 0 ; Now get sums ndrmsr = sqrt(reform(ndrms(0,*,*)^2+ndrms(3,*,*)^2)) nptsr = reform(ndflg(0,*,*)+ndflg(3,*,*)) ndrmsr = ndrmsr/(nptsr>1) ndrmsl = sqrt(reform(ndrms(1,*,*)^2+ndrms(2,*,*)^2)) nptsl = reform(ndflg(1,*,*)+ndflg(2,*,*)) ndrmsl = ndrmsl/(nptsl>1) ; Combine into a single array ndrms = ndrms(0:1,*,*) ndrms(0,*,*) = ndrmsr ndrms(1,*,*) = ndrmsl ; Set 0's back to NaN bad = where(ndrms eq 0,nbad) IF (nbad gt 0) THEN ndrms(bad) = !values.f_nan onnd = temporary(ndavg) onvar = temporary(ndrms) ; Now work on OFF-SOURCE ;--------------------------------------------------------------------- ; Obtain information and statistics on the noise diode increment, from ; the data where the noise diode is on. IF (gettpnd(tpoff,tpofflg,ndfac,ndavg,ndrms) NE 0) THEN BEGIN errstr = 'TPCALCHEK: Error returned from routine GETTPND.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) ENDIF ; Now that we have the noise diode information, we can eliminate the ; data where the noise diode is on, and consider only the data where ; the noise diode is off tpgood = reform(tpoff(0,*,*,*,*)) tpflag = reform(tpofflg(0,*,*,*,*)) ; Take TPGOOD array, of size (4,2,85,NRPT), determine ON-OFF difference ; (RR-RL and LL-LR), and average over the NRPT independent measurements, ; returning the result of size (2,2,85) in TPAVG. The RMS deviation ; of the NRPT independent measurements is returned in TPRMS. IF (tpnorm(tpgood,tpflag,ndfac,ndavg,tpoff,offrms) NE 0) THEN BEGIN errstr = 'TPCALCHEK: Error returned from routine TPNORM.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) ENDIF ; The noise diode increments for RR and LR (i.e. feed select R) are ; essentially identical. Likewise, those for LL and RL (feed select L) ; are also identical except for digitization noise of 1 or 2 DN. We ; can therefore average these. We do this in such a way that if one or ; the other of the two measurements is missing, the single good value is ; used. ; Form flag array ndflg = ndavg*0 good = where(finite(ndavg),ngood) IF (ngood gt 0) THEN ndflg(good) = 1 ;Zero the NaNs in NDAVG bad = where(finite(ndavg) eq 0,nbad) IF (nbad gt 0) THEN ndavg(bad) = 0 ; Now get sums ndavgr = reform(ndavg(0,*,*)+ndavg(3,*,*)) nptsr = reform(ndflg(0,*,*)+ndflg(3,*,*)) ndavgr = ndavgr/(nptsr>1) ndavgl = reform(ndavg(1,*,*)+ndavg(2,*,*)) nptsl = reform(ndflg(1,*,*)+ndflg(2,*,*)) ndavgl = ndavgl/(nptsl>1) ; Combine into a single array ndavg = ndavg(0:1,*,*) ndavg(0,*,*) = ndavgr ndavg(1,*,*) = ndavgl ; Set 0's back to NaN bad = where(ndavg eq 0,nbad) IF (nbad gt 0) THEN ndavg(bad) = !values.f_nan ; Do all of this for RMS also ; Form flag array ndflg = ndrms*0 good = where(finite(ndrms),ngood) IF (ngood gt 0) THEN ndflg(good) = 1 ;Zero the NaNs in NDRMS bad = where(finite(ndrms) eq 0,nbad) IF (nbad gt 0) THEN ndrms(bad) = 0 ; Now get sums ndrmsr = sqrt(reform(ndrms(0,*,*)^2+ndrms(3,*,*)^2)) nptsr = reform(ndflg(0,*,*)+ndflg(3,*,*)) ndrmsr = ndrmsr/(nptsr>1) ndrmsl = sqrt(reform(ndrms(1,*,*)^2+ndrms(2,*,*)^2)) nptsl = reform(ndflg(1,*,*)+ndflg(2,*,*)) ndrmsl = ndrmsl/(nptsl>1) ; Combine into a single array ndrms = ndrms(0:1,*,*) ndrms(0,*,*) = ndrmsr ndrms(1,*,*) = ndrmsl ; Set 0's back to NaN bad = where(ndrms eq 0,nbad) IF (nbad gt 0) THEN ndrms(bad) = !values.f_nan offnd = temporary(ndavg) offvar = temporary(ndrms) ; We now have the variables: ; TPON, ONRMS, ONND, ONVAR ; TPOFF, OFFRMS, OFFND, OFFVAR ; representing the total power and noise diode on-source and off-source. tpsolve2,tpon,onrms,onnd,onvar,tpoff,offrms,$ offnd,offvar,f,tls,oldtp,tpfac,tprms ; Create a TPCAL segment for inclusion into the raw data files. newtp = make_tpcal2(tls,tpfac,tprms,onnd,oldtp,nd_eff) ; Plot the result compared with the previous result, to check quality. plot_tpcal,newtp,oldtp ; Write the resulting TPCAL segment to the location saved within the ; REFTPCAL scan. result = write_tpseg(filename,hrec,!SEGM.TPCAL,newtp) IF (keyword_set(cmdfile)) THEN free_lun,msglun return END