;+ ; NAME: ; DLACHEK ; PURPOSE: ; Main routine to analyze delay center data. Responds to scans ; of type TEST DATA (scan code 1). ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; dlachek[,filename][,hrec | ,after=after][,debug][,/cmdfile] ; INPUTS: ; filename the name of the file containing PNTCAL data. If ; omitted, the file DAILY.ARC in directory ; !DEFAULTS.WORKDIR is assumed. ; hrec the record number of the SOLAR 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] after which to start ; looking for a valid PNTCAL or TEST DATA 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 \ACAL.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, getdata, lasthrec, tl_decode, decode, acfit ; OUTPUTS: ; COMMENTS: ; There will have to be modifications to this routine when ; additional antennas are available. ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 16-Jun-1998 by Dale Gary ; 21-Jun-1998 DG ; Changed from reading a(rec) to getdata(rec,a) throughout, ; so that routine should work under UNIX. ; 25-Aug-1998 DG ; Added check for new TRAJECTORY segments (but contents ignored ; for now.) ; 03-Nov-1998 DG ; Removed NaNs before using MEDIAN routine at the very end. ; 04-Feb-1999 DG ; Converted to use NEWSCAN routine ;- pro dlachek,filename,hrec,after=after,debug=debug,cmdfile=cmdfile ; 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) ; Verify that the specified record has either a DLASCAN or TEST DATA ; scan code IF (header.tls.scancode NE !SCAN.TEST AND header.tls.scancode NE !SCAN.DLASCAN) THEN BEGIN errstr = 'ACALCHEK: Header found is of wrong type.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/ERROR) return ENDIF ; The PNTCAL scan has NBLK blocks of data, each of duration NOBS seconds. ; These parameters are set by the TRAJECTORY. Currently this hardcodes ; the values for 1 TRAJECTORY, which are 24 blocks, 48 s per block. nblk = 24 nobs = 48 npersec = 1000/cfg.sampintms ; Number of samples per second koff = 2 ; Number of samples before DOSEQ starts ; We allow NAQC seconds to acquire the source, and allow NTOL seconds ; on either end for "sloppy timing." nacq = 12 ntol = 2 ; Length of observing sequence lseq = obseq.nx nant = header.nant nndw = header.ndbw + header.nfgw + 1 ; Number of non-data words nws = header.nws ; Number of words per sample nsr = header.nsr ; Number of samples per record nsamp = nblk*nobs*npersec ; Number of samples in entire traj bigdata = fltarr(nant,nsamp) ; Big enough to hold all of the data ; The number of records to expect should be nrecs=nsamp/header.nsr ; (24*48*10/33 = 349). Note that there should be some samples in the 350th ; record, but CPC doesn't send it. nrecs = nsamp/header.nsr ; rec = rec + 1 datoff = header.rdoff+nndw ndat = nant ; A late start may mean 1 fewer records, so we must be able to withstand ; a short scan. Use ON_IOERROR ON_IOERROR,EOD ; Read data from the file into a single large array, BIGDATA, of ; size (3,24*48*10), which is just a long list of 3 total power channels ; Loop over the 349 records for i = 0,nrecs-1 do begin data = getdata(rec,a) ; Read a record ; Loop over the 33 samples in each record for j = 0, header.nsr-1 do begin soff = j*header.nws+2 ; Additional offset for each sample ; +2 to skip over 27-m antennas bigdata(*,i*header.nsr+j) = data(datoff+soff:datoff+soff+ndat-1) endfor rec = rec + 1 ; Increment record number endfor EOD: ON_IOERROR,NULL ; Determine the start and stop indexes into the bigdata array where ; good data are expected. n = indgen(24) k1 = (n*nobs+nacq+ntol)*npersec - 1 +koff ; Array of end indexes of bad data k2 = [0,((n+1)*nobs-ntol)*npersec + 1]+koff ; Array of start indexes of bad data ; Set bad data parts to zero for i = 0, nblk-1 do begin bigdata(*,k2(i):k1(i)) = 0 endfor ; Create parallel array of flags indicating where the good data are bigflag = bigdata bigflag(where(bigdata ne 0)) = 1 k1 = (n*nobs+nacq+ntol)*npersec + koff ; Array of start indexes of good data k2 = ((n+1)*nobs -ntol)*npersec + koff ; Array of end indexes of good data ; Create array to hold the averaged data avdat = fltarr(nant,lseq/2,nblk) ; Loop over the 24 blocks (pointings) of data for i = 0, nblk-1 do begin ; Zero the accumulation arrays sum = fltarr(nant,lseq) nsum = intarr(nant,lseq) ; Determine start and end cycle numbers for this block. In general ; both the start and end cycles will be partially full, but where ; the data are outside of the good data interval they have already ; been set to zero above, so they will not contribute to the sums. istart = floor((k1(i)-koff)/float(lseq)) iend = ceil((k2(i)-koff)/float(lseq)) ; Accumulate data for this block for j = istart,iend-1 do begin ; Determine start and end indexes of first half of cycle j1 = j*lseq + koff j2 = (j+1)*lseq + koff - 1 - lseq/2 ; Accumulate first half of cycle sum = sum + bigdata(*,j1:j2) nsum = nsum + bigflag(*,j1:j2) ; Accumulate second half, which is just a repeat of the ; first half since R and L entries are both lin polarization sum = sum + bigdata(*,j1+lseq/2:j2+lseq/2) nsum = nsum + bigflag(*,j1+lseq/2:j2+lseq/2) endfor ; Determine the average for this block and store it avdat(*,*,i) = sum/nsum endfor ; AVDAT now contains the desired averages, but includes blocks ; 2, 4, 21, and 23, which are not wanted (these were blocks when ; there were extra long slew times, so should be discarded). AVDAT ; also contains entries in the observing sequence that are bad due ; to switching oscillators or rotating feeds. We must compress ; AVDAT to contain only the good data. ; First determine good entries in one half of observing sequence gdseq = where((*obseq.pharm)(0:lseq/2-1) ne -1) ; Make list of good blocks gdblks = [0,2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,21,23] ; Select only the good values--the statement below does it all! avdat = avdat(*,gdseq,gdblks) ; Perform fit to data contained in AVDAT, returning the peak level ; PK(nant,nf), the pointing offsets in degrees PO(nant,nf,2), ; the half-power-beamwidth HPBW(nant,nf,2), and rms error of the ; fit RMSERR(nant,nf,2) acfit,avdat,pk,po,hpbw,rmserr ; Get best offset values print,"Pointing offsets for this location are:" hoff = (doff = fltarr(nant)) for i = 0, nant-1 do begin poff = reform(po(i,*,0)) ; This code needed to eliminate poff = poff(where(finite(poff))) ; NaNs (median routine dies) hoff(i) = median(poff) poff = reform(po(i,*,1)) poff = poff(where(finite(poff))) doff(i) = median(poff) ant = header.aatab(i+2) print,ant,hoff(i),ant,doff(i),format='(i1,"HO = ",f7.3," ",i1,"DO = ",f7.3)' endfor free_lun,lun return END