;+ ; NAME: ; BASFIT ; PURPOSE: ; Given a list of sources and corresponding 1 GHz phases ; determined by BASEPHZ routine, determines phase change for ; slightly different baseline lengths in a 3D matrix, applies ; the phase change to data from each baseline, and locates ; the set of baseline shifts dBx, dBy, dBz that give the ; smallest residuals. ; CATEGORY: ; OVRO APC INSTRUMENT CALIBRATION ; CALLING SEQUENCE: ; resid = basfit(infile,amfile,minphz[,stepsize][,blcor=blcor][,nsteps=nsteps]) ; INPUTS: ; infile the name of a file that contains the list of source ; data. Consists of one line for each source measurement. ; Each line has filename, scan header record, time, and list ; of 1 GHz phases (from BASEPHZ routine) and chi-square, ; one pair for each baseline (14 numbers all together). ; amfile the routine has been changed to ignore most of the contents ; of infile (although the filename, hrec, and time are ; still important), and instead use the output of AMPHIT ; (specifically WAMPHIT) to fit all of the frequencies ; instead of simply a best fit line. ; stepsize the size of the step [cm] between adjacent points in the ; 3D search space. 3 cm is a good size to start with. ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; blcor an array of corrections from a previous run of BASFIT ; (as returned by VIS_BASFIT routine), to apply before ; starting the search in phase space. Of size (3,NBL) ; where the first index if for the three axes Bx, By, ; and Bz, and NBL is the number of baselines. ; nsteps either a single value or a 3-element array giving the ; number of steps IN EACH DIRECTION to include in the ; search. The actual number of steps is NSTEPS*2 + 1 in ; each dimension. If a single value, the same number ; of steps is used in each axis. If a 3-element array ; is given, they are interpreted as NX,NY,NZ. If omitted, ; NX = NY = NZ = 8 is used as the default. ; ; ROUTINES CALLED: ; OUTPUTS: ; resid the giant array of residuals, of size ; (NSTEPS*2+1,NSTEPS*2+1,NSTEPS*2+1,NBL), ; where the indexes are BX offset, BY offset, BZ offset (all ; ranging from -NSTEPS to NSTEPS in units of the stepsize), ; and the baseline number. ; minphz an array of size (NANT,NANT,NMEAS) containing both the ; original phase information (in the bottom diagonal part of ; the NANT x NANT matrix) and the phases for the set of ; baselines with the minimum residuals (upper diagonal part ; of the NANT x NANT matrix) ; COMMENTS: ; The residuals and a comparison of minimum phases with input phases ; can be viewed using VIS_BASFIT routine. ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 01-Jul-1999 by Dale E. Gary ; 05-July-2000 DG ; Several changes to make it work with any number of antennas (maybe), and ; for a settable number of steps in each axis. Also now uses all of the ; phases directly rather than the phase slope fit by BASEPHZ (via the AMFILE ; input file). ;- ;----------------------------------------------------------------------- function basfit_resid,phz,wt n = n_elements(phz[0,0,*,0]) nf = n_elements(phz[0,0,0,*]) phzdif = ((phz + 18180.) mod 360.) - 180. ; Calculate the phase differences between adjacent sources (allows ; for a solution with a fixed or slowly varying offset) phzdif = ((shift(phz,0,0,1,0) - phz + 18180.) mod 360.) - 180. ; Eliminate first phzdif and corresponding weight, since it involves ; a difference of phases at the ends of the observation phzdif = phzdif[*,*,1:n-1,*] resid = total(total(phzdif*phzdif*wt,3)/(total(wt,3)>1),3)/nf ; resid = total(phzdif*phzdif,3)/n return,resid end ;----------------------------------------------------------------------- function basfit,infile,amfile,minphz,stepsz,blcor=blcor,nsteps=nsteps NMAX = 100 ; Maximum number of source measurements allowed ; Set the number of steps from NSTEPS keyword. Default is 8,8,8. If ; only one value given, all three are assumed to be the same. If three ; values are given then the three can be different. if (not keyword_set(nsteps)) then nsteps = 8 if (n_elements(nsteps) eq 3) then begin nstepsx = nsteps[0] nstepsy = nsteps[1] nstepsz = nsteps[2] endif else nstepsx = (nstepsy = (nstepsz = nsteps)) f = [1.8,2.6,3.8,4.4,5.0,7.2] ; Frequency list nf = n_elements(f) ; Open input file CATCH,err_stat if (err_stat eq 0) then begin openr,/get_lun,lun,infile CATCH,/CANCEL endif else begin junk = dialog_message('BASFIT: Could not open input file '+infile,/ERROR) return,0 endelse ; Read contents of file. Contains 2 header lines and NMEAS source ; measurement lines, where NMEAS is arbitrary. ON_IOERROR,done head = strarr(2) line = ' ' filename = strarr(NMAX) hrec = fltarr(NMAX) tim = fltarr(NMAX) readf,lun,head ; Determine number of antennas from first line of header (this means the ; header must be correct!) nant = (n_elements(str_sep(strtrim(strcompress(head[0])),' ')) - 3)/2 nbl = 2*nant-3 linedat = fltarr(nbl*2) srcphz = fltarr(nant,nant,NMAX) srcchi = fltarr(nant,nant,NMAX) nmeas = 0 icnt = -1 while (not EOF(lun)) do begin readf,lun,line icnt = icnt + 1 ; Skip line if first char is ';' if (strmid(line,0,1) NE ';') then begin if (n_elements(ikeep) gt 0) then ikeep = [ikeep,icnt] else ikeep = icnt reads,strmid(line,12),hrc,t,linedat filename[nmeas] = strmid(line,0,12) hrec[nmeas] = hrc tim[nmeas] = t k = 0 for i = 0,1 do begin for j = i+1,nant-1 do begin srcphz[i,j,nmeas] = linedat[k] srcphz[j,i,nmeas] = linedat[k] srcchi[i,j,nmeas] = linedat[k+1] srcchi[j,i,nmeas] = linedat[k+1] k = k+2 endfor endfor nmeas = nmeas+1 if (nmeas eq NMAX) then begin ans = dialog_message(['BASFIT: Reached maximum number of sources.',$ 'No further lines will be read'],/info) goto, done endif endif endwhile done: free_lun,lun ON_IOERROR,NULL ; Truncate the arrays to only the size read in if (n_elements(ikeep) eq 0) then begin print,'BASFIT: Error, no valid entries?' return,0 endif filename = !DEFAULTS.DATADIR+filename[0:nmeas-1] ; Append data directory hrec = hrec[0:nmeas-1] tim = tim[0:nmeas-1] srcphz = srcphz[*,*,0:nmeas-1]/5.0 ; Convert to 1 GHz srcchi = srcchi[*,*,0:nmeas-1]/5.0 ; Convert to 1 GHz ; Loop over the input lines, opening each file in turn and reading ; the header and geometry information, to store for later use. id = progmeter(/init,button='Abort',label='Reading file(s). Percent Complete') for i = 0, nmeas-1 do begin ; Update progress meter ans = progmeter(id,(i+1.)/nmeas) if (ans eq 'Cancel') then begin free_lun,flun ans = progmeter(id,/destroy) return,0 endif skipit = 0 ; Only open the file if different than the previous one if (i eq 0) then begin flun = openarc(filename[i],a,nrec) endif else if (filename[i] ne filename[i-1]) then begin free_lun,flun flun = openarc(filename[i],a,nrec) endif else begin ; If this is the same file as the previous one, and the header ; record is the same, then we do not need to redo the geometry, ; so set a flag to skip NEWSCAN call. if (hrec[i] eq hrec[i-1]) then skipit = 1 endelse if (flun eq 0) then begin ans = DIALOG_MESSAGE('BASFIT: Error, data file '+filename[i]+' not found.',/error) return,0 endif ; Run the procedure to read header information from the scan and ; declare pointers to data arrays (unless we have already read it ; and can skip it (when SKIPIT flag is non-zero) if (skipit eq 0) then rec = newscan(a,nrec,hrec[i],header,obseq,cfg,traj,refcal,gparm,geometry,pRecent,pOld,cycle) if (rec eq -1) then begin ; NEWSCAN returned an error ans = DIALOG_MESSAGE('BASFIT: Error in file structure or missing segment.',/error) free_lun,flun return,0 endif if (i eq 0) then begin geo = replicate(geometry,nmeas) srcname = strarr(nmeas) endif srcname[i] = traj.srcname geo[i].bx = geometry.bx geo[i].by = geometry.by geo[i].bz = geometry.bz geo[i].htdiff = geometry.htdiff geo[i].a_ra = geometry.a_ra geo[i].a_dec = geometry.a_dec geo[i].year = geometry.year geo[i].doy = geometry.doy geo[i].dut = geometry.dut endfor free_lun,flun ans = progmeter(id,/destroy) ; Read in the AMF file containing the data ON_IOERROR,datdone openr,extlun,/get_lun,amfile readf,extlun,nf1,nbl1,nmeas1 in = {hed:strarr(4),data:fltarr(6,nf1)} datin = fltarr(6,nf1,nbl1,nmeas1) bl = intarr(nbl) for j = 0, nmeas1-1 do begin for i = 0, nbl1-1 do begin readf,extlun,in bl[i] = fix(strmid(in.hed[0],0,2)) datin[*,*,i,j] = in.data endfor endfor datdone: free_lun,extlun ON_IOERROR,NULL srcphz = fltarr(nant,nant,nmeas,nf) srcchi = srcphz for i = 0, 1 do begin for j = i+1,nant-1 do begin for k = 0, nf-1 do begin srcphz[i,j,*,k] = datin[3,k,i+j-1,ikeep] srcchi[i,j,*,k] = datin[2,k,i+j-1,ikeep]/datin[1,k,i+j-1,ikeep] ; 1/SN srcphz[j,i,*,k] = datin[3,k,i+j-1,ikeep] srcchi[j,i,*,k] = datin[2,k,i+j-1,ikeep]/datin[1,k,i+j-1,ikeep] ; 1/SN endfor endfor endfor ; Calculate residuals for baseline offsets in a 3D grid id = progmeter(/init,button='Abort',label='Calculating Residuals. Percent Complete') ;ifrq = [0,2,3,4] ;srcphz = srcphz[*,*,*,ifrq] ;srcchi = srcchi[*,*,*,ifrq] ;nf = n_elements(ifrq) ; These are the indexes of unused values in the nant x nant arrays blah = intarr(nant,nant)*0 + 1 blah[0,0] = (blah[1,1] = (blah[2:*,2:*] = 0)) zeroes = where(blah eq 0) ; Declare storage for phases after adjustment for offset nuphz = fltarr(nant,nant,nmeas,nf) ; Declare storage for output residuals resid = fltarr(nstepsx*2+1,nstepsy*2+1,nstepsz*2+1,nbl) ; Calculate weights, which are not a function of baseline offset ; Elimination of the first measurement is due to requirement of ; routine BASFIT_RESID, where the weights are used. wt = fltarr(nant,nant,nmeas-1,nf) wt = srcchi[*,*,1:nmeas-1,*] ; The weight is 1/chi-square, so determine which are the unused ; values (which are zero) to avoid taking 1/0. nonzero = where(wt ne 0.) wt[nonzero] = 1./wt[nonzero] ; Variables for finding and saving minimum residual info resmin = 100000.*replicate(1,nbl) ; For convenience, generate some indexes into NANT x NANT arrays indx = 0 jndx = 1 for i = 0,1 do begin for j = 2,nant-1 do begin indx = [indx,i] jndx = [jndx,j] endfor endfor minphz = srcphz ; Loop over all grid positions nstep = 0L ; Make baseline corrections if BLCOR keyword is set bxcor = (bycor = (bzcor = fltarr(nant,nant))) if (keyword_set(blcor)) then begin bxcor[indx,jndx] = (bxcor[jndx,indx] = blcor(0,*)) bycor[indx,jndx] = (bycor[jndx,indx] = blcor(1,*)) bzcor[indx,jndx] = (bzcor[jndx,indx] = blcor(2,*)) endif ; Convert input stepsz in cm to units of 1 GHz fringe separation if (n_elements(stepsz) ne 0) then stepsize = stepsz/30. else stepsize=0.1 ; Kernel (BC baselines are opposite sign from others) ;kern = [[0,1,1,1,1],[1,0,-1,-1,-1],[1,-1,0,1,1],[1,-1,1,0,-1],[1,-1,1,-1,0]] kern = replicate(1,nant,nant) for ibz = -nstepsz,nstepsz do begin iz = ibz+nstepsz ; Current Z baseline offsets (convert existing BZCOR in cm to nsec) bz = float(kern*ibz)*stepsize + bzcor/30. bz[zeroes] = 0. ; Zero the unused values for iby = -nstepsy,nstepsy do begin iy = iby+nstepsy ; Current Y baseline offsets (convert existing BYCOR in cm to nsec) by = float(kern*iby)*stepsize + bycor/30. by[zeroes] = 0. ; Zero the unused values for ibx = -nstepsx,nstepsx do begin ix = ibx+nstepsx ; Current X baseline offsets (convert existing BXCOR in cm to nsec) bx = float(kern*ibx)*stepsize + bxcor/30. bx[zeroes] = 0. ; Zero the unused values ; Update progress meter nstep = nstep + 1L ans = progmeter(id,nstep/((nstepsx*2+1)*(nstepsy*2+1.)*(nstepsz*2+1.))) if (ans eq 'Cancel') then begin free_lun,flun ans = progmeter(id,/destroy) return,0 endif ; Loop over source measurements for i = 0, nmeas-1 do begin ; Insert source geometry into geometry structure struct_assign,geo[i],geometry ; Override baselines with baseline offsets geometry.bx = bx geometry.by = by geometry.bz = bz geometry.htdiff = geometry.htdiff*0 ; Set height difference to zero ; Determine the delay for these offsets, in nsec delay,tim[i]*1000.,10000.,tau,dtau,geometry tau[zeroes] = 0. ; Zero the unused values ; Since light travels exactly 1 wavelength at 1 GHz ; in 1 nsec, the phaseshift we want, in degrees, is ; just the delay times 360 degrees ; Convert to 5 GHz reference by multiplying tau by 5.0 for j = 0,nf-1 do begin nuphz[*,*,i,j] = reform(srcphz[*,*,i,j]) - tau*360.*f[j] endfor endfor ; Now all we have to do is calculate how close these ; new phases are to zero, with appropriate weighting res = basfit_resid(nuphz,wt) resid[ix,iy,iz,*] = res[jndx,indx] ; Determine if any of these are the "best yet" ilt = where(reform(resid[ix,iy,iz,*]) lt resmin,nlt) if (nlt ne 0) then begin ; At least one baseline is the "best yet", so save ; some relevant info resmin[ilt] = resid[ix,iy,iz,ilt] minphz[indx[ilt],jndx[ilt],*,*] = nuphz[indx[ilt],jndx[ilt],*,*] endif endfor endfor endfor ans = progmeter(id, /destroy) return,resid end