;+ ; NAME: ; SURVANAL ; PURPOSE: ; Reduce SURVEY observations and produce an overview plot ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; survanal,filename,hrec ; 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. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; ROUTINES CALLED: ; openarc, newscan, getcycle, flagool, gaincor ; OUTPUTS: ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 16-Jun-2001 by Dale E. Gary ;- pro survanal,filename,hrec ; 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 idx = ovsa_get_index(a,nrec,/create) rec = findseg(a,hrec,nrec,!SEGM.EPHEM,idx=idx) if (rec le 0) then begin ephem = get_eph_struct() ; Get an empty structure endif else begin data = getdata(rec-1,a) ephem = decode(data,!SEGM.EPHEM) endelse srcname = traj.srcname if (srcname ne 'SURVEY ') then begin print,'SURVANAL: Source name not consistent with SURVEY scan.' return endif nsrc = traj.nx ; Number of lines in Trajectory ; First 2 lines are not source lines, so exclude them xoff = (*traj.pha)[2:nsrc-1]*3.6 ; X pointing offset, in arcsec yoff = (*traj.pdec)[2:nsrc-1]*3.6 ; Y pointing offset, in arcsec cmds = (*traj.pcmd)[2:nsrc-1] dur = (*traj.pdur)[2:nsrc-1] nsrc = nsrc-2 ; Do some checks to ensure we are not making a mistake test = where(strtrim(cmds) eq 'TRACK',ntest) if (ntest ne nsrc) then begin print,'SURVANAL: Trajectory commands not consistent with SURVEY scan.' return endif test = where(dur eq 120,ntest) if (ntest ne nsrc) then begin print,'SURVANAL: Trajectory durations not consistent with SURVEY scan.' return endif if (xoff[nsrc-1] ne 0.0 or yoff[nsrc-1] ne 0.0) then begin print,'SURVANAL: Last pointing in trajectory is not Sun Center.' return endif ; Compare pointing offsets in TRAJECTORY with active region offsets ; as a means of identifying active regions. ; First convert all EPHEM solar coordinates to x,y nregn = ephem.regn.nregn ; mjd = cvdoy(ephem.regn.epoch.yr,ephem.regn.epoch.doy,/jul) rad = ephem.solar.rad[0]*0.36 ; Convert radius from 0.1 millideg to arcsec pa = ephem.solar.pa[0]/10000. ; Covert p-angle from 0.1 millideg to deg b0 = ephem.solar.b0[0]/10000. ; Covert b-angle from 0.1 millideg to deg ddif = 16./24 rx = fltarr(nregn) ry = fltarr(nregn) for i = 0,nregn-1 do begin lat = ephem.regn.info[i].lat lng = -ephem.regn.info[i].lng ; Rotate from 0 UT to 16 UT lng = solrotate(lng,lat,ddif) cnv_h2xy,pa,b0,rad,lat,lng,x,y rx[i] = x ry[i] = y endfor ; Now identify each active region by finding the closest regname = strarr(nsrc) for i = 0, nsrc-2 do begin dx = xoff[i]-rx dy = yoff[i]-ry dr = sqrt(dx^2+dy^2) rmin = min(dr,imin) print,'Min distance is ',rmin,' arcsec, for AR ',ephem.regn.info[imin].noa regname[i] = 'AR'+string(ephem.regn.info[imin].noa,format='(i4)') endfor regname[nsrc-1] = 'SUN-CENTER' ; 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) tcycle = 0L cycle = get_cycle(a,data,rec,header,obseq,cfg,tcycle) ; 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 onidx = where(*obseq.pnd eq 1,non) 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 nant = header.nant nchan = nant*nant ncycx = nsrc*60000./(cfg.sampintms*obseq.nx) ; Expected number of cycles out = fltarr(nchan,nf,npol,ncycx) id = progmeter(/init,button='Abort') 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) out[*,*,*,i] = smldata i = i + 1 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,endtime if (cycle eq !SEGM.HEADER) then goto,endtime goto,skiprec endif endwhile endtime: abort: free_lun,lun res = progmeter(id,/destroy) ncyc = i if (ncyc ne ncycx) then begin print,'SURVANAL: Warning, expected ',ncycx,' data points, but got ',ncyc out = out[*,*,*,0:ncyc-1] endif ; Average over samples for each region k1 = ceil(7.5*indgen(nsrc))>0 k2 = (ceil(7.5*(indgen(nsrc)+1))-1)