;+ ; NAME: ; PCLOSE ; PURPOSE: ; Force phase closure in general case. ; CATEGORY: ; OVRO APC CALIBRATION ; CALLING SEQUENCE: ; pout = pclose(pin,wts[,cor=cor][,/nolobe]) ; INPUTS: ; pin A float array of phases of size (nb) where NB=7 is the number of ; baselines for which the phase closure is to be calculated. ; wts A float array of weight factors (basically 1/sn where sn = signal ; to noise ratio), that indicates how accurately the phases were ; determined. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; cor A keyword to be set to a named variable, which will on output ; contain the nb-element phase correction array. ; nolobe A switch that, if set, causes lobe adjustment (forcing ; closure result between -180 and 180) to be skipped. This is ; necessary when phase closure for phase slope is desired, ; since a closure to 360 degrees is very different from a ; closure to 0 degrees for phase slope! ; ROUTINES CALLED: ; OUTPUTS: ; pout An array of the same size as pin, but with phases corrected ; for phase closure. ; COMMENTS: ; The algorithm for phase closure is a bit complex. The task is to ; adjust the phases within the errors to take advantage of the fact ; that we know the phases on each triad of baselines must close. Since ; the large-large baseline is much more sensitive that the large-small, ; and the small-small baselines are not used, the algorithm needs to ; first adjust the large-large (one baseline of each triad) to be ; mutually consistent, then adjust each of the large-small according ; to their errors. ; ; Define: ; closure phase: ; c_1 = p_14 - p_24 - p_12 - 225 (124 triad) ; c_2 = p_15 - p_25 - p_12 - 225 (125 triad) ; c_3 = p_16 - p_26 - p_12 - 225 (126 triad) ; weights: ; w_12 = 1/sn_12 (inverse of S/N in amplitude, baseline 12) ; w_2j = 1/sn_2j (BC baseline for triad j) ; w_3j = 1/sn_3j (CA baseline for triad j) ; fractions: ; f_1j = w_12 / (w_12 + w_2j + w_3j) (triad j = 1, 3) ; ; Then we can calculate three phase corrections for baseline 12: ; C12_j = f_1j*c_j (triad j = 1, 3) ; but we only want one, so take the average ; C12 = Sum(C12_j) / 3 (triad j = 1, 3) ; and correct the closure phase accordingly ; c_j' = c_j + C12 (triad j = 1, 3) ; ; Now determine fractions: ; f_2j = w_2j / (w_2j + w_3j) ; f_3j = w_3j / (w_2j + w_3j) ; and determine final phase corrections ; C14 = f_21*c_1' ; C15 = f_22*c_2' ; C16 = f_23*c_3' ; C24 = f_31*c_1' ; C25 = f_32*c_2' ; C26 = f_33*c_3' ; ; Finally, adjust all of the input phases according to these ; corrections, applying the appropriate signs ; p12' = p12-C12 ; p1j' = p1j-C1j ; p2j' = p2j+C2j ; ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 01-Jul-1999 by Dale E. Gary ; 08-Jul-1999 DG ; Attempt to deal with missing data (indicated by NaN) by artificially ; setting the weight to a large value. ; 06-Apr-2000 DG ; Treat case of phase being non-NaN, but weight being NaN. ; 09-May-2000 DG ; Added NOLOBE keyword, so skip lobe adjustment (used when the closure ; is sought for phase slope). ; 16-Jun-2000 DG ; Made the code general for 5, 6, or 7 antennas. May work with 4, but ; probably not for 3 antennas. ; 27-Sep-2000 DG ; Changed routine to reflect new, simplified definition of BC baselines. ; The sign of the BCsin component is reversed in GAINCOR, which allows ; us to not worry about special treatment of BC baselines here. ;- function pclose,pzin,wtin,offset,cor=Cor,nolobe=nolobe wts = wtin ; Make copies of variables that I can change without changing pin = pzin ; the inputs nbl = n_elements(pzin) Cor = fltarr(nbl) ; Deal with missing data, as indicated by NaN bad = where(finite(pin) eq 0,nbad) if (nbad gt 0) then begin ; Simply set the weights of such data to a large number (100) and ; overwrite the NaN with 0. NB: This will work fine for one missing ; phase of a triad. For two or three missing phases, the error in ; phase will be distributed among the missing values, which means ; neither (or none) is likely to be correct. wts[bad] = 100. pin[bad] = 0. endif else begin bad = where(finite(wts) eq 0, nbad) ; If phases are okay, but wts are NaN, set weights to large number (100). if (nbad gt 0) then wts[bad] = 100. endelse ; This section of code should be made more general--will break for 7 dishes ; Determine array ranges for large-small baselines n1s = 1 ; Start index of 1-n baselines (always 1) n1e = (nbl-1)/2 ; End index of 1-n baselines n2s = n1e+1 ; Start index of 2-n baselines (one more than n1e) n2e = nbl-1 ; End index of 2-n baselines nsmall = n1e ; Number of small baselines for each large antenna ; Closure phases (array) c = reform(-pin[n2s:n2e] + pin[n1s:n1e] - pin[0] - offset) ; Ensure that the closure phases are between -180 and 180, unless ; NOLOBE is set if (not keyword_set(nolobe)) then c = lobe(c,/mid) ; Fractions f = fltarr(3,nsmall) for i = 0, nsmall-1 do begin f[0,i] = wts[0]/total(wts[[0,n1s+i,n2s+i]]) f[1,i] = wts[n1s+i]/(wts[n1s+i]+wts[n2s+i]) f[2,i] = wts[n2s+i]/(wts[n1s+i]+wts[n2s+i]) endfor ; Baseline 12 correction (average of measurements) Cor[0] = total(reform(f[0,*])*c)/float(nsmall) ; Correct the closure phases for C12 correction cp = c + Cor[0] ; Determine other baseline corrections, with proper sign on baselines ; with antenna 1. Cor[n1s:n1e] = reform(f[1,*])*cp Cor[n2s:n2e] = -reform(f[2,*])*cp ; Correct the phases for all of the baselines (with special sign ; convention for 6-element mux mode--see above) pout = pin - Cor return,pout end