pro tmpphz,blin bl = [12,14,15,16,24,25,26] ibl = where(bl eq blin,nbl) if (nbl eq 0) then begin print,'Baseline entry, ',blin,' does not exist.' return endif ibl = ibl(0) file = dialog_pickfile(path=!defaults.workdir,filter='*.amf',$ title='Select multiscan .AMF file') ; Read the amphit file data for the R,L,I data ; in = {hed:strarr(3),data:fltarr(6,43)} ; openr,lun,/get_lun,file1 ; for i = 0, ibl(0) do readf,lun,in ; free_lun,lun ; dat = in.data ; Determine which entries have s/n better than 3, ; and save only those entries ; igood = where(dat[1,*]/dat[2,*] gt 2,ngood) ; data = fltarr(6,ngood) ; for i = 0, ngood-1 do data[*,i] = dat[*,igood[i]] ; Read the "current" calibration data from !DEFAULTS.EPHEMDIR+'AMPCAL.REC' and ; C:\EPHEM\PHZCAL.REC (must read the former in order to have the ; S/N for quality control purposes openr,lun,/get_lun,!defaults.ephemdir+'ampcal.rec' a = assoc(lun,uintarr(1024)) ampdat = [(a(0))[36:1015],(a(1))[12:*]] fampdat = ampdat*0. good = where(ampdat ne 'FFFF'x,ngood) if (ngood gt 0) then fampdat[good] = ampdat[good]*0.01 bad = where(ampdat eq 'FFFF'x,nbad) if (nbad gt 0) then fampdat[bad] = !values.f_nan free_lun,lun openr,lun,/get_lun,!defaults.ephemdir+'phzcal.rec' a = assoc(lun,uintarr(1024)) phzdat = (a(0))(36:*) fphzdat = phzdat*0. good = where(phzdat ne 'FFFF'x,ngood) if (ngood gt 0) then fphzdat[good] = phzdat[good]*0.01 bad = where(phzdat eq 'FFFF'x,nbad) if (nbad gt 0) then fphzdat[bad] = !values.f_nan free_lun,lun dat = fltarr(6,86) idx = indgen(86) dat[0,*] = (idx+5)*0.2 ; Frequency list case 1 of IBL lt 4: begin dat[1,*] = fampdat[idx*20]*fampdat[idx*20+3+ibl] dat[2,*] = fampdat[idx*20+10]*fampdat[idx*20+13+ibl] dat[3,*] = fphzdat[idx*10+3+ibl]-fphzdat[idx*10] end IBL ge 4: begin dat[1,*] = fampdat[idx*20+3]*fampdat[idx*20+ibl] dat[2,*] = fampdat[idx*20+13]*fampdat[idx*20+13+ibl] dat[3,*] = fphzdat[idx*10+3]-fphzdat[idx*10+3+ibl-3]-135.0 end endcase ; Determine which entries have s/n better than 2, ; and save only those entries igood = where(finite(dat[1,*]),ngood) dat = dat[*,igood] igood = where(dat[1,*]/dat[2,*] gt 2,ngood) ; data = fltarr(6,ngood) ; for i = 0, ngood-1 do data[*,i] = dat[*,igood[i]] data = dat[*,igood] idat = fltarr(6,10,5) ; Read the amphit file data for the 5 scans in = {hed:strarr(3),data:fltarr(6,10)} openr,lun,/get_lun,file for i = 0, ibl(0) do readf,lun,in idat[*,*,0] = in.data for j = 1, 4 do begin for i = 0, 6 do readf,lun,in idat[*,*,j] = in.data endfor free_lun,lun for k = 0, 4 do begin j = -1 f = (p = (dp = fltarr(10))) for i = 0, 9 do begin ifnd = where(abs(data[0,*] - idat[0,i,k]) lt 0.001,nfnd) if (nfnd eq 1) then begin j = j + 1 f[j] = idat[0,i,k] p[j] = ((idat[3,i,k]-data[3,ifnd]+18180.) mod 360.) - 180. sn = idat[1,i,k]/idat[2,i,k] dp[j] = (1/(sn*!dtor))>5. endif endfor if (j lt 2) then begin print,'Baseline ',blin,' scan ',k+1,' has too few frequencies.' goto,skip endif f = f[0:j] p = p[0:j] dp = dp[0:j] nf = n_elements(f) junk = '' manual = 1 window,0,xsiz=450,ysiz=600 junk = '' outstr = '' manual = 0 nit = 0 ap = fltarr(nf) again: nit = nit + 1 res = linfit(f,p+ap*360,sdev=dp,chi=chi,sigma=sigma) ; res2 = linfit(f,p+ap*360,sdev=dp,chi=chi2) if (manual) then begin ovsa_ploterr,f,p+ap*360.,dp,psym=-2,xran=[0,20],xsty=1,$ tit=' Chisq = '+string(chi,format='(f7.1)') oplot,[0,f,20],res(0)+res(1)*[0,f,20] ; oplot,[0,f,12],res2(0)+res2(1)*[0,f,12],linesty=2 read,'Adjustments (all 0=>cont): ',ix,dx x = x*0 x[ix:nf-1] = dx adj = where(x ne 0,nadj) ap = ap + x if (nadj gt 0) then goto, again endif else begin x = intarr(nf) diff = [p+ap*360,0] - [0,p+ap*360] ; See if we need to lower any phases ilower = where(diff[0:nf-1] gt 180,nlower) if (nlower ne 0) then begin ; Yes, we have to lower at least one, so lower first x[ilower[0]:nf-1] = -1 endif else begin ; We did not have to lower any, so see if we ; have to raise any. iraise = where(diff[0:nf-1] lt -180,nraise) if (nraise ne 0) then x[iraise[0]:nf-1] = 1 endelse adj = where(x ne 0,nadj) ap = ap + x if (nadj gt 0 and nit lt 20) then goto, again ovsa_ploterr,f,p+ap*360.,dp,psym=-2,xran=[0,20],xsty=1,$ tit=' Chisq = '+string(chi,format='(f7.1)') oplot,[0,f,20],res(0)+res(1)*[0,f,20] ; oplot,[0,f,12],res2(0)+res2(1)*[0,f,12],linesty=2 ; read,'CR to continue: ',junk manual = 1 ; If excessive chi-square, ask for manual verification if (chi gt 30) then goto, again endelse ; print,strmid(file,8,5),tim,res[0]+res[1],chi,' '+bl(i),format='(A5,f8.0,f8.1,f8.0,A7)' outstr = string(res[0],res[1],chi,format='(2f7.1,f6.0)') ; printf,outlun,strmid(file,8,5),tim,outstr,format='(A5,f8.0,A)' print,blin,k+1,outstr,format='("Baseline ",I2," scan ",I2," ",A)' ; free_lun,outlun read,'CR to continue: ',junk skip: endfor return end