;+ ; NAME: ; HIGHRES ; PURPOSE: ; Create an "image" array of high resolution data (currently taken ; with 23 DOSEQ, but may be generalized later). Each row represents ; 4 s of data at 0.1 s per pixel, and the rows are stacked in time ; order. ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; highres,filename,hrec,range[,/time][,outfile] ; INPUTS: ; filename the name of the file containing the data. If ; omitted, the file DAILY.ARC in directory ; !DEFAULTS.WORKDIR is assumed. ; 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. ; range a two-element array containing the start and end ; range in records or if /TIME keyword is set, the ; time range in seconds. ; outfile optional output file name (default !DEFAULTS.WORKDIR+GENERIC.SAV) ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; time a switch. If 0, the RANGE is interpreted as a ; record range. If 1, the RANGE is interpreted as ; a time range in s. ; ROUTINES CALLED: ; openarc, getdata, lasthrec, tl_decode, decode, acfit ; OUTPUTS: ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; Initially written for 23 DOSEQ. May be generalized later. ; MODIFICATION HISTORY: ; Written 26-May-1999 by Dale Gary ; 05-Feb-2000 DG ; Removed now unnecessary sin and cos lookup tables in GAINCOR call. ; 14-Apr-2000 DG ; Get rid of C:\TEMP references (use !Defaults.workdir) ;- pro highres,filename,hrec,range,time=time,outfile=outfile 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 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) ; Read header info, etc. rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cyc_str) - 1 ; Check that this is indeed 23 DOSEQ data if (obseq.id ne 23) then begin ans = dialog_message('This scan was not taken with sequence 23 DOSEQ.',/error) return endif ; 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) cycletime = (obseq.nx*cfg.sampintms) ; Duration of cycle, in msec ; If the TIME switch was set, interpret the range as a time range if (time eq 1) then begin trange = range*1000L erec = 0 endif else begin ; Otherwise interpret the range as a record range and set the ; time range to all day. trange = [0L,24*3600000L] 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) endwhile ; 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] junk = where(*obseq.ppol eq 12,ni) ; NI contains number of I polarizations ; 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 = header.nant^2 h14 = where((*obseq.pharm)(offidx) eq 14,n14) out = fltarr(nchan,n14,len) ; t = lonarr(ncyc2avg) ; tavg = fltarr(navg) j = -1 ; avg = (sigma = fltarr(nchan,nf,navg)) while(n_elements(cycle) eq n) do begin 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,outdata) ; Save the H14 data j = j + 1 out(*,*,j) = outdata(*,h14) ; ; If NI is zero (no I polarizations) then save the R data for chans 5,6 ; if (ni eq 0) then out(5:6,*,i) = reform(smldata(5:6,*,0)) ; t(i) = tcycle ; if (ncyc2avg eq 1) then begin ; j = j + 1 ; print,j ; avg(*,*,j) = out ; tavg(j) = t(i) ; sigma(*,*,j) = !values.f_nan ; 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 ; print,j ; 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,*)),nfinite) ; if (nfinite eq 1) then begin ; avg(k1,k2,j) = out(k1,k2,good(0)) ; sigma(k1,k2,j) = !values.f_nan ; endif else begin ; mom = moment(reform(out(k1,k2,*)),/nan) ; avg(k1,k2,j) = mom(0) ; sigma(k1,k2,j) = mom(1) ; endelse ; tavg(j) = (moment(t))(0) ; endfor ; endfor ; endif ; endelse if (j eq len-1) then begin cycle = 0 endif else begin ; Read another cycle cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) 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: free_lun,lun save,file=outfile,out return end