;+ ; NAME: ; pcal ; PURPOSE: ; Phase calibration data fitting as a linear function of frequency. ; CATEGORY: ; scbong OVRO ; CALLING SEQUENCE: ; pcal, files, savefile[, pfold = pfold][, /pclose] ; INPUTS: ; files String vector of file names to be read. If not set, the user is prompted for files. ; savefile File name to be saved. If not set, "pcal.sav" is used. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; antset Antenna numbers of sets to be shown. ; This is [3, nset] dimensional array where nset is the number of sets. ; A set consists of three antennas which should satisfy the phase closure condition. ; The default value is [[1, 2, 4], [1, 2, 5], [1, 2, 6]]. ; chisq If pclose keyword is set, chi square value is returned. ; nowin Set this keyword if plotting is not wanted. ; pclose If set, fit will satisfy the phase closure condition so that, ; for any three antennas a, b, c in increasing antenna number order, ; ab - ac + bc = 0 where ab is the phase of baseline a-b and so on. ; If pfold keyword is not given properly, fitting line is found in the whole range. ; If pfold keyword is given with the right dimension, the fitting line ; is found in the vicinity of given points. ; pfold Phase fold factors. This should have [nbl, nfrq, nfl] dimension, ; where nbl = nant * (nant - 1) / 2 and nant is the number of antennas, ; nfrq is the number of frequencies observed and nfl is the number of files to be read. ; The fit will be done with the phases to which pfold times 360 degrees is added. ; The default value is 0 for all the elements. ; If this has a wrong dimension and pclose is not set, default value is returned to this. ; If pclose is set, the values found in phase closure is returned to this. ; zero If set, fit will be done with zero offset which is what it should be. ; ROUTINES CALLED: ; plot_pcal ; OUTPUTS: ; ; COMMENTS: ; plot_pcal will show the same number of windows as that of files used in pcal. ; A window has 0 marked points which are the points of 0 phase fold. ; In the last row plots, three dashed lines are fitted lines of ; each labeled baseline index and the solid line is the sum of them. ; dashed lines and the data points are flipped if the sign next to the label is '-', ; so that the sum of them becomes 0 according to the phase closure condition. ; SIDE EFFECTS: ; Common block named pcal_pcfunc_data is set. ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 22-Jan-2000 by Su-Chan Bong ;- function pcal_linfunc, f, npm if npm eq 1 then return, [f] else return, [f, 1.] end function pcal_pcfunc, pm common pcal_pcfunc_data, pdat, sdat, f, pfdat, nant1, nant2, zero, pma, pmb, iant, jant, one, express pmb[1 : *] = pm[0 : nant2] if express then begin b = pmb[jant] - pmb[iant] adelp = (((b # f - pdat) mod 360.) + 360.) mod 360. chisq = total(((adelp < (360. - adelp)) / sdat)^2, /nan) endif else begin if not zero then pma[1 : *] = pm[nant1 : *] b = pmb[jant] - pmb[iant] a = pma[jant] - pma[iant] delp = a # one + b # f - pdat good = where(finite(pdat) and finite(sdat), ngood) if ngood gt 0 then pfdat[good] = round(delp[good] / 360.) chisq = total(((delp - 360. * pfdat) / sdat)^2, /nan) endelse return, chisq end pro pcal, files_, savefile_, pfold = pfold, pclose = pclose_, antset = antset, nowin = nowin, $ zero = zero_, chisq = chisq common pcal_pcfunc_data, pdat, sdat, f, pfdat, nant1, nant2, zero, pma, pmb, iant, jant, one, express if keyword_set(files_) then files = files_ else files=dialog_pickfile(path=!defaults.workdir,filter='*.sav',/multi) if (files[0] eq '') then return if keyword_set(savefile_) then savefile = savefile_ else savefile = 'pcal.sav' if keyword_set(pclose_) then pclose = 1b else pclose = 0b if keyword_set(zero_) then zero = 1b else zero = 0b nfl = n_elements(files) nfl1 = nfl - 1 ;Restore the first file to determine the dimention. ;tavg, avg, sigma, f, tls, refcal, uvt, traj, geometry, polcode are restored. restore, files[0] refcal_phase = *refcal.pphase savg = size(avg) nch = savg[1] nant = fix(sqrt(nch)) nant1 = nant - 1 nant2 = nant - 2 nchbl = nch - nant nchbl1 = nchbl - 1 nbl = nant * (nant1) / 2 nbl1 = nbl - 1 nfq = savg[2] nfq1 = nfq - 1 ntm = savg[3] dflt = 1b if keyword_set(pfold) then begin spfold = size(pfold) if spfold[0] eq 3 then $ if (spfold[1] eq nbl and spfold[2] eq nfq and spfold[3] eq nfl) then dflt = 0b endif if dflt then pfold = fltarr(nbl, nfq, nfl) pshift = pfold * 360. ichbl = indgen(nchbl) + nant avgavg = fltarr(nchbl, nfq, nfl) varavg = avgavg ttab = dblarr(nfl) ;Average over the time and obtain the phase. for ifl = 0, nfl1 do begin if ifl gt 0 then restore, files[ifl] ntag = n_tags(refcal) for i = 0, ntag - 1 do free_pointer, refcal.(i) ntag = n_tags(traj) for i = 0, ntag - 1 do free_pointer, traj.(i) ;Because time returns to 0 at 24h, UT, time correction should be done before average. itavg_nd = where(tavg lt 4.32d7, nitavg_nd) if (nitavg_nd gt 0) and (nitavg_nd lt n_elements(tavg)) then begin tavg[itavg_nd] = tavg[itavg_nd] + 8.64d7 ttab[ifl] = (moment(tavg, /double, maxmoment = 1))[0] mod 8.64d7 endif else ttab[ifl] = (moment(tavg, /double, maxmoment = 1))[0] for ifq = 0, nfq1 do for ich = 0, nchbl1 do begin ;moment requires more then 1 finite values. avg2mnt = avg[ichbl[ich], ifq, *] iavg2mnt_f = where(finite(avg2mnt), niavg2mnt_f) case niavg2mnt_f of 0: begin avgavg[ich, ifq, ifl] = !values.f_nan varavg[ich, ifq, ifl] = !values.f_nan end 1: begin avgavg[ich, ifq, ifl] = avg2mnt[iavg2mnt_f[0]] varavg[ich, ifq, ifl] = !values.f_nan end else: begin mnt = moment(avg2mnt[iavg2mnt_f], maxmoment = 2) avgavg[ich, ifq, ifl] = mnt[0] varavg[ich, ifq, ifl] = mnt[1] end endcase endfor endfor ichsin = indgen(nbl) * 2 ichcos = ichsin + 1 avgsin = avgavg[ichsin, *, *] avgcos = avgavg[ichcos, *, *] ptab = atan(avgsin, avgcos) * !radeg ;Define standard deviation of the phases as 1/(a*sp)^2 = (sinp/sx)^2 + (cosp/sy)^2, ;where a is the amplitude, p is the phase, and s is the standard deviation, ;so that sp = 1 / sqrt((y/sx)^2 + (x/sy)^2) stab = 1. / sqrt(avgsin^2 / varavg[ichcos, *, *] + avgcos^2 / varavg[ichsin, *, *]) p000 = reform(refcal_phase[2, *, *]) for ifl = 0, nfl1 do ptab[*, *, ifl] = ptab[*, *, ifl] - p000 ptab = (ptab + 720.) mod 360. p2fit = ptab + pshift ;Fit with the given points. if zero then npm = 1 else npm = 2 pm = fltarr(nbl, nfl, npm) for ifl = 0, nfl1 do for ibl = 0, nbl1 do begin p2f = p2fit[ibl, *, ifl] ip2f = where(finite(p2f), nip2f) if nip2f ge 3 then pm[ibl, ifl, *] = svdfit(f[ip2f], p2f[ip2f], npm, $ function_name = 'pcal_linfunc', weights = 1. / stab[ibl, ip2f, ifl]^2) $ else pm[ibl, ifl, *] = !values.f_nan endfor btab = pm[*, *, 0] if zero then atab = fltarr(nbl, nfl) else atab = pm[*, *, 1] ;Phase closure using 1-b baselines. ;From 1-b baselines, we obtain the phases, of the preferred linear form, ;for each antenna except for antenna 1 which is set to 0. ;Using these parameters as the initial values, the optimal parameters ;which minimize the chi sqare value is obtained. ;Then from these, we reproduce the parameters of the baselines. if pclose then begin bscl = 180. / f[nfq1] pmmax = 200. pmmin = -pmmax pmini = fltarr(nant1) + pmmin pmrange = 2. * pmmax pmrange2 = pmrange^2 pmb = fltarr(nant) pma = pmb pmscl = fltarr(nant1) + bscl if not zero then pmscl = [pmscl, fltarr(nant1) + 30.] nant3 = nant - 3 ipmmax = nant - 4 > 0 ntry = fix(pmrange / bscl) ntry1 = ntry - 1 chisqtry = dblarr(ntry, ntry) pmtry = findgen(ntry) * bscl - pmmax lablend = '/' + strtrim(string(nfl), 1) + ')' iant = intarr(nbl) jant = iant ibl = 0 for i = 0, nant2 do for j = i + 1, nant1 do begin iant[ibl] = i jant[ibl] = j ibl = ibl + 1 endfor one = fltarr(nfq) + 1. chisq = dblarr(nfl) for ifl = 0, nfl1 do begin pdat = ptab[*, *, ifl] sdat = stab[*, *, ifl] pfdat = pfold[*, *, ifl] if dflt then begin ;Find the initial parameters for the phase closure. ;Parameters are [b1, ..., bnant1, a1, ..., anant1] where bi is the 1st order factor ;of antenna index i and ai is the zeroth order factor of antenna index i. ;If zero keyword is set, a's are omitted. ;For the initial parameters, we find only b's in the local depression of chi sqaure where ;global minimum value is located. express = 1b pm = pmini pm0 = pm ipm = ipmmax go = 1b chisqmin = pcal_pcfunc(pm) proglable = 'Phase Closure (' + strtrim(string(ifl + 1), 1) + lablend id = progmeter(/init, label = proglable, buttontext = 'Cancel') while go do begin for itry = 0, ntry1 do for jtry = 0, ntry1 do begin pm[nant2] = pmtry[jtry] pm[nant3] = pmtry[itry] chisqtry[jtry, itry] = pcal_pcfunc(pm) endfor chisq0 = min(chisqtry, itry) if chisq0 lt chisqmin then begin chisqmin = chisq0 pm0 = pm pm0[nant2] = pmtry[itry mod ntry] pm0[nant3] = pmtry[itry / ntry] endif pm[ipm] = pm[ipm] + bscl if (pm[ipm] gt pmmax) then begin while (pm[ipm] gt pmmax) and go do begin pm[ipm] = pmmin ipm = ipm - 1 if ipm lt 0 then begin go = 0b ipm = ipmmax endif else pm[ipm] = pm[ipm] + bscl endwhile if (progmeter(id, (pm[0] < pmmax) / pmrange + 0.5 + $ bscl * ((pm[1] < pmmax) + pmmax) / pmrange2) eq 'Cancel') then begin val = progmeter(id, /destroy) print, 'Phase closure was canceled by user.' goto, finish endif ipm = ipmmax endif endwhile if not zero then pm0 = [pm0, fltarr(nant1)] express = 0b val = progmeter(id, /destroy) endif else begin ;Set initial parameters from the fit lines. pm0 = btab[0 : nant2, ifl] if not zero then pm0 = [pm0, atab[0 : nant2, ifl]] ipm0_bad = where(1b - finite(pm0[0 : nant2]), nipm0_bad) if nipm0_bad gt 0 then begin pm0[ipm0_bad] = 0. if not zero then pm0[ipm0_bad + nant1] = 0. endif express = 0b endelse ;Find optimal parameters and obtain baseline fit parameters. pm = amoeba(0.0001, function_name = 'pcal_pcfunc', p0 = pm0, scale = pmscl) if n_elements(pm) gt 1 then begin chisq[ifl] = pcal_pcfunc(pm) pmb[1 : *] = pm[0 : nant2] if not zero then pma[1 : *] = pm[nant1 : *] btab[*, ifl] = pmb[jant] - pmb[iant] atab[*, ifl] = pma[jant] - pma[iant] pfold[*, *, ifl] = pfdat endif else print, 'Phase closure failure for file ' + strtrim(string(ifl), 1), '.' endfor endif finish: save, ttab, atab, btab, ptab, f, pfold, filename = savefile ;Plot results. if not keyword_set(nowin) then plot_pcal, savefile, antset = antset end