;+ ; NAME: ; MAKE_APCAL ; PURPOSE: ; Interim routine to generate the AMPCAL and PHZCAL segment records from AMPHIT data ; obtained from Reference amplitude and phase calibration observations. The segment ; records are written into each scan by ARCHIVE. ; CATEGORY: ; OVRO APC DATA CALIBRATION ; CALLING SEQUENCE: ; make_apcal[,/s273] ; INPUTS: ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; S273 If set, the calibrator is assumed to be 3C273 and the 3C273 fluxes are ; applied. Otherwise the calibrator is assumed to be 3C84. ; ROUTINES CALLED: ; getparms, tlsnow, get_tl_struct, tl_encode, dayofyr ; OUTPUTS: ; COMMENTS: ; The values written to the file are in scaled, 2-byte integer form. To convert ; the amplitude values back to units/SFU, read the integer value, multiply by 0.01, ; and apply to baseline data by multiplying the factors for BOTH antennas involved. ; To convert the phase values to degrees, read the integer value and multiply by 0.01. ; To apply to baseline data, use the phase difference for the two antennas. ; SIDE EFFECTS: ; Data are written to the files !DEFAULTS.EPHEMDIR+'AMPCAL.REC' and ; !DEFAULTS.EPHEMDIR+'PHZCAL.REC'. ; RESTRICTIONS: ; Provision is made for 2-m calibration factors, but the actual 2-m calibration is not ; yet set, so this routine just leaves the space for these values blank. ; MODIFICATION HISTORY: ; Written 23-Jul-1999 by Dale E. Gary ; 24-Jul-1999 DG ; Made several significant changes to the scaling of the output. All values are now ; scaled by a factor of 100. The amplitudes and errors ; 26-Nov-1999 DG ; Slight change to reflect the new form of GNP structure (introduced due ; to creation of GAINPARM segment). ; 11-Jan-2000 DG ; Eliminated hardwired directory locations. ; 06-Apr-2000 DG ; Eliminate now unneeded scaling to noise diode, since that is done in ; AMPHIT on writing the output results. Update 3C84 spectrum, and ; generally fix problems. ; 23-Jul-2000 DG ; Changed to use the function TLSNOW (no functional change). ; 01-Apr-2001 DG ; Looks like including RERR in formal errors is the wrong thing to do, ; so eliminate it by setting it to zero. ; 06-Apr-2001 DG ; Add /S273 switch to indicate that the calibrator source is 3C273, and ; add 3C273 flux measurements as of 4 Apr 2001. ;- pro make_apcal,s273=s273 nf = 43 nsect = 10 ; Number of sections (8 linear poln and 2 circular) amfdat = fltarr(6,nf,nsect) in = {hed:strarr(3),data:fltarr(6,nf)} ; Read in the R/L/I data file = dialog_pickfile(path=!defaults.workdir,filter='*.amf',title='Select the R/L amphit file.') if (file eq '') then return openr,lun,/get_lun,file readf,lun,nf1,nb1,ns1 while (not EOF(lun)) do begin readf,lun,in if (strpos(in.hed[0],'12 baseline') ge 0) then begin CASE strmid(in.hed[0],strlen(in.hed[0])-1,1) OF 'R': amfdat[*,*,0] = in.data 'L': amfdat[*,*,1] = in.data ELSE: amfdat[*,*,2] = in.data ENDCASE endif endwhile free_lun,lun ; Read in the 7-baseline I data file = dialog_pickfile(path=!defaults.workdir,filter='*.adj',title='Select the I amphit file.') if (file eq '') then return openr,lun,/get_lun,file readf,lun,nf2,nb2,ns2 for i = 3, 9 do begin readf,lun,in amfdat(*,*,i) = in.data endfor free_lun,lun ; Read frequency (assumed to be the same for all of the baselines) frq = amfdat(0,*,0) ; Convert bad values to NaN bad = where(amfdat eq -99.9,nbad) if (nbad ne 0) then amfdat[bad] = !values.f_nan !p.multi=[0,2,3,0,0] window,0,xsiz=300,ysiz=450 ; List of indexes into frequency array, indicating the frequencies ; used for each of the multi-scan I scans. Harmonics must only ; appear once in this array, so although 3.6 GHz (index 10) actually ; occurs in all 5 scans, it is listed only in the first. That is ; the scan that is used to normalize all of the phases to 3.6 GHz, ; so it is the appropriate one to contain it. iflist =[[0,5,10,15,20,25,30,35,40], $ [1,6,11,16,21,26,31,36,41], $ [2,7,12,17,22,27,32,37,42], $ [3,8,13,18,23,28,33,38,-1], $ [4,9,14,19,24,29,34,39,-1]] nscn = 5 poff = fltarr(nscn) pslp = fltarr(nscn) ; Loop over multi-scan I scan, fitting a straight line to the phase ; differences between the multi-scan phases and the RLI I scan phases. ; This removes residual phase differences, ensuring that the R, L, and ; I phase calibration are mutually consistent. for i = 0, nscn-1 do begin ; Use only the first 7 frequencies to fit the straight line ; (to avoid bad fits due to noisy data). ifuse = iflist[0:4,i] f = frq[ifuse] p = (amfdat[3,ifuse,3]-amfdat[3,ifuse,2] + 3600) mod 360 ; oplot,f,p,psym=syms[i] calout = amfdat[*,ifuse,3:*] calout[3,*,0] = p transfr = {dx:0,ix:-2} phzfit,transfr,ap,calout,0,phz,res,sigma oplot,f,res[0]+res[1]*f poff[i] = res[0] pslp[i] = res[1] endfor po = median(poff) ps = median(pslp) ; Correct RLI phases by *adding* median slope and offset for i = 0, 2 do begin amfdat[3,*,i] = amfdat[3,*,i] + (po[0]+ps[0]*frq) endfor ; Correct individual multi-scan I phases by subtracting differences for i = 0, nscn-1 do begin ; Correct the phases in this scan for the fitted line, using ; all of the frequencies in IFLIST, but removing -1 first ilst=where(iflist[*,i] ne -1) ifuse = (iflist[ilst,i]) amfdat[3,ifuse,3] = amfdat[3,ifuse,3] - ((poff[i]-po[0]) + (pslp[i]-ps[0])*frq[ifuse]) endfor p = (amfdat[3,*,3]-amfdat[3,*,2] + 3600) mod 360 plot,frq,p,psym=1,yran=[-180,180] oplot,frq,p-360,psym=1 !p.multi=0 ans = dialog_message('Does this result look correct? [No => exit]: ',/question) if (ans eq 'No') then return ; Eliminate now extraneous I section from RLI scan amfdat[*,*,2:8]=amfdat[*,*,3:9] nsect = 9 ; Normalize the amplitudes to 3C84 flux spectrum ;obsolete ; We need to measure the 3C84 spectrum relatively often, but for now we will use a ;comment ; flat spectrum, based on the plot from the University of Michigan (see the URL: ; ; http://www.astro.lsa.umich.edu/obs/radiotel/umrao.html). From this plot it appears ; ; that using 21 Jy across the entire spectrum will be within 5% of the true flux. if (keyword_set(s273)) then begin ; 3C273 flux spectrum based on 04 Apr 2001 measurements flux = [63.1969,60.5378,58.0073,55.6016,53.3166,49.0930,$ 47.1468,45.3059,41.9254,40.3785,38.9225,37.5538,$ 36.2690,35.0647,33.9377,32.8846,31.9021,30.9873,$ 30.1369,29.3478,28.6172,27.3193,26.7464,25.7385,$ 24.8968,24.2006,23.6299,23.1657,22.7900,22.4854,$ 22.2356,21.9300,21.6642,21.3940,21.0796,20.6853,$ 20.0000,20.0000,20.0000,20.0000,20.0000,20.0000,20.0000] endif else begin ; 3c84 flux spectrum based on 14 Jan 2000 measurements flux = [21.6311,23.9239,25.5927,26.5496,26.9153,26.5253,$ 26.0855,25.6211,24.6414,24.1305,23.6082,23.0765,$ 22.5382,22.0049,21.4925,21.0173,20.5957,20.2446,$ 19.9722,19.7685,19.6220,19.5219,19.4199,19.4046,$ 19.4370,19.4802,19.5104,19.5277,19.5347,19.5334,$ 19.5255,19.5048,19.4775,19.4474,19.4159,19.3832,$ 19.3376,19.2897,19.2388,19.1841,19.1241,19.0567,$ 18.9785] endelse ; Adjust amplitudes to 3C84 flux, to form units/SFU calibration ; factors. Also adjust errors to the same units. fac = 20000./flux ; 9 factors to convert to SFU at ND_eff=500 for j = 0, nsect-1 do begin ; Invert to get units/SFU good = where(amfdat[1,*,j] ne -99.9,ngood) if (ngood ne 0) then amfdat[1,good,j] = amfdat[1,good,j]*fac[j] good = where(amfdat[2,*,j] ne -99.9,ngood) if (ngood ne 0) then amfdat[2,good,j] = amfdat[2,good,j]*fac[j] endfor ; At this stage, all of the amplitudes and phases are properly normalized, etc., so ; all we need to do is combine into feed-based values and write it all out to a file. ; First we need to average the 3 measurements of the quantity I1C/I2C. ratio = fltarr(nf) rerr = fltarr(nf) n = fltarr(nf) for i = 0, 2 do begin ; See which entries have both 1C and 2C okay good = where(finite(amfdat[1,*,3+i]) and finite(amfdat[1,*,6+i]),ngood) if (ngood ne 0) then begin ; Increment the number of good points found, and add them to ratio n[good] = n[good] + 1 ratio[good] = ratio[good] + reform(amfdat[1,good,3+i]/amfdat[1,good,6+i]) rerr[good] = rerr[good] + (amfdat[2,good,3+i]/amfdat[1,good,3+i])^2 $ + (amfdat[2,good,6+i]/amfdat[1,good,6+i])^2 endif endfor bad = where(ratio eq 0,nbad) if (nbad gt 0) then ratio[bad] = !values.f_nan ; Square of relative error in the ratio rerr = rerr/(ratio^2) ; Average ratio ratio = ratio/(n>1) ; Rename some arrays for ease of notation R12 = reform(amfdat[1,*,0]) & PR12 = reform(amfdat[3,*,0]) L12 = reform(amfdat[1,*,1]) & PL12 = reform(amfdat[3,*,1]) I12 = reform(amfdat[1,*,2]) & PI12 = reform(amfdat[3,*,2]) I14 = reform(amfdat[1,*,3]) & PI14 = reform(amfdat[3,*,3]) I15 = reform(amfdat[1,*,4]) & PI15 = reform(amfdat[3,*,4]) I16 = reform(amfdat[1,*,5]) & PI16 = reform(amfdat[3,*,5]) I24 = reform(amfdat[1,*,6]) & PI24 = reform(amfdat[3,*,6]) I25 = reform(amfdat[1,*,7]) & PI25 = reform(amfdat[3,*,7]) I26 = reform(amfdat[1,*,8]) & PI26 = reform(amfdat[3,*,8]) ; Corresponding squares of relative errors dR12 = (reform(amfdat[2,*,0]/amfdat[1,*,0]))^2 dL12 = (reform(amfdat[2,*,1]/amfdat[1,*,1]))^2 dI12 = (reform(amfdat[2,*,2]/amfdat[1,*,2]))^2 dI14 = (reform(amfdat[2,*,3]/amfdat[1,*,3]))^2 dI15 = (reform(amfdat[2,*,4]/amfdat[1,*,4]))^2 dI16 = (reform(amfdat[2,*,5]/amfdat[1,*,5]))^2 dI24 = (reform(amfdat[2,*,6]/amfdat[1,*,6]))^2 dI25 = (reform(amfdat[2,*,7]/amfdat[1,*,7]))^2 dI26 = (reform(amfdat[2,*,8]/amfdat[1,*,8]))^2 ; Calculate the feed-based amplitudes A2R = R12/sqrt(ratio*I12) A1L = L12*sqrt(ratio/I12) A1I = sqrt(I12*ratio) A2I = sqrt(I12/ratio) A4I = sqrt(I14*I24/I12) A5I = sqrt(I15*I25/I12) A6I = sqrt(I16*I26/I12) ifrq = reform(nint(frq/0.2)-5) ; Converts frequency array to array of indexes amp = fltarr(10,86) + !values.f_nan ; Generate an array of NaN values ; Overwrite NaN with these feed-based amplitude factors, scaled by 100 amp[0,ifrq] = A1I amp[1,ifrq] = A1L amp[2,ifrq] = A2R amp[3,ifrq] = A2I amp[4,ifrq] = A4I amp[5,ifrq] = A5I amp[6,ifrq] = A6I amp = amp*100 ; Including RERR appears to be the wrong thing to do, since these errors ; divide out exactly when the corresponding baseline combinations are constructed. ; Rather than change the code to eliminate RERR, I set it to zero for now... rerr = 0 ; Calculate the corresponding errors dA2R = 0.5*A2R*sqrt(4*dR12 + rerr + dI12) dA1L = 0.5*A1L*sqrt(4*dL12 + rerr + dI12) dA1I = 0.5*A1I*sqrt(dI12 + rerr) dA2I = 0.5*A2I*sqrt(dI12 + rerr) dA4I = 0.5*A4I*sqrt(dI14 + dI24 + dI12) dA5I = 0.5*A5I*sqrt(dI15 + dI25 + dI12) dA6I = 0.5*A6I*sqrt(dI16 + dI26 + dI12) damp = fltarr(10,86) + !values.f_nan ; Generate an array of NaN values ; Overwrite NaN with these feed-based amplitude factors, scaled by 100 damp[0,ifrq] = dA1I damp[1,ifrq] = dA1L damp[2,ifrq] = dA2R damp[3,ifrq] = dA2I damp[4,ifrq] = dA4I damp[5,ifrq] = dA5I damp[6,ifrq] = dA6I damp = damp*100 ; Calculate the feed-based phases P2R = (PR12 + 18000.) mod 360. P1L = (PI12 - PL12 + 18000.) mod 360. P1I = replicate(0.,nf) P2I = (PI12 + 18000.) mod 360. P4I = (PI14 + 18000.) mod 360. P5I = (PI15 + 18000.) mod 360. P6I = (PI16 + 18000.) mod 360. phz = fltarr(10,86) + !values.f_nan ; Generate an array of NaN values ; Overwrite NaN with these feed-based phase factors, scaled by 100 phz[0,ifrq] = P1I phz[1,ifrq] = P1L phz[2,ifrq] = P2R phz[3,ifrq] = P2I phz[4,ifrq] = P4I phz[5,ifrq] = P5I phz[6,ifrq] = P6I phz = phz*100 ; Get an appropriate TIME/LABEL structure containing the system ; date/time. tl_struct = tlsnow(!SEGM.AMPCAL) ; (For now, just ask for date. This will be in the AMF file soon) tls1 = get_tl_struct() read,'Enter year, day of year, and start HHMM of I file [YYYY DDD HHMM]: ',year,doy,hhmm tls1.year = year tls1.day = doy tls1.yrday = string(year,doy,format='(I4,".",I3.3)') tls1.msec = hms2sec(hhmm*100L)*1000L tls1.timstr = string(fix(hhmm/100),hhmm mod 100,format='(" ",I2.2,":",I2.2,":00.000")') read,'Enter year, day of year, and start HHMM of R/L file [YYYY DDD HHMM]: ',year,doy,hhmm tls2 = tls1 tls2.year = year tls2.day = doy tls2.yrday = string(year,doy,format='(I4,".",I3.3)') tls2.msec = hms2sec(hhmm*100L)*1000L tls2.timstr = string(fix(hhmm/100),hhmm mod 100,format='(" ",I2.2,":",I2.2,":00.000")') ; Convert tl_struct to binary data tldata = tl_encode(tl_struct) tlsi = tl_encode(tls1) ; Mostly empty TLS structure corresponding to the I data tlsrl = tl_encode(tls2) ; Mostly empty TLS structure corresponding to the R/L data ; Declare storage for entire amplitude record (2 records) data = intarr(1024) data2 = intarr(1024) ; Fill in record data[0:11] = tldata ; Time/Label field data[12:23] = tlsi ; Time/Label field of 27-m data epoch data[24:35] = tlsrl ; Time/Label field of 2-m data epoch ; For the second record, update the NRS field tl_struct.nrs = 2 tldata = tl_encode(tl_struct) data2[0:11] = tldata ; Time/Label field ; Loop over frequency. Replace all bad (NaN) data with 'FFFF'x for i = 0, 48 do begin bad = where(finite(amp[*,i]) eq 0,nbad) ampdat = uint(reform(amp[*,i])) if (nbad ne 0) then ampdat[bad] = 'FFFF'x bad = where(finite(damp[*,i]) eq 0,nbad) rmsdat = uint(reform(damp[*,i])) if (nbad ne 0) then rmsdat[bad] = 'FFFF'x ; AMPDAT now contains the 7 scaled calibration factors for feeds 1I, 1L, 2R, 2I, 4I, 5I, and 6I. ; RMSDAT now contains the 7 scaled rms values corresponding to these feeds. ; To get back the actual calibration factors in units/SFU, read the integer data, ; multiply by 0.01, and apply the factors to the baseline data by multiplying by ; the amp factors for BOTH antennas. wrdoff = 36 + i*20 data[wrdoff] = ampdat data[wrdoff+10] = rmsdat endfor for i = 49, 85 do begin bad = where(finite(amp[*,i]) eq 0,nbad) ampdat = uint(reform(amp[*,i])) if (nbad ne 0) then ampdat[bad] = 'FFFF'x bad = where(finite(damp[*,i]) eq 0,nbad) rmsdat = uint(reform(damp[*,i])) if (nbad ne 0) then rmsdat[bad] = 'FFFF'x ; AMPDAT now contains the 7 scaled calibration factors for feeds 1I, 1L, 2R, 2I, 4I, 5I, and 6I. ; RMSDAT now contains the 7 scaled rms values corresponding to these feeds. ; To get back the actual calibration factors in units/SFU, read the integer data, ; multiply by 0.01, and apply the factors to the baseline data by multiplying by ; the amp factors for BOTH antennas. wrdoff = 12 + (i-49)*20 data2[wrdoff] = ampdat data2[wrdoff+10] = rmsdat endfor cd,!defaults.ephemdir,current=cwd ; Open the output file to write the record. This will overwrite any ; existing file, or create a new one if none exists. openw,lun,/get_lun,'ampcal.rec',bufsiz=2048 ; Create an associated variable for the file b = assoc(lun,intarr(1024)) ; Write the data b(0) = data b(1) = data2 ; Close the file and free the unit number free_lun,lun ; The amplitude file has been written. Now write the phase file. tl_struct = tlsnow(!SEGM.PHZCAL) ; Convert tl_struct to binary data tldata = tl_encode(tl_struct) tlsi = tl_encode(tls1) ; Mostly empty TLS structure corresponding to the I data tlsrl = tl_encode(tls2) ; Mostly empty TLS structure corresponding to the R/L data ; Declare storage for entire amplitude record (2 records) data = intarr(1024) ; Fill in record data[0:11] = tldata ; Time/Label field data[12:23] = tlsi ; Time/Label field of 27-m data epoch data[24:35] = tlsrl ; Time/Label field of 2-m data epoch ; Loop over frequency. Replace all bad (NaN) data with 'FFFF'x for i = 0, 85 do begin bad = where(finite(phz[*,i]) eq 0,nbad) phzdat = uint(reform(phz[*,i])) if (nbad ne 0) then phzdat[bad] = 'FFFF'x ; PHZDAT now contains the 7 scaled phases for feeds 1I, 1L, 2R, 2I, 4I, 5I, and 6I. ; To get back the actual phases in degrees, read the integer data, multiply by 0.01, and ; adjust baseline phase by the difference of phases on the two antennas. wrdoff = 36 + i*10 data[wrdoff] = phzdat endfor ; Open the output file to write the record. This will overwrite any ; existing file, or create a new one if none exists. openw,lun,/get_lun,'phzcal.rec',bufsiz=2048 ; Create an associated variable for the file b = assoc(lun,intarr(1024)) ; Write the data b(0) = data ; Close the file and free the unit number free_lun,lun cd,cwd return end