;+ ; NAME: ; ACALCHEK ; PURPOSE: ; Main routine to analyze antenna pointing data. Responds to scans ; of type PNTCAL (scan code 8) or TEST DATA (scan code 1). ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; acalchek[,filename][,hrec | ,after=after][,output][,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: ; output a 2*nant-element string array containing formatted lines of ; output of the following form. The two header lines are ; not included. These lines are also appended to ; the file boresite.pnt in the directory !DEFAULTS.DBDIR. ; ; Source Name Date Time Rec# HA(deg) DEC(deg) Ant Axis Offset ; ------------ ---------- -------- ----- -------- -------- --- ---- ------ ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 4 HO 0.075 ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 4 DO -0.017 ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 5 HO 0.027 ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 5 DO 0.118 ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 6 HO 0.020 ; SUN 2000/03/08 19:09:40 1402 -13.531 -4.540 6 DO 0.037 ; ; 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 ; 07-Jun-1999 DG ; Began work on plotting diagnostic info, for evaluating the results ; 13-Mar-2000 DG ; Changed to work with more than 5 antennas. ; 04-Nov-2000 DG ; Major changes to make it work with VARDB, etc. Now applies ; proper gain calibration before determining pointing. ; 08-Dec-2003 DG ; Fixed bug that caused a crash when one antenna was missing (due to ; plotting with too few points). ; 01-Sep-2005 DG ; Several changes to plot frequency-dependent offsets (likely will be ; commented out except for special tests). ;- pro acalchek,filename,hrec,output,after=after,debug=debug,cmdfile=cmdfile ; If the /cmdfile switch is set, open a "new" ACAL.MSG before doing ; anything. If it exists it will be cleared. msglun = 0 ; Indicates ACAL.MSG does not exist IF (keyword_set(cmdfile)) THEN $ openw,msglun,/get_lun,!defaults.workdir+'acal.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,geometry,pRecent,pOld,cyc_str) - 1 ; Verify that data were taken with the correct observing sequence. If not, ; issue a warning (since other observing sequence numbers could still work). if (obseq.id ne 21) then begin errstr = 'ACALCHEK: Incorrect observing sequence '+string(obseq.id,format='(I3)')+$ '. Sequence 21 expected.' IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $ ELSE ans = widget_message(errstr,/INFO) endif ; Verify that the specified record has either a PNTCAL or TEST DATA ; scan code IF (header.tls.scancode NE !SCAN.TEST AND header.tls.scancode NE !SCAN.PNTCAL) 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 n27m = 2 nndw = header.ndbw + header.nfgw + 1 + n27m ; Number of non-data words (excludes 27 m) nant = header.nant - n27m ; Number of antennas excludes 27 m ; 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 ; Number of cycles expected ncyc = nblk*nobs*1000L/(cfg.sampintms*obseq.nx) ; Read all of the data into the array ACALDATA ; Explicitly copy "CYCLE structure" arrays bigdata = temporary(*cyc_str.pbigdata) pol = temporary(*cyc_str.ppol) ; Get a cycle of data cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) acaldata = fltarr(nant,obseq.nx,ncyc+1) n = n_elements(cycle) i = 0 while(n_elements(cycle) eq n) do begin ; Get a cycle of data bigdata(0) = cycle ; Convert CYCLE to channelized floating array flagool,bigdata,header,obseq ; Apply flags for notrack, OOL ; Correct for complex gain and attenuation smldata = gaincor(bigdata,pol,header,cfg,obseq,tcycle,gparm,geometry,uv) acaldata[*,*,i] = bigdata[nndw:nndw+nant-1,*] ; keep only 2 m channels i = i + 1 ; Get a new cycle of data skiprec: cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) if (n_elements(cycle) eq 1) then begin ; Check if this is a header record, meaning we have reached ; the end of the scan, otherwise, continue on if (cycle eq !SEGM.EOS) then goto,endloop if (cycle eq !SEGM.HEADER) then goto,endloop if (cycle eq 0) then goto,endloop ; Case of end of data goto, skiprec endif endwhile endloop: 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 frq = (*obseq.phord)*0.2 nant = header.nant - 2 ; Number of antennas excludes 27 m 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 = ceil(nblk*nobs*npersec+koff) ; Number of samples in entire traj bigdata = reform(acaldata,nant,obseq.nx*(ncyc+1)) ; Big enough to hold all of the data ; Determine the start and stop indexes into the bigdata array where ; good data are expected. n = indgen(nblk) k1 = nint((n*nobs+nacq+ntol)*npersec) - 1 +koff ; Array of end indexes of bad data k2 = [0,nint(((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 = nint((n*nobs+nacq+ntol)*npersec) + koff ; Array of start indexes of good data k2 = nint(((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)) if (iend*lseq gt nsamp) then iend = iend-1 ; 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;,/debug ; Get best offset values, integrated over frequency ; print,"Pointing offsets for this location are:" hoff = (doff = fltarr(nant)) ant = intarr(nant) fbad = [4,19,21] po[*,fbad,*] = !values.f_nan ; Eliminate 2.0 and 6.0 GHz (interference) for i = 0, nant-1 do begin poff = reform(po(i,*,0)) good = where(finite(poff),ngood) if (ngood ne 0) then begin ; This code needed to eliminate poff = poff(good) ; NaNs (median routine dies) hoff(i) = median(poff) endif else hoff(i) = -10000. poff = reform(po(i,*,1)) good = where(finite(poff),ngood) if (ngood ne 0) then begin ; This code needed to eliminate poff = poff(good) ; NaNs (median routine dies) doff(i) = median(poff) endif else doff(i) = -10000. ant[i] = header.aatab(i+2) ; print,ant[i],hoff(i),ant,doff(i),format='(i1,"HO = ",f7.3," ",i1,"DO = ",f7.3)' endfor ; Issue the following line (outside of this routine) before using the commented part ; below... ; window,20 & !p.multi=[0,5,2] & for i = 0, 9 do plot,/nodata,[1,18],[-20,20],/xsty,/ysty ; Uncomment the next 10 lines for diagnostics of frequency-dependent feed pointing ;wset,20 ;k = 2*nant ;for i = 0, nant-1 do begin ; for j = 0, 1 do begin ; !p.multi[0]=k ; plot,[1,18],[-20,20],/nodata,/xsty,/ysty,/noerase ; oplot,frq,(po[i,*,j]-([hoff[i],doff[i]])[j])*60. ; k = k - 1 ; endfor ;endfor ; Duration of measurement is 19 min, 12 sec, so offset time in the ; header by 1/2 of this in [msec], to get time of center of observation dtmsec = ((19*60. + 12)*1000.)/2 tmsec = header.tls.msec+dtmsec ; Get Right Ascension and Declination at this time tarray = [1D,tmsec,tmsec^2] ra = total(geometry.a_ra*tarray) ; [msec] dec = total(geometry.a_dec*tarray) ; [masec] ; Calculate HA ha = lst(header.tls.year,header.tls.day,tmsec) - ra ; [msec] ha = ha*15./3600.D3 if (ha lt -180.) then ha = 360.+ha else if (ha gt 180.) then ha = ha-360. dec = dec/3600.D3 nf = n_elements(frq) msav = !p.multi csav = !p.charsize window,2,xsiz=220*nant,ysiz=550 device,decompose=0 loadct,39 !p.multi=[0,nant,3,0,0] !p.charsize=1.5 for i = 0, nant-1 do begin good = where(finite(pk[i,*]),ngood) if (ngood gt 1) then begin plot,frq,pk(i,*),psym=-1,xtit='Frequency [GHz]',ytit='Amplitude [raw units]',symsize=0.5 oplot,frq(0:20),pk(i,0:20),psym=-1,color=250,symsize=0.5 oplot,frq(20:35),pk(i,20:35),psym=-1,color=180,symsize=0.5 oplot,frq(35:nf-1),pk(i,35:nf-1),psym=-1,color=60,symsize=0.5 oplot,frq,pk(i,*)*sqrt(rmserr(i,*,0)^2+rmserr(i,*,1)^2)/2,psym=3 endif else begin plot,[0,1],[0,1],psym=-1,xtit='Frequency [GHz]',ytit='Amplitude [raw units]',symsize=0.5 xyouts,0.5,0.5,'No data',align=0.5 endelse endfor ; window,3,xsiz=650,ysiz=180 for i = 0,nant-1 do begin plot,[-0.5,0.5],[0,0],xsty=1,ysty=1,xtit='HA Offset [deg]',ytit='Dec Offset [deg]',$ xran=[-0.8,0.8],yran=[-0.8,0.8] oplot,[0,0],[-0.5,0.5] oplot,po(i,0:20,0),po(i,0:20,1),psym=1,color=250,symsize=0.5 oplot,po(i,20:35,0),po(i,20:35,1),psym=1,color=180,symsize=0.5 oplot,po(i,35:nf-1,0),po(i,35:nf-1,1),psym=1,color=60,symsize=0.5 plots,hoff(i),doff(i),psym=1,color=255,symsize=3 endfor bsiz = (46.5/60./frq)*13.5 for i = 0,nant-1 do begin plot_io,frq,hpbw(i,*,0),psym=-1,xran=[0,20],yran=[0.3,10],$ xtit='Frequency [GHz]',ytit='Beam Width [degrees]',xsty=1,ysty=1,symsize=0.5 oplot,frq(0:20),hpbw(i,0:20,0),psym=-1,color=250,symsize=0.5 oplot,frq(20:35),hpbw(i,20:35,0),psym=-1,color=180,symsize=0.5 oplot,frq(35:NF-1),hpbw(i,35:nf-1,0),psym=-1,color=60,symsize=0.5 oplot,frq(0:20),hpbw(i,0:20,1),psym=-6,color=250,symsize=0.5 oplot,frq(20:35),hpbw(i,20:35,1),psym=-6,color=180,symsize=0.5 oplot,frq(35:NF-1),hpbw(i,35:nf-1,1),psym=-6,color=60,symsize=0.5 oplot,frq, bsiz endfor !p.multi=msav !p.charsize=csav olun = 0 openu,olun,/get_lun,!defaults.dbdir+'boresite.pnt',/append ; Write out the header lines if this is a new file (CUR_PTR = 0) if ((FSTAT(olun)).CUR_PTR eq 0) then begin out = 'Source Name Date Time Rec# HA(deg) DEC(deg) Ant Axis Offset' printf,olun,out,format='(a)' out = '------------ ---------- -------- ----- -------- -------- --- ---- ------' printf,olun,out,format='(a)' endif on_ioerror,NULL line = string(traj.srcname,format='(A12)')+' '+header.tls.date ; Get best offset values print,traj.srcname ; Set output to a dummy (blank) string output = ' ' ; Write out the lines for the NANT 2m antennas for i = 0, nant-1 do begin out = line+strmid(msec2str(tmsec),0,9)+string(hrec,format='(I6)')$ +string(ha,dec,format='(2F9.3)')+string(ant[i],format='(I3)') out2 = out+' HO '+string(hoff[i],format='(F6.3)') printf,olun,out2 output = [output,out2] out2 = out+' DO '+string(doff[i],format='(F6.3)') printf,olun,out2 output = [output,out2] endfor ; Remove the initial blank string from output output = output[1:2*nant] free_lun,olun ; Free memory associated with structures freescan,obseq,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str free_lun,lun return END