;+ ; NAME: ; CHECK_12AP ; PURPOSE: ; APCAL_ADD_27M_RL support routine. Checks 12 baseline data for ; multiple .AMF files, plots some comparison data, and returns ; appropriately averaged amplitudes and phases. ; CATEGORY: ; OVRO APC DATA CALIBRATION ; CALLING SEQUENCE: ; check_12ap,f,amp,damp,phz,dphz ; INPUTS: ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; ROUTINES CALLED: ; OUTPUTS: ; f Array of size [nf] containing the frequencies used. NB: must be the ; same for all .AMF files (61 DOSEQ) ; amp Array of size [nf,npol] containing amplitudes for 12 baseline ; damp Array of size [nf,npol] containing amplitude rms deviation for 12 baseline ; phz Array of size [nf,npol] containing phases for 12 baseline ; dphz Array of size [nf,npol] containing phase rms deviation for 12 baseline ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 26-Jul-2005 by Dale E. Gary ; 05-Sep-2005 DG ; Used LOBE function to fix problem in which output phases (PHZ array) were ; ranging from -180 to 180. ; 20-Jan-2014 DG ; Added filelist keyword so that this can be automated. ;- pro check_12ap,f,amp,damp,phz,dphz,filelist=filelist default,filelist,'' if (filelist eq '') then begin ; Request list of AMF files from user filelist = dialog_pickfile(path=!defaults.datadir,filter='*.amf',/multi,/read,$ title='Select all .AMF files you wish to include.') endif if (filelist[0] eq '') then return filelist = filelist[sort(filelist)] nfile = n_elements(filelist) ; Read the first AMF file to get sizes, and create some arrays for ; processing the AMF data amf = readamf(filelist[0],info=info) f = reform(amf[0,*,0,0]) nf = n_elements(f) npol = n_elements(amf[0,0,0,*]) if (npol ne 2) then begin print,'CHECK_12AP: Error, NPOL=',npol,' but this routine only works for R,L data.' return endif p12 = fltarr(nf,npol,nfile) dp12 = fltarr(nf,npol,nfile) a12 = fltarr(nf,npol,nfile) da12 = fltarr(nf,npol,nfile) bc = bytarr(nfile) ab = bytarr(nfile) ; Initialize the phase plot windows window,0 plot,[1,18],[-180,180],/nodata,title='Phase Difference R-L',/ysty,/xsty window,1 plot,[1,18],[-180,180],/nodata,title='Phases Relative to first file',/ysty,/xsty psyms = [1,2,4,5,6] ; Loop over AMF files for i = 0, nfile-1 do begin ; Read the AMF file amf = readamf(filelist[i],info=info) ; Convert bad values to NaN bad = where(amf eq -99.9,nbad) if (nbad ne 0) then amf[bad] = !values.f_nan bc[i] = info[0].BE eq 'BC' ; 1 where baseline is BC, 0 elsewhere ab[i] = info[0].BE eq 'AB' ; 1 where baseline is AB, 0 elsewhere ; Put data into working arrays, correcting the phases for CHANPHZ and ORDER a12[*,*,i] = reform(amf[1,*,0,*]) da12[*,*,i] = reform(amf[2,*,0,*]) p12[*,*,i] = lobe((reform(amf[3,*,0,*])-info[0].chanphz)*info[0].order) dp12[*,*,i] = reform(amf[4,*,0,*]) pdif = lobe(p12[*,0,i] - p12[*,1,i],/mid) ; Window 0 contains R-L phases, as a sanity check. They should agree ; for all measurementes wset,0 oplot,f,pdif,psym=psyms[i] if (i eq 0) then begin ; Use the first file in the list as reference phases prefr = p12[*,0,i] prefrerr = amf[4,*,0,0] prefl = p12[*,1,i] preflerr = amf[4,*,1,0] endif else begin ; Window 1 contains 12 phase of current file ; minus 12 phase of reference (first) file, overplotted ; with best fit straight line wset,1 ; Do RCP p = lobe(p12[*,0,i]-prefr,/mid) perr = amf[4,*,0,0] meas = sqrt(perr^2+prefrerr^2)/2. oplot,f,p,psym=psyms[i] good = where(finite(p) and finite(meas),ngood) if (ngood ne 0) then begin res = linfit(f[good],p[good],meas=meas[good],chisq=chisq) oplot,[1,18],res[0]+res[1]*[1,18] ; Correct phase for slope and offset p12[*,0,i] = lobe(p12[*,0,i]-res[0]-res[1]*f) endif ; Repeat for LCP, but without plotting p = lobe(p12[*,1,i]-prefl,/mid) perr = amf[4,*,1,0] meas = sqrt(perr^2+preflerr^2)/2. good = where(finite(p) and finite(meas),ngood) if (ngood ne 0) then begin res = linfit(f[good],p[good],meas=meas[good],chisq=chisq) ; Correct phase for slope and offset p12[*,1,i] = lobe(p12[*,1,i]-res[0]-res[1]*f) endif endelse endfor ; We have the corrected phases for each file, so now generate a single median phase phz = fltarr(nf,npol) dphz = fltarr(nf,npol) if (nfile eq 1) then begin phz = reform(p12[*,*,0]) dphz = reform(dp12[*,*,0]) endif else begin for i = 0, nf-1 do begin for j = 0, npol-1 do begin ; Take median of phase at each frequency ; Convert to x,y (using unit amplitude), find median, then convert back to phase x = cos(p12[i,j,*]*!dtor) y = sin(p12[i,j,*]*!dtor) xmed = median(x,/even) ymed = median(y,/even) phz[i,j] = lobe(atan(ymed,xmed)/!dtor) ; Take internal fluctuation as dphz ; Convert to x,y, find fluctuation in each axis, convert to total fluctuation ; then convert back to phase (using unit amplitude) xres = moment(x-xmed,/nan) yres = moment(y-ymed,/nan) sigx = sqrt(xres[1]) sigy = sqrt(yres[1]) sig = sqrt((sigx*sin(phz[i,j]*!dtor))^2 + (sigy*cos(phz[i,j]*!dtor))^2) dphz[i,j] = atan(sig,1)/!dtor endfor endfor endelse ; Now handle amplitudes window,2 plot,[1,18],[0,1.5],/nodata,title='Amplitudes',/ysty,/xsty window,3 plot,[1,18],[0,1.5],/nodata,title='Amplitudes',/ysty,/xsty ; I discovered that the GCAL NDEFF only applies to 12 on the BC channels, ; since that is the MUX mode when it is measured. However, 125,126 and 128 ; reference calibration data use 12 on the AB channels, and there is a ; difference in the NDEFF. The code below determines the scale factor ; between BC and AB baselines, and applies it to give consistent results ; for which the median can be determined. ibc = where(bc,nbc) iab = where(ab,nab) fac = replicate(1,nfile) if (nbc ne 0) then begin amp_bc = total(reform(a12[*,0,ibc]),2)/nbc ; Average of BC amplitudes if (nab ne 0) then begin ab_fac = (total(reform(a12[*,0,iab]),2)/nab)/amp_bc ; Ratio of AB to BC amps. good = where(finite(ab_fac),ngood) if (ngood ne 0) then ab_fac = total(ab_fac,/nan)/ngood ; FAC is an array of length NFILE, with AB factors for AB baselines ; and unity for BC baselines, e.g. [1.0,0.68,0.68,1.0,0.68] fac = ab_fac*ab + bc endif endif ; Loop over files, applying the factors for i = 0, nfile-1 do begin wset,2 a12[*,0,i] = a12[*,0,i]/fac[i] oplot,f,a12[*,0,i],psym=psyms[i] wset,3 a12[*,1,i] = a12[*,1,i]/fac[i] oplot,f,a12[*,1,i],psym=psyms[i] endfor ; The amplitudes are now consistent, so determine median amp = fltarr(nf,npol) damp = fltarr(nf,npol) for i = 0, nf-1 do begin for j = 0, npol-1 do begin amp[i,j] = median(a12[i,j,*],/even) res = moment(a12[i,j,*],/nan) ; For DAMP, use greater of internal error and formal error in AMF files damp[i,j] = max([min(da12[i,j,*],/nan),sqrt(res[1])]) endfor endfor ; Overplot the amplitudes wset,2 oploterr,f,amp[*,0],damp[*,0] wset,3 oploterr,f,amp[*,1],damp[*,1] return end