;+ ; NAME: ; ANALYZE ; PURPOSE: ; Reduce observations by applying the phase and gain corrections ; and optionally performing a vector average over a specified ; time average and writing out an IDL save file of the results. ; Only one polarization at a time is currently handled. ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; analyze,filename,hrec,timavg,range[,time=time][,outfile=outfile] $ ; [,pshift=pshift][,/notrack][,xpol] ; INPUTS: ; filename the name of the file containing the data. ; hrec the record number of the scan header. Data ; will be processed up to the next EOS segment. ; If this argument is given, then AFTER keyword ; is ignored. ; timavg a single integer or float value: the time, in sec, over ; which to integrate. Must correspond to an integral number ; of cycles. Specify the cycle time for no integration. ; range a 2-element array indicating the range of data wanted, as ; either record numbers or times (the latter is ; indicated by the /TIME keyword). Times are specified ; as HHMMSS. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; time if set, specifies that the RANGE is to be interpreted as ; a time range. If unset or omitted, the RANGE is interpreted ; as a record range. ; outfile an optional string giving the name of the output SAVE file. If ; omitted, the data are placed in file !DEFAULTS.WORKDIR+GENERIC.SAV ; pshift an optional 2-element array specifying the distance to shift ; the phase center, in arcsec, as [x,y] with x positive W ; and y positive N. ; notrack an optional switch to inform FLAGOOL to disregard tracking ; (i.e. bad tracking data are NOT flagged bad). ; 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 in NEWSCAN. ; ROUTINES CALLED: ; openarc, newscan, getcycle, flagool, gaincor, freescan, ovsa_adjtime ; OUTPUTS: ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 07-Jan-1999 by Dale Gary ; 29-Jan-1999 DG ; Add outfile keyword ; 21-Jul-1999 DG ; Add output of TLS structure to .SAV file ; 28-Aug-1999 DG ; Add POLCODE, and change input time range from seconds to more ; convenient HHMMSS. ; 28-Mar-2000 DG ; Fix problem with routine getting into infinite loop when it hit ; the end of the file without an end of file record. ; 31-Mar-2000 DG ; Added saving of original datafile name as o_datfile ; 02-Apr-2000 DG ; Now works for all polarizations. Removed POLCODE keyword. ; 06-Apr-2000 DG ; Added saving of NDFAC, since this is needed for calibration. ; 14-Apr-2000 DG ; Get rid of C:\TEMP references (use !Defaults.workdir) ; 28-Jul-2000 DG ; Added HEADER to SAV file output, to provide NANT and AATAB ; Added code to extract uv values, and added some more outputs ; for Jeongwoo's mapping requirements. ANALYZE needs to be ; considerably revamped for "production" purposes soon. ; 04-Nov-2000 DG ; Added call to FREE_SCAN to free memory ; 06-Mar-2001 DG ; Added solar info to GEOMETRY structure (in NEWSCAN) and changed ; to write a single file instead of three. ; 18-Mar-2001 DG ; Previous change was a disaster for other code that assumed only ; one array name to hold the data, so all three polarizations, if ; they exist, are now contained in a single array. The variable ; NPOL gives info about which polarizations exist. NPOL=1 => I, ; NPOL=2 => R and L, NPOL=3 => R, L and I ; 15-Dec-2001 DG ; Changed TAVG from a FLOAT to a LONG array. Added OVSA_ADJTIME ; call, and save new variables MSOFF (time offsets for each data ; sample, in msec), TTP2M (time array for 2m double time resolution) ; and TP2M (total power for 2m at double time resolution). ; 08-Jun-2002 DG ; Added PSHIFT keyword to allow specifying a shift of the phase center. ; 17-Dec-2002 DG ; Added output of CFG structure to save file, and use /COMPRESS ; keyword. ; 23-Dec-2002 DG ; Added NOTRACK keyword ; 07-Nov-2003 DG ; Added XPOL keyword, to handle cross-polarization. ; 10-Jan-2004 DG ; Added error return in case of OPENARC failure ;- pro analyze,filename,hrec,timavg,range,time=time,outfile=outfile,$ pshift=pshift,notrack=notrack,xpol=xpol if (not keyword_set(outfile)) then outfile = !defaults.workdir+'generic.sav' if (not keyword_set(time)) then time = 0 else if (time ne 1) then time = 0 if (not keyword_set(pshift)) then pshift = [0.,0.] if (not keyword_set(xpol)) then xpol = 0 ; 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 (lun eq 0) then return ; Read header info, etc. rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str,pshift=pshift,xpol=xpol) - 1 ; Transfer needed part of cyc_str structure to more convenient variable names, then ; free the memory associated with those parts of cyc_str. bigdata = temporary(*cyc_str.pbigdata) pol = temporary(*cyc_str.ppol) nant = header.nant ; If the TIME switch was set, interpret the range as a time range if (time eq 1) then begin trange = hms2sec(range)*1000L if (range[1] lt range[0]) then trange[1] = trange[1]+hms2sec(240000L)*1000L erec = 0 endif else begin ; Otherwise interpret the range as a record range and set the ; time range to all day. trange = [0L,48*3600000L] ; Use 48 h in case there is a day change during scan if (range[0] gt hrec) then rec = range[0] erec = range(1) endelse ; Get cycles of data until start time is reached tcycle = 0L while (tcycle le trange(0)) do begin cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) ; tcycle = tcycle+3000 endwhile ; Want to average over TAVG seconds so determine number of cycles ; in TAVG s. If cycle time is not comensurate with TAVG s then ; averaging time can be slightly less, so that an integer number ; of cycles is averaged. cycletime = (obseq.nx*cfg.sampintms) ; Duration of cycle, in msec ; Check that timavg is defined sensibly, or set to 240. if not. if (n_elements(timavg) eq 0) then timavg = 240. $ else if (timavg lt cycletime/1000L) then timavg = cycletime/1000. ncyc2avg = fix(timavg*1000L/cycletime) ; Number of cycles to average ; Determine the number of good frequencies in the cycle nf = n_elements(*obseq.phord) offidx = where(*obseq.pnd eq 2) ; Indexes of obseq where noise diode is off nh = n_elements(*obseq.phord) ; Number of good harmonics f = (*obseq.phord)*0.2 ; Frequency list [GHz] ; Loop through all cycles in scan n = n_elements(cycle) i = 0 ; LEN is the number of cycles in the file if (erec gt 0) then len = (erec-hrec)*long(header.nsr)/obseq.nx $ else len = fix((trange(1)-trange(0))/cycletime)+1 navg = len/ncyc2avg nchan = nant^2 npol = 3 out = fltarr(nchan,nf,npol,ncyc2avg) t = lonarr(ncyc2avg) tavg = lonarr(navg) ; Declare storage for uv information uvout = fltarr(nant,nant,ncyc2avg) uvt = fltarr(nant,nant,navg) id = progmeter(/init,button='Abort') j = -1 avg = (sigma = fltarr(nchan,nf,navg,npol)) while(n_elements(cycle) eq n) do begin bigdata(0) = cycle ; Convert CYCLE to channelized floating array ; Apply flags for notrack, OOL if (keyword_set(notrack)) then flagool,bigdata,header,obseq,/notrack $ else flagool,bigdata,header,obseq ; Correct for complex gain and attenuation smldata = gaincor(bigdata,pol,header,cfg,obseq,tcycle,gparm,geometry,uv) ; Save the data into a large array out[*,*,*,i] = smldata uvout[*,*,i] = uv ; Save uv values, but no need for multi-frequency t(i) = tcycle if (ncyc2avg eq 1) then begin j = j + 1 res = progmeter(id,float(j)/navg) if (res eq 'Cancel') then goto,abort avg[*,*,j,*] = out tavg[j] = t[i] sigma(*,*,j,*) = !values.f_nan uvt[*,*,j] = uvout[*,*,i] i = -1 endif else begin if (i eq (ncyc2avg-1)) then begin ; We have read NCYC2AVG cycles, so time to average the data i = -1 j = j + 1 res = progmeter(id,float(j)/navg) if (res eq 'Cancel') then goto,abort ; Average the uv coordinates to represent the uv values for this ; averaged data sample. This may not be precisely correct in the ; case of missing/bad data, but since the missing data will depend ; on frequency, the correct way to do it would require the uvt ; array to depend on frequency as well. for iant = 0, nant-1 do begin for jant = 0, nant-1 do begin mom = moment(reform(uvout[jant,iant,*]),/nan) uvt[jant,iant,j] = mom[0] endfor endfor for k0 = 0, npol-1 do begin for k1 = 0, nchan-1 do begin for k2 = 0, nf-1 do begin ; Average the NCYC2AVG points for this channel and frequency ; Note that the MOMENT() function is screwed up--doesn't work ; when there is only one value, but works for none or any number ; other than 1. We have to write extra code to handle the ; one-value case good = where(finite(out(k1,k2,k0,*)),nfinite) if (nfinite eq 1) then begin avg[k1,k2,j,k0] = out[k1,k2,k0,good[0]] sigma[k1,k2,j,k0] = !values.f_nan endif else begin mom = moment(reform(out(k1,k2,k0,*)),/nan) avg[k1,k2,j,k0] = mom[0] sigma[k1,k2,j,k0] = mom[1] endelse tavg[j] = (moment(t))[0] endfor endfor endfor endif endelse if (j eq navg-1) then begin cycle = 0 endif else begin ; Read another cycle skiprec: cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) ; tcycle = tcycle+3000 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,endtime if (cycle eq !SEGM.HEADER) then goto,endtime if (cycle eq 0) then goto,endtime ; Case of end of data goto, skiprec endif if (tcycle gt trange(1)) then goto,endtime if (erec gt 0 and rec gt erec) then goto,endtime i = i + 1 endelse endwhile endtime: abort: free_lun,lun res = progmeter(id,/destroy) ; Get rid of any data at the end with zero times, since this signifies ; that our integer math was off slightly in calculating the number of ; points to expect. good = where(tavg ne 0,ngood) if (ngood gt 0) then begin avgg = temporary(avg[*,*,good,*]) sigm = temporary(sigma[*,*,good,*]) tavg = tavg[good] uvt = uvt[*,*,good] endif else begin print,'ANALYZE: No good data? No file written.' return endelse tls = header.tls ; Get calibrator source name, RA, and Dec from structures calinfo = {caleph} calinfo.name = traj.srcname calinfo.ra = long(geometry.a_ra[0]*15./360.) ; Convert msec to 0.1 mdeg calinfo.dec = long(geometry.a_dec[0]/360.) ; Convert masec to 0.1 mdeg ; Get the original datafile name for saving break_file,filename,disk,dir,file,ext o_datfile = file+ext ; Get Effective Noise Diode factors (these are to be divided into the ; data, in order to convert to ND_eff = 500). ndfac = (*gparm.pgnp).ndfac avgr = -1 avgl = -1 avgi = -1 sigmar = -1 sigmal = -1 sigmai = -1 ; Determine which polarizations exist, and eliminate missing ones from variables goodr = where(finite(avgg[*,*,*,0]),ngoodr) goodl = where(finite(avgg[*,*,*,1]),ngoodl) goodi = where(finite(avgg[*,*,*,2]),ngoodi) ; NPOL = 1 if only I poln exists, ; 2 if only R and L exist, and ; 3 if all three exist if (ngoodi eq 0) then begin ; If no I poln, assume only R and L exist avg = temporary(reform(avgg[*,*,*,0:1])) sigma = temporary(reform(sigm[*,*,*,0:1])) npol = 2 endif else if (ngoodr eq 0 and ngoodl eq 0) then begin ; If no R and no L, assume only I exists avg = temporary(reform(avgg[*,*,*,2])) sigma = temporary(reform(sigm[*,*,*,2])) npol = 1 endif else begin ; Must be all three polarizations avg = temporary(avgg) sigma = temporary(sigm) npol = 3 endelse msoff = temporary(*cyc_str.pmsoff) ovsa_adjtime,tavg,avg,msoff,ttp2m,tp2m,'SPLINE' if (n_elements(ttp2m) ne 0) then begin save,file=outfile,tavg,avg,sigma,npol,f,tls,ttp2m,tp2m,cfg,$ refcal,calinfo,o_datfile,ndfac,header,uvt,traj,geometry,msoff,/compress endif else begin save,file=outfile,tavg,avg,sigma,npol,f,tls,cfg,$ refcal,calinfo,o_datfile,ndfac,header,uvt,traj,geometry,msoff,/compress endelse print,'Successful Completion.' ; Free memory associated with structures freescan,obseq,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str return end