pro analyze,filename,hrec,timavg,range,time=time,outfile=outfile,pshift=pshift,notrack=notrack,auto=auto,error=error
;This is the code to handle the errors
error=''
CATCH, Error_status
;This statement begins the error handler.
IF Error_status NE 0 THEN BEGIN
error='ANALYZE:'+!ERR_STRING
PRINT, 'Error:', error
CATCH,/Cancel
if n_elements(lun) ne 0 then if lun ne 0 then free_lun,lun
return
ENDIF
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 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,auto=auto)
if lun eq 0 then message, 'OPENARC:cannot open this file'
; Read header info, etc.
rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str,pshift=pshift) - 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 = n_elements(uniq(*obseq.ppol,sort(*obseq.ppol)))
if (header.tls.yrday lt 2004.231) then npol = 3
if (header.tls.yrday ge 2004.231) then npol = 2
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 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
msoff = temporary(*cyc_str.pmsoff)
if (header.tls.yrday lt 2004.231) then begin
; 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
ovsa_adjtime,tavg,avg,msoff,ttp2m,tp2m,'SPLINE',/old
endif else begin
avg = temporary(avgg)
sigma = temporary(sigm)
ovsa_adjtime,tavg,avg,msoff,ttp2m,tp2m,'SPLINE'
endelse
;If outfile is not specified define the standard filename for for the outputfile
;based on date and start time
if n_elements(outfile) eq 0 then begin
junk = str_sep(header.tls.date,'/',/remove_all)
outfile = !defaults.workdir+junk[0]+junk[1]+junk[2]+$
'_'+strmid(header.tls.timstr,1,2)+strmid(header.tls.timstr,4,2)+'.sav'
end
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
en