;+ ; NAME: ; REFINE_OFF ; PURPOSE: ; Routine containing much of the code of GCFIT, to refine the ; channel offset values for a given model of nonlinearity. ; CATEGORY: ; OVRO APC DATA ANALYSIS ; CALLING SEQUENCE: ; newoff = refine_off(gcdata,gcflag,nlparm,oldoff) ; INPUTS: ; gcdata the GCAL data array, of size (2,7,5,86), where ; the first index is ND off/on, second index is the ; attenuation state (0,5,10,15,20,35,20*DB), the ; third index is the total power channel, and the ; last index is the frequency ; gcflag a corresponding array containing 0 where the data are ; bad and 1 where the data are good. ; nlparm a 5-element array of nonlinearity parameters, one for ; each data channel. ; oldoff a 5-element array of offsets, one for each channel ; OPTIONAL (KEYWORD) INPUT PARAMETERS: ; ROUTINES CALLED: ; OUTPUTS: ; gcparm a 5-element array of revised offsets ; COMMENTS: ; SIDE EFFECTS: ; RESTRICTIONS: ; MODIFICATION HISTORY: ; Written 23-Jan-1998 by Dale E. Gary ;- pro flagpair,gccopy,c0,c1 ; Break out arrays for the two columns gc0 = gccopy(*,c0,*,*) gc1 = gccopy(*,c1,*,*) ; Find which entries are bad in column c1 and flag column c0 ; accordingly bad = where(gc1 EQ 0) IF (bad(0) NE -1) THEN gc0(bad) = 0 ; Find which entries are bad in column c0 and flag column c1 ; accordingly bad = where(gc0 EQ 0) IF (bad(0) NE -1) THEN gc1(bad) = 0 ; Put the flagged values back gccopy(*,c0,*,*) = gc0 gccopy(*,c1,*,*) = gc1 return end ;------------------------------------------------------------------------- function refine_off,gcdata,gcflag,nlparm,oldoff gcparm = fltarr(7,5) dbnom = [1., 1., .31623, .1, .01, .01] ; Make copy of data to work with gcd = gcdata gcdflg = gcflag ; Subtract old offsets and flag data where inconsistent for i = 0, 4 do begin gcd(*,*,i,*) = gcd(*,*,i,*) - oldoff(i) endfor bad = where(gcd LT 0.001) ; Flag all negative values (for taking log) good = where(gcd GT 0) IF (bad(0) NE -1) THEN BEGIN gcd(bad) = -1 gcdflg(bad) = 0 ENDIF IF (good(0) NE -1) THEN BEGIN gcdflg(good) = 1 ENDIF ; Apply nonlinearity parameters for each channel for i = 0, 4 do begin gcd(*,*,i,*) = gcd(*,*,i,*)*(1 + nlparm(i)*gcd(*,*,i,*)/2047.) endfor ; Save a linearized copy gclin = gcd ; Take logarithm of all good data, and set bad data to zero (so ; that bad data do not contribute to sums origood = where(gcdflg EQ 1) gcd(origood) = alog(gcd(origood)) origbad = where(gcdflg EQ 0) IF (origbad(0) NE -1) THEN gcd(origbad) = 0 ; Make another copy to work with gccopy = gcd gcfcpy = gcdflg ; For 5 DB calculation, make sure missing values are flagged ; in pairs in columns 0,1 (0:5 DB) and columns 2,3 (10:15 DB) flagpair,gccopy,0,1 flagpair,gcfcpy,0,1 flagpair,gccopy,2,3 flagpair,gcfcpy,2,3 ; Accumulate sums. Inner TOTAL sums over frequency, second TOTAL ; sums over ND on/off, and outer TOTAL sums over the attn settings ; listed. Resultant sums are 5 element arrays, one for each total ; power channel sum1 = reform(total(total(total(gccopy(*,[0,2],*,*),4),1),1)) ; 0 DB & 10 DB n1 = reform(total(total(total(gcfcpy(*,[0,2],*,*),4),1),1)) sum2 = reform(total(total(total(gccopy(*,[1,3],*,*),4),1),1)) ; 5 DB & 15 DB n2 = reform(total(total(total(gcfcpy(*,[1,3],*,*),4),1),1)) FOR i = 0, 4 DO BEGIN IF (n1(i) NE n2(i)) THEN print,'Flagging in pairs is not complete!' ENDFOR attn5 = exp(sum2/n2-sum1/n1) gcparm(2,*) = attn5/dbnom(2) ; For 10 DB calculation, make sure missing values are flagged ; in pairs in columns 0,2 (0:10 DB) and columns 1,3 (5:15 DB) gccopy = gcd gcfcpy = gcdflg flagpair,gccopy,0,2 flagpair,gcfcpy,0,2 flagpair,gccopy,1,3 flagpair,gcfcpy,1,3 ; Accumulate sums. Inner TOTAL sums over frequency, second TOTAL ; sums over ND on/off, and outer TOTAL sums over the attn settings ; listed. Resultant sums are 5 element arrays, one for each total ; power channel sum1 = reform(total(total(total(gccopy(*,[0,1],*,*),4),1),1)) ; 0 DB & 5 DB n1 = reform(total(total(total(gcfcpy(*,[0,1],*,*),4),1),1)) sum2 = reform(total(total(total(gccopy(*,[2,3],*,*),4),1),1)) ; 10 DB & 15 DB n2 = reform(total(total(total(gcfcpy(*,[2,3],*,*),4),1),1)) FOR i = 0, 4 DO BEGIN IF (n1(i) NE n2(i)) THEN print,'Flagging in pairs is not complete!' ENDFOR attn10 = exp(sum2/n2-sum1/n1) gcparm(3,*) = attn10/dbnom(3) ; For 20 DB calculation, make sure missing values are flagged ; in pairs in columns 2,4 (10:20 DB) gccopy = gcd gcfcpy = gcdflg flagpair,gccopy,2,4 flagpair,gcfcpy,2,4 ; Accumulate sums. Inner TOTAL sums over frequency, second TOTAL ; sums over ND on/off, and outer TOTAL sums over the attn settings ; listed. Resultant sums are 5 element arrays, one for each total ; power channel sum1 = reform(total(total(total(gccopy(*,2,*,*),4),1),1)) ; 10 DB n1 = reform(total(total(total(gcfcpy(*,2,*,*),4),1),1)) sum2 = reform(total(total(total(gccopy(*,4,*,*),4),1),1)) ; 20 DB n2 = reform(total(total(total(gcfcpy(*,4,*,*),4),1),1)) FOR i = 0, 4 DO BEGIN IF (n1(i) NE n2(i)) THEN print,'Flagging in pairs is not complete!' ENDFOR attn20 = exp(sum2/n2 - sum1/n1)*attn10 gcparm(4,*) = (attn20/dbnom(4)) ; From the original linearized data, separate the data for the first three ; attenuation settings. The resulting arrays have size (2,5,86) gc0db = reform(gclin(*,0,*,*)) ; 0 DB gc5db = reform(gclin(*,1,*,*)) ; 5 DB gc10db = reform(gclin(*,2,*,*)) ; 10 DB ; Determine receiver and noise diode values for each frequency. Use ; order of preference 0 DB, 5 DB, 10 DB. ; First correct the 5 DB and 10 DB columns to 0 DB attenuation. FOR i = 0, 4 DO BEGIN ; Loop over the 5 total power channels gc5db(*,i,*) = gc5db(*,i,*)/attn5(i) ; 5 DB gc10db(*,i,*) = gc10db(*,i,*)/attn10(i) ; 10 DB ENDFOR ; Determine missing values in 5 DB column and replace with 10 DB values bad = where(gc5db LE 0) if (bad(0) NE -1) THEN gc5db(bad) = gc10db(bad) ; Determine missing values in 0 DB column and replace with 5 DB values bad = where(gc0db LE 0) if (bad(0) NE -1) THEN gc0db(bad) = gc5db(bad) ; GC0DB now contains the best data possible, with missing data filled ; in from the next higher attenuation setting, suitably adjusted to ; 0 DB. Data values for a given frequency are now missing only if ; they are missing in all three attenuation settings, 0, 5 and 10 DB. ; Now get receiver values for each frequency. GRCVR has size (5,86). grcvr = reform(gc0db(0,*,*)) ; Get the noise diode increments by subtracting the GRCVR values. ; Values with ND entries of zero are automatically zeroed. Values ; with RCVR entries of zero are zeroed manually. gnd = (reform(gc0db(1,*,*)) - grcvr)>0 bad = where(grcvr LE 0) IF (bad(0) NE -1) THEN gnd(bad) = 0 attn15 = attn5*attn10 attn35 = attn15*attn20 a5 = attn5#replicate(1,86) a10 = attn10#replicate(1,86) a15 = attn15#replicate(1,86) a20 = attn20#replicate(1,86) a35 = attn35#replicate(1,86) ; Now that we have a best guess at the 35 DB attenuation, and the ; receiver and ND increments at each frequency, we can estimate ; the amount of true signal we subtracted when determining offsets. ; Subtract this small amount from the original 35 DB input data, and ; determine new offset values. gcoff = reform(gcdata(*,5,*,*)) gcoflg = reform(gcflag(*,5,*,*)) offsig = grcvr*a35 onsig = (grcvr+gnd)*a35 gcoff(0,*,*) = gcoff(0,*,*) - offsig gcoff(1,*,*) = gcoff(1,*,*) - onsig ; Flag data where inconsistent bad = where(gcoff LT 0) ; Flag all negative values good = where(gcoff GT 0) IF (bad(0) NE -1) THEN BEGIN gcoff(bad) = 0 gcoflg(bad) = 0 ENDIF IF (good(0) NE -1) THEN BEGIN gcoflg(good) = 1 ENDIF ; Finally, determine the new offset values and return. ; The inner TOTAL function sums over ND on/off, the outer one ; sums over frequency. newoff = total(total(gcoff,1),2)/((total(total(gcoflg,1),2))>1) return,newoff end