Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.stsci.edu/~robberto/Main/Software/IDL4pipeline/ramp_basiclinfit.pro Дата изменения: Thu Feb 27 21:56:50 2014 Дата индексирования: Sun Mar 2 09:52:26 2014 Кодировка: IBM-866

; docformat='rst'
;+
; This function performs a basic linear fit, returns intercept and slope for each pixel
;-

;+
; Basic linear fit, returns intercept and slope for each pixel
;
; :categories:
; JWST Data Processing - NIRCam Cryo tests
; :examples:
; IDL> ramp_BasicLinFit, data0, fit_params
; :Params:
; data0: in, required, type="fltarr(2048,2048,Ngroup)"
; ramp datacube to be fit
; fit_params: out, required, type="fltarr(2048,2048,2)"
; intercept => [2048,2048,0]
; slope => [2048,2048,1]
; sigma_intercept => [2048,2048,2]
; sigma_slope => [2048,2048,3]
;
; :Author:
; - Based on clever code exploiting IDL matrix operators from R. Arendt (GSFC)
; - Adapted by M. Robberto, 25 Apr. 2013
; - Added errors on intercept and slope, MR. 25 Feb. 2014
;-
PRO RAMP_BasicLinFit, data0, fit_params
st=systime(1)
print,'begin ramp fitting: ', 0.000,' sec ',memory(/high)/1e9,' GB'
;
;remove linear trend per pixel;
;first get the size
s = size(data0)
nn = lindgen(s[3]) ;an array of 108 intgers, 0 to 107
;
;then reform data0 to arrange the reads in a time series.
;The dimension of data0 is now (4194304,108), i.e. pixels values are now ordered in
;a 2-d matrix with 108 rows and 4M columns
;
data0 = reform(data0,s[1]*s[2],s[3],/over)
;
; Here starts the ramp fitting procedure, elegantly done exploiting the IDL capabilities of dealing with matrixes.
;
sxy = [nn##data0[0:1999,*],nn##data0[2000:*,*]] ; nn##data
;
; For obscure IDL reasons, this calculation is faster when done in two parts.
;
;sxy is 4M vector, the result of a [88] rows x [4M,88] matrix multiplication.
;Each element, for each pixel, is the sum of the (i x y_i) products, where i=1..88; the other 3 terms are then obviousтАж
;
sx = total(nn)
sxx = total(nn^2)
sy =[(nn*0+1)##data0[0:1999,*],(nn*0+1)##data0[2000:*,*]] ; total(data,2) OK
ymean = rebin(sy,s[1]*s[2],s[3])/s[3] ;datacube
ssyy =[(nn*0+1)##(data0[0:1999,*]-ymean[0:1999,*])^2,(nn*0+1)##(data0[2000:*,*]-ymean[2000:*,*])^2] ; total(data,2) ok
;
;the slope and intercept can now be calculated using the 4 sums (equations for the linear regression can be found in any book on data analysis, see also a WFC3 ISR by Robberto for derivation and discussion)
;
a = (sy*sxx-sxy*sx)/(s[3]*sxx-sx^2) ;intercept
b = (s[3]*sxy-sx*sy)/(s[3]*sxx-sx^2) ;slope
;
; m is the slope and b is the intercept, of each pixel. These are still 1d vectors each 4M values
;

;restore the sizes at 2048x2048
data0 = reform(data0,s[1],s[2],s[3],/over)
a = reform(a,s[1],s[2],/over)
b = reform(b,s[1],s[2],/over)

;Errors
;in the case you care, here are the residual vs. the fit, (this is pretty good IDL programming, let me say....)
; model = rebin(a,2048,2048,s[3]) + rebin(b,2048,2048,s[3])*rebin(reform(indgen(s[3]),1,1,s[3]),2048,2048,s[3])
; fit_residuals = data0-model
SSxx = Sxx - Sx^2/s[3]
SSyy = reform(ssyy,2048,2048,/over);;(residuals^2,3)
SSxy = rebin(Sxy-Sx*Sy/s[3],2048,2048)
SSS = SQRT((SSyy-b*SSxy)/(s[3]-2))

sigma_a = SSS*SQRT(1./s[3]+(sx/s[3])^2/SSxx) ;error on the intercept
sigma_b = SQRT(( SSYY/SSxx-b^2)/(s[3]-2) )
;
fit_params = [[[a]],[[b]],[[sigma_a]],[[sigma_b]]] ;intercept, slope, sigma_intercept, sigma_slope
print,' end ramp fitting: ',systime(1)-st,' sec ',memory(/high)/1e9,' GB'
END