;+ ; NAME: ; AMPHIT ; PURPOSE: ; Routine to fit best straight line to calibrator phases and ; return the average amplitude and phase for the entire timerange ; of data, with optional plotting and saving of data to a file. ; This is the main routine for creating amplitude and phase ; calibration. The handles only one baseline at a time and must ; be called multiple times to do all baselines. ; CATEGORY: ; OVRO APC DATA CALIBRATION ; CALLING SEQUENCE: ; amphit,frq,data,bl,tavg,cfg,out,[,pslope][,ifrq=ifrq] $ ; [,fitslope=fitslope][,debug=debug][,/save,lun=lun] $ ; [,pol=pol][,ndfac=ndfac] ; INPUTS: ; frq Array of frequencies contained in the data [GHz] ; data Data array of size (NCHAN,NFRQ,NMEAS), where NCHAN is the ; number of channels, NFRQ is the number of frequencies, ; and NMEAS is the number of measurements (times). ; bl A number representing the baseline to be fit. Currently one ; of [12,14,15,16,24,25,26]. ; tavg The array of times, of length (NMEAS), in [msec] corresponding ; to the data ; cfg The CONFIG structure from the scan that was analyzed. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; pslope A guess at the frequency-dependent slope that will remove ; temporal drifts in the phases. If omitted, zero is assumed. ; ifrq Specifies an index within the FRQ array, for handling a single ; frequency at a time. ; fitslope A switch declaring whether a frequency-dependent phase drift ; should be allowed in the fit. If set, a phase drift is ; allowed (and reported in the output). If not set, zero ; phase drift is assumed (after possible correction for ; the value of PSLOPE. ; debug A switch that turns plotting on or off. If set, plots of the ; phases as a function of time at each frequency are shown, ; along with the straight line fit. ; lun The logical unit number of the file opened by the routine on ; an initial call. This should be undefined on the first call, ; and unchanged on subsequent calls. ; pol A single string character specifying the polarization state ; 'R' -> RCP, 'L' -> LCP, 'I' -> Linear polarization. If ; omitted, 'I' is the default. Note that this information is ; merely written to the header line of the output. ; ndfac An array of size (NANT), giving the factors for each antenna ; that are needed to normalize the data to an effective noise ; diode value of 500. If omitted, the values are set to unity. ; ROUTINES CALLED: ; OUTPUTS: ; lun The logical unit number of the file opened by the routine on ; the initial call. ; out The result of the fit, in 6 columns for each frequency. The ; columns are Freq, Amp, dAmp, Pha, dPha, and Pslope in ; degrees/hr/GHz. ; COMMENTS: ; SIDE EFFECTS: ; Writes the solution to the output log or terminal. If the SAVE keyword ; is set, the solution is also written to a file (whose name is requested ; via widget dialog). When all of the baselines have been written, the ; file must be explicitly closed by the caller with FREE_LUN,LUN. ; RESTRICTIONS: ; When all of the baselines have been written, the file must be ; explicitly closed by the caller with FREE_LUN,LUN. The baselines ; should be written in the order 12, 14, 15, 16, 24, 25, 26 in order to ; make a valid calibration file. ; MODIFICATION HISTORY: ; Written 04-Jul-1999 by Dale E. Gary ; 01-Apr-2000 DG ; Added IFRQ keyword, for handling a single frequency at a time. Also ; now writes out frequency in OUT array. ; 06-Apr-2000 DG ; Added POL keyword, whose sole purpose is to allow writing of polarization ; state to the header lines of the output. Also added NDFAC keyword, to ; normalize the output to correspond to an effective noise diode value of ; 500. ; 15-Jun-2000 DG ; Fix problem with fitting non-5-element data. ; 02-Jul-2000 DG ; Fix a problem with 3-element data. ; 22-Jul-2000 DG ; Change averaging scheme to do algebraic average of amplitudes if S/N ; is greater than 2.5 ; 15-Oct-2002 DG ; No longer sets frequency column to NaN when there are no good data. ; 17-Dec-2002 DG ; Made a number of changes to reflect new calibration scheme. Now ; both the HEADER and CONFIG structures are included in the arg list, ; and information is written to the output file giving the antennas, ; their order, which backend channels were used, and what the relevent ; backend channel phase offset is. ; 03-Jan-2003 DG ; I had to resort to a "trick" to deal with a problem with the 12RLI ; data. If the logical unit number is not defined, some header lines ; are written to the beginning of the file. The first line contains ; information assuming this is 12RLI data. This is okay, since WAMPHIT, ; which is used for non-12RLI data, will never call AMPHIT with an ; undefined LUN. ; 11-Jan-2003 DG ; Changes to allow use of non-zero PHI_AB value. ;- pro amphit,frq,data,bl,tavg,header,cfg,out,pslope,ifrq=ifrq,fitslope=fitslope,$ debug=debug,save=save,lun=lun,pol=pol,ndfac=ndfac nant = header.nant ; Get index of active antennas IANTS (e.g. [1,2] => [0,1], [4,5,6] => [2,3,4] antlist = [1,2,4,5,6,7,8] ; Translation from antenna number to index aa = header.aatab(where(header.aatab ne 255)) iants = intarr(nant) for i = 0, nant-1 do begin iants[i] = where(antlist eq aa[i]) ; Antenna index endfor antn = aa[iants] ; Antenna number ; Get indexes of the two antennas, for determination of noise diode factor ant1 = nint(bl/10) ; Antenna number of first antenna iant1 = where(aa eq ant1,na1) ; Index of first antenna if (na1 eq 0) then begin print,'AMPHIT: Invalid baseline specified. Ant '+string(ant1,format='(i1)')+' not in antenna list.' return endif ant2 = bl mod 10 ; Antenna number of second antenna iant2 = where(aa eq ant2,na2) ; Index of second antenna if (na2 eq 0) then begin print,'AMPHIT: Invalid baseline specified. Ant '+string(ant2,format='(i1)')+' not in antenna list.' return endif ; Read the channel phases from cfg.chanphz for i = 0, 7 do cfg.muxarr[i,i] = -2 ; Flag none-baseline (total power) AB = where(cfg.muxarr eq 0, nab) if (nab eq 0) then begin ab = 0. print,'AMPHIT: Warning, file contains no AB baseline?' endif else phi_ab = cfg.chanphz[ab[0]]/!dtor BC = where(cfg.muxarr eq -1,nbc) if (nbc eq 0) then begin bc = 0. print,'AMPHIT: Warning, file contains no BC baseline?' endif else phi_bc = cfg.chanphz[bc[0]]/!dtor CA = where(cfg.muxarr eq 1,nca) if (nca eq 0) then begin ca = 0. print,'AMPHIT: Warning, file contains no CA baseline?' endif else phi_ca = cfg.chanphz[ca[0]]/!dtor ; Create the string containing the order of antennas, which backend channels ; were used, and current backend channel phase offset order = cfg.bldir[iants[iant1],iants[iant2]] ; 1 = direct, -1 = inverse bechan = cfg.muxarr[iants[iant1],iants[iant2]] ; 0 = AB, 1 = CA, -1 = BC chphz = cfg.chanphz[iants[iant1],iants[iant2]]/!dtor ; convert to degrees bechanstr = ['BC','AB','CA'] infostr = 'Baseline: '+bechanstr[bechan+1]+' Order: '+string(order,format='(I2)')+' Chan Phase: '+string(chphz,format='(F6.2)') nndw = 1 + header.ndbw + header.nfgw ; Number of non-data words in sample idex = header.lchoff[iants[iant1],iants[iant2]]-nndw ; CASE nant OF ; 3: begin ; bl_list = [12, 14, 24] ; A guess ; iant2 = ((bl mod 10) - 2)>1 ; Index of second antenna (1=Ant2, 2-5=Ant4-Ant7) ; if (iant2 ne 1) then bl_list = [12,10+iant2+2,20+iant2+2] ; Fix the guess ; if (iant2 gt 1) then iant2 = 2 ; end ; 6: begin ; bl_list = [12, 14,15,16,17, 24,25,26,27] ; iant2 = ((bl mod 10) - 2)>1 ; end ; 7: begin ; bl_list = [12,13,14,15,16,17,23,24,25,26,27] ; iant2 = (bl mod 10) - 1 ; end ; ELSE: begin ; bl_list = [12, 14,15,16, 24,25,26] ; iant2 = ((bl mod 10) - 2)>1 ; Antenna index of second antenna ; end ; ENDCASE ; idex = (where(bl_list eq bl,ndex))[0] ; if (ndex eq 0) then begin ; print,'Invalid baseline specified.' ; return ; endif if (n_elements(pslope) eq 0) then pslope = 0 if (keyword_set(ifrq)) then f = frq(ifrq-1) else f = frq if (not keyword_set(pol)) then pol = 'I' ; If no NDFAC was given, set to nominal unity. if (n_elements(ndfac) eq 0) then ndfac = replicate(1.0,nant) ; ND factor for this baseline nd = sqrt(ndfac[iant1]*ndfac[iant2]) ; remove any points with zero time (which can occur at the end ; of the scans) good = where(tavg ne 0,ngood) if (ngood gt 0) then begin tavg = tavg(good) data = data(*,*,good) endif junk = ' ' nf = n_elements(f) out = fltarr(6,nf) time = tavg/1000. ; convert to seconds tref = total(time/60.)/n_elements(time) if (not keyword_set(fitslope)) then tref = -tref if (keyword_set(debug)) then begin if (debug eq 1) then begin nx = 2 if (nf gt 10) then nx = ceil(nf/5) ny = ceil(float(nf)/nx) psav = !p !p.multi=[0,nx,ny,0,0] print,bl,' baseline [',pslope,']',pol,format='(I2,A11,F5.1,A1,1X,A1)' print,infostr print,' GHz Amp dAmp Pha dPha pslope' print,'---- ------- ------- ------ ------ ------' endif endif if (keyword_set(save)) then begin if (n_elements(lun) eq 0) then begin filename = dialog_pickfile(path=!defaults.workdir,/write) ans = 'No' if (filename ne '') then ans = oktowrite(filename) if (ans eq 'No') then return openw,/get_lun,lun,filename printf,lun,nf,3,1,phi_ab,phi_bc,phi_ca endif printf,lun,bl,' baseline [',pslope,']',pol,format='(I2,A11,F5.1,A1,1X,A1)' printf,lun,infostr printf,lun,' GHz Amp dAmp Pha dPha pslope' printf,lun,'---- ------- ------- ------ ------ ------' endif for i = 0, nf-1 do begin if (keyword_set(ifrq)) then begin x = reform(data(idex+1,ifrq-1,*)) y = reform(data(idex,ifrq-1,*)) endif else begin x = reform(data(idex+1,i,*)) y = reform(data(idex,i,*)) endelse good = where(finite(y),ngood) if (ngood ge 3) then begin xp = x(good) yp = y(good) t = time(good)/60. a = sqrt(xp^2+yp^2) p = atan(yp,xp)/!dtor - pslope*((t-abs(tref))/60.)*f(i) xp = a*cos(p*!dtor) yp = a*sin(p*!dtor) res = rvav(xp,yp,t,tref,aref,pref,pdot,siga,sigp) ; RVAV does a vector average of amplitudes, but if S/N of ; individual measurements is greater than about 3, a straight ; algebraic average is more appropriate. avg = total(a)/ngood s2n = avg/siga/sqrt(ngood) if (s2n gt 2.5) then aref = avg ; Ensure that phases are between 0 and 360 pref = lobe(pref) aref = aref/nd ; Scale amplitudes to ND factor siga = siga/nd ; Scale amp S. DEV. to ND factor if (keyword_set(debug)) then begin plot,t,lobe(p),psym=3,yran=[0,360] oplot,t,lobe(pref+pdot*(t-tref)/60.) endif if (keyword_set(save)) then begin printf,lun,f(i),aref,siga,pref,sigp,pdot/f(i),format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))' endif if (keyword_set(debug)) then begin if (debug eq 1) then print,f(i),aref,siga,pref,sigp,pdot/f(i),$ format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))' endif out(*,i) = [f(i),aref,siga,pref,sigp,pdot/f(i)] endif else begin if (keyword_set(save)) then begin printf,lun,f(i),-99.9,-99.9,-99.9,-99.9,-99.9,format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))' endif if (keyword_set(debug)) then begin if (debug eq 1) then print,f(i),-99.9,-99.9,-99.9,-99.9,-99.9,$ format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))' endif out(*,i) = [f[i],replicate(!values.f_nan,5)] endelse endfor if (keyword_set(debug)) then if (debug eq 1) then !p = psav return end