analyze_dark_ramp.pro (JWST Pipeline routines in IDL)

;-
; docformat = 'rst'
;+
; Wrapper to correct NIRCam dark ramps for 1/f noise and derive dark map and RON
;-    
   
;+
; Wrapper to correct NIRCam dark ramps for 1/f noise and derive dark map and RON
; :Categories: 
;       JWST Data Processing - NIRCam Cryo tests
; :Examples:
;       IDL> analyze_dark_ramp
; :uses:
;       Common blocks::
;       
;          COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask;
;       
; :History: created by M. Robberto,  25 Apr. 201
;-
pro analyze_dark_ramp
COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
st=systime(1)
;on robberto2
DIR_OUT = '~/FUNCTIONAL/JWST/NIRCAM/Studies/2013/2013.NIRCAM.Noise/Results/'
;on witserv
;DIR_OUT = '/witserv/data15/data_analysis_results/'
;
;set detector parameters
SET_DETECTOR_PARAMS
;
;get ready to find files
 DATA_PATH= '/Volumes/witserv'
 DATA_PATH = DATA_PATH+'/witserv/data14/nrc/nrc_cryo2' 
 OBS_ID = 'NRC_2D6-DARK'
 DETECTOR = '486' 
 FIND_FILES,DATA_PATH,OBS_ID,DETECTOR,FILES
;
;DATA_PATH = '~/FUNCTIONAL/JWST/2011/RauscherReadOut/NIRSPEC_data'
;files = FILE_SEARCH(DATA_PATH+'/*.fits',count=nfiles)
;
;read the first file and process
 N_FILES = N_ELEMENTS(FILES)
 FOR jf = 0, N_files-1 do begin
   print,'begin reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
   data0 = float(readfits(files[jf],/silent))*gain
   print,'  end reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
   ;
   ;
   ;>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
   ;EXTRACT THE DIR AND FILE STRINGS AND SET THE DIRECTORY IN OUTPUT
   ;first the position of the bar in the filename
   bar = STREGEX(files[jf],'/')      
   ;
   ;then the last character, before the .fits extension
   dfits = STREGEX(files[jf],'.fits')
   ;
   ;now get the directory name
   filedirname  = strmid(files[jf],0,bar-1)
   ;
   ;and the filebasename
   filebasename = strmid(files[jf],bar+1,dfits-bar-1)
   ;
   ;put dire and filebasename in a single variable for the output
   fileinfo = [DIR_OUT,filedirname, filebasename]
   ;done with the dir and file strings
   ;<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
;

;run the analysis
STOP
   RAMP_COMMON_LEVEL, data0
   RAMP_FIXSECTORS, data0
   correct_Evenodd_ramp, data0
   ramp_OPTIMAL_FILTER, data0
   ;make dark datacube
   RAMP_BasicLinFit, data0, darkcube
   RAMP_PlotBasicLinFit, darkcube, dark_averages, fileinfo   
   RAMP_GETRON, data0, darkcube, RON
   RAMP_PLOT_4histo,RON, RON_AVERAGES, MASK, fileinfo 
   STOP
   PRINT_RESULTS,fileinfo,darkcube,ron, dark_averages,ron_averages
 ENDFOR  
END

;;-------------------------------------------------------------
;;+
;; NAME:
;;       DETECTOR_PARAMS
;; PURPOSE:
;;       Specifies main detector parameters
;; CATEGORY: 
;;       JWST Data Processing
;; CALLING SEQUENCE:
;;       detector_params
;; INPUTS: 
;;       none
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       none
;; COMMON BLOCKS:
;;       COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
;; NOTES:
;;       The detector parameters are transmitted out through the common block
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO DETECTOR_PARAMS
;COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
;gain=2
;h =lindgen(2048,2048) 
;h1mask = h[0:511,*]
;h2mask = rotate(h[512:1023,*],5)
;h3mask = h[1024:512*3-1,*]     
;h4mask = rotate(h[512*3:*,*],5) 
;END

;;-------------------------------------------------------------
;;+
;; NAME:
;;       FIND_FILES
;; PURPOSE:
;;       Select a set of fits files to be processed, all in a given directory
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       find_files, data_path, obs_id, detector, files
;; INPUTS: 
;;       data_path: Directory where the fits files are stored
;;       obs_id: Initial String common to all filenames, e.g. 'NRC_2D6-DARK' 
;;       detector: Detector ID number, e.g. '485' (string)
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       files: string array of fits files with ames matching the obs_id and detector strings
;; COMMON BLOCKS:
;;       COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask;
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       The detector parameters are transmitted out through the common block
;;       Example:
;;         DATA_PATH = '/Volumes/witserv/witserv/data14/nrc/nrc_cryo2' 
;;         OBS_ID = 'NRC_2D6-DARK'
;;         DETECTOR = '486' 
;;         FIND_FILES,DATA_PATH,OBS_ID,DETECTOR,FILES
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 201
;;-
;;-------------------------------------------------------------
;PRO FIND_FILES,DATA_PATH,OBS_ID,DETECTOR,FILES
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
;print,'begin reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;cd, DATA_PATH
;;for test purposes
;; OBS_ID = 'NRC_2D6-DARK'
;; DETECTOR = '486' 
;files = FILE_SEARCH(OBS_ID+'*'+'/*'+DETECTOR+'*.fits',count=nfiles)
;;
;;;read datacube, assume gain=2
;;;jf = 0
;;  data0 = float(readfits(file,/silent))*gain
;;print,'  end reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;print,'  end reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;END

;;-------------------------------------------------------------
;;+
;; NAME:
;;       READ_RAMP
;; PURPOSE:
;;       Read a single fits file in multiaccum format
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       READ_RAMP, data0
;; INPUTS: 
;;       none
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       data0: ramp datacube, typical format [2048,2048,Ngroups]
;; COMMON BLOCKS:
;;       COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask;
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       This routine is used to read a single file, identified inside (editing needed); 
;;       used for debugging
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO READ_RAMP,data0
;;to test: find and read only one file
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
;st=systime(1)
;print,'begin reading ramp: ',st,' sec     ',memory(/high)/1e9,' GB'
;Cross_correlation = 0
;
;DATA_PATH= '/Volumes/witserv'
;DATA_PATH = DATA_PATH+'/witserv/data14/nrc/nrc_cryo2' 
;cd, DATA_PATH
; OBS_ID = 'NRC_2D6-DARK'
; DETECTOR = '486' 
;file = FILE_SEARCH(OBS_ID+'*'+'/*'+DETECTOR+'*.fits',count=nfiles)
;;
;;read datacube, assume gain=2
;gain=2
;jf = 0
;  data0 = float(readfits(file[jf],/silent))*gain
;print,'  end reading ramp: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;END

;;-------------------------------------------------------------
;;+
;; NAME:
;;       RAMP_FIXSECTORS
;; PURPOSE:
;;       Correct for the amplifier reset of NIRCam ASICs by matching the top/bottom reference pixes of adjacent groups
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       ramp_fixsectors, data0
;; INPUTS: 
;;       data0; ramp datacube to be corrected, assumes format [2048,2048,Ngroup]
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       data0: ramp datacube after correction, assumes format [2048,2048,Ngroups]
;; COMMON BLOCKS:
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       This routine assumes that the average of the 4 bottom rows of group N+1 is the same of the average of the 
;;       top rows of group N+1, and equalizes accordingly.  
;;       The comparison is performed after subtraction of the average level of each pixel.
;;       => Clearly, this assumption becomes rather weak/false if the groups are "real groups", i.e. if the frame is an average
;;       of 2,4,8 single frames.
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO RAMP_FIXSECTORS,data0
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;print,'begin fixing sectors: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;;
;;x  because of the amplifier reset, need to homogenize the signal ramps
;;x  h =lindgen(2048,2048) 
;;x  h1mask = h[0:511,*]
;;x  h2mask = rotate(h[512:1023,*],5)
;;x  h3mask = h[1024:512*3-1,*]     
;;x  h4mask = rotate(h[512*3:*,*],5) 
;
;;
;;Find the global averages of the bottom and top pixels, as we will work on the differences from the average
;mbottom1=MEAN(data0[0:511,0:3,*],dim=3) & mtop1=MEAN(data0[0:511,2044:2047,*],dim=3)
;mbottom2=MEAN(data0[512:1023,0:3,*],dim=3) & mtop2=MEAN(data0[512:1023,2044:2047,*],dim=3)
;mbottom3=MEAN(data0[1024:1535,0:3,*],dim=3) & mtop3=MEAN(data0[1024:1535,2044:2047,*],dim=3)
;mbottom4=MEAN(data0[1536:2047,0:3,*],dim=3) & mtop4=MEAN(data0[1536:2047,2044:2047,*],dim=3)
;;
;;perform the correction of the ramp
;for i=1,107 do begin & $
;data0[0:511,*,i]=data0[0:511,*,i]-MEAN(data0[0:511,0:3,i]-mbottom1)+MEAN(data0[0:511,2044:2047,i-1]-mtop1) & $
;data0[512:1023,*,i]=data0[512:1023,*,i]-MEAN(data0[512:1023,0:3,i]-mbottom2)+MEAN(data0[512:1023,2044:2047,i-1]-mtop2) & $
;data0[1024:1535,*,i]=data0[1024:1535,*,i]-MEAN(data0[1024:1535,0:3,i]-mbottom3)+MEAN(data0[1024:1535,2044:2047,i-1]-mtop3) & $
;data0[1536:2047,*,i]=data0[1536:2047,*,i]-MEAN(data0[1536:2047,0:3,i]-mbottom4)+MEAN(data0[1536:2047,2044:2047,i-1]-mtop4) & $
;endfor
;print,'  end fixing sectors: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;END

;;-------------------------------------------------------------
;;+
;; NAME:
;;       RAMP_COMMON_LEVEL
;; PURPOSE:
;;       Match the four sectors of the detector to a common, average DC level
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       common_level, data0
;; INPUTS: 
;;       data0; ramp datacube to be corrected, assumes format [2048,2048,Ngroup]
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       data0: ramp datacube after correction, assumes format [2048,2048,Ngroups]
;; COMMON BLOCKS:
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       This routine uses the average of the top and bottom reference pixels to establish the global supermean level of the ramp.  
;;       Then the reference pixels (top&bottom) of each sector&frame are tuned to match the global supermean.
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO RAMP_COMMON_LEVEL,data0
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;print,'begin setting common level: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;s = size(data0) 
;nn = lindgen(s[3])   ;an array of 108 intgers, 0 to 107
;;
;;REFPIXDOWN
;;refpix1down = data0[0:511,0:3,*]      ;isolate the datacubes of ref pixels
;refpix2down = data0[512:1023,0:3,*] ;& refpix2down = refpix2down[rev,*,*]
;refpix3down = data0[1024:1535,0:3,*]
;refpix4down = data0[1536:2047,0:3,*]  ;& refpix4down = refpix4down[rev,*,*]
;;REFPIXUP
;refpix1up = data0[0:511,2044:2047,*]      ;isolate the datacubes of ref pixels
;refpix2up = data0[512:1023,2044:2047,*] ;& refpix2up = refpix2up[rev,*,*]
;refpix3up = data0[1024:1535,2044:2047,*]
;refpix4up = data0[1536:2047,2044:2047,*] ; & refpix4up = refpix4up[rev,*,*]
;;
;;HERE ONE CAN CORRECT FOR THE EVEN AND ODD PIXE:S. here is a template
;;rp1 = [(refpix1down+refpix1up)/2.] & rp1=avg(rp1,1) & rp1=avg(rp1,1) ;array of 512 elements
;;rp1even=rp1[indgen(256)*2] & rp1odd=rp1[indgen(256)*2+1] 
;;rp1_DELTA = AVG(rp1even)-AVG(rp1odd)
;;data0[indgen(256)*2,*,*]-=rp1_DELTA/2. & data0[indgen(256)*2+1,*,*]+=rp1_DELTA/2.
;;
;;rp2 = [(refpix2down+refpix2up)/2.] & rp2=avg(rp2,1) & rp2=avg(rp2,1) ;array of 512 elements
;;rp2even=rp2[indgen(256)*2+512] & rp2odd=rp2[indgen(256)*2+1+512]
;;rp2_DELTA = AVG(rp2odd)-AVG(rp2even)
;;data0[indgen(256)*2+1+512,*,*]+=rp2_DELTA/2. & data0[indgen(256)*2+512,*,*]+=rp2_DELTA/2.
;;etc...
;;stop
;;
;;ALLTOGETHER
;alltogether = [data0[*,0:3,*],data0[*,2044:2047,*]]
;;
;;SUPERMEAN
;resistant_mean,alltogether,3,supermean ;the average of the average is the DC level of the output channel
;;stop
;;CORRECT FOR THE OFFSETS
;for iread=0,s[3]-1 do begin
;;Q1
;    refpix_quad = [refpix1down[*,*,iread],refpix1up[*,*,iread]]
;    resistant_mean,refpix_quad,3,mean_quad ;the average of the average is the DC level of the output channel 
;    data0[0:511,*,iread] = data0[0:511,*,iread] - mean_quad + supermean
;;Q2
;    refpix_quad = [refpix2down[*,*,iread],refpix2up[*,*,iread]]
;    resistant_mean,refpix_quad,3,mean_quad ;the average of the average is the DC level of the output channel 
;    data0[512:1023,*,iread] = data0[512:1023,*,iread] - mean_quad + supermean
;;Q3
;    refpix_quad = [refpix3down[*,*,iread],refpix3up[*,*,iread]]
;    resistant_mean,refpix_quad,3,mean_quad ;the average of the average is the DC level of the output channel 
;    data0[1024:1535,*,iread] = data0[1024:1535,*,iread] - mean_quad + supermean
;;Q4
;    refpix_quad = [refpix4down[*,*,iread],refpix4up[*,*,iread]]
;    resistant_mean,refpix_quad,3,mean_quad ;the average of the average is the DC level of the output channel 
;    data0[1536:2047,*,iread] = data0[1536:2047,*,iread] - mean_quad + supermean
;endfor
;print,'end setting common level: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;;image data0 is now on a even ground
;end

;;-------------------------------------------------------------
;;+
;; NAME:
;;       OPTIMAL_FILTER
;; PURPOSE:
;;       Perform an optimal filtering of the vertical reference pixel to reduce 1/f noise (horizontal stripes) 
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       optimal filter, data0
;; INPUTS: 
;;       data0; ramp datacube to be corrected, assumes format [2048,2048,Ngroup]
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       data0: ramp datacube after correction, assumes format [2048,2048,Ngroups]
;; COMMON BLOCKS:
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       Calls the routing optimal_smooth_fft.pro
;; MODIFICATION HISTORY:
;;       M. Robberto, created 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO OPTIMAL_FILTER,data0
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;print,'begin optimal filtering: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;s = size(data0) 
;nn = lindgen(s[3])   ;an array of 108 intgers, 0 to 107
;;-extract the side reference pixels
;refpixleft = data0[0:3,*,*]
;refpixright = data0[2044:2047,*,*]
;;- set to zero the average of ALL pixels
;AVEdata0=avg(data0,2)
;;-set to zero the average of each reference pixel
;mrefpixleft = avg(refpixleft,2)
;mrefpixright = avg(refpixright,2)
;for igroup=0,s[3]-1 do begin & $
;  refpixleft[*,*,igroup]=refpixleft[*,*,igroup]-mrefpixleft[*,*] & $
;  refpixright[*,*,igroup]=refpixright[*,*,igroup]-mrefpixright[*,*] & $
;endfor
;refpix = (MEAN(refpixleft,DIM=1)+MEAN(refpixright,DIM=1) ) / 2.
;;
;optimal_smooth_fft, refpix, 10.6*510, Smoothed_refpix
;for i=0,s[3]-1 do begin & $
;    optimal_smooth_fft, refpix[0+(i*2048):2047+(i*2048)], 10.6*510, smrp & $ 
;    Smoothed_Refpix[i*2048:2047+(i*2048)]=smrp & $
;    wait,0.1
;endfor
;wait,1
;for ix=0,2047 do begin
;  print,ix
;  data0[ix,*,*]-=smoothed_refpix
;endfor
;print,'end optimal filtering: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;end


;;-------------------------------------------------------------
;;+
;; NAME:
;;       RAMP_BasicLinFit
;; PURPOSE:
;;       Basic linear fit, returns intercept and slope for each pixel
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       ramp_BasicLinFit, data0, fit_params
;; INPUTS: 
;;       data0; ramp datacube to be corrected, assumes format [2048,2048,Ngroup]
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       fit_params: datacube with dark information, size [2048,2048,2]
;;         intercept => [2048,2048,0]
;;         slope     => [2048,2048,1]
;; COMMON BLOCKS:
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       none
;;; MODIFICATION HISTORY:
;;       - Based on clever code exploiting matrix operators from GSFC 
;;       - Adapted by M. Robberto, 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO RAMP_BasicLinFit, data0, fit_params
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;print,'begin darkcube creation: ',systime(1)-st,' 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)
;;
;;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)
;; 
; m = (s[3]*sxy-sx*sy)/(s[3]*sxx-sx^2) ;slope
; b = (sy*sxx-sxy*sx)/(s[3]*sxx-sx^2)  ;intercept
;;
;; 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)
;m = reform(m,s[1],s[2],/over)
;b = reform(b,s[1],s[2],/over)
;fit_params=[[[b]],[[m]]]   ;intercept & slope
;print,'  end darkcube creation: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;END


;;-------------------------------------------------------------
;;+
;; NAME:
;;       RAMP_PlotBasicLinFit
;; PURPOSE:
;;       Plot the histograms of the twi linear fit parameters, returning alsi gaussian fit results
;; CATEGORY: 
;;       JWST Data Processing - NIRCam Cryo tests
;; CALLING SEQUENCE:
;;       RAMP_PlotBasicLinFit,darkcube,dark_averages
;; INPUTS: 
;;       fit_params: datacube with dark information, size [2048,2048,2]
;;         intercept => [2048,2048,0]
;;         slope     => [2048,2048,1]
;; KEYWORD PARAMETERS:
;;       none
;; OUTPUTS:
;;       gaussfit_results: datacube with gaussfit results [5,2]
;;         intercept => [*,0]
;;         slope     => [*,1]
;;Let us remember the 5 values returned by gaussfit,array,aa 
;;      aa = coefficients of the Gaussian fit: Height, mean, sigma,...
;;               aa[0]= the height of the Gaussian
;;               aa[1]= the mean
;;               aa[2]= the standard deviation
;;               aa[3]= the half-width of the 95% conf. interval of the standard
;;                     mean
;;               aa[4]= 1/(N-1)*total( (y-mean)/sigma)^2 ) = a measure of 
;;                       normality
;; COMMON BLOCKS:
;;       COMMON SYSTEM_PARAMS,st,DIR_OUT
;; NOTES:
;;       none
;;; MODIFICATION HISTORY:
;;       - Adapted by M. Robberto, 25 Apr. 2013
;;-
;;-------------------------------------------------------------
;PRO RAMP_PlotBasicLinFit,fit_params,dark_averages 
;  window,1
;  !p.multi = [0,1,2]
;  histogauss,fit_params[4:2043,4:2043,0],aa1   ;intercept HISTOGAUSS is in the IDLASTRO library at http://idlastro.gsfc.nasa.gov
;  histogauss,fit_params[4:2043,4:2043,1],aa2   ;slope
;  gaussfit_results = [[aa1],[aa2]] ;intercept, slope fit parameters
;end

;PRO PLOT_RONMAP,RON,RON_AVERAGES, MASK, file
;COMMON SYSTEM_PARAMS,st,DIR_OUT
;COMMON DETECTOR_PARAMS,gain,h,h1mask,h2mask,h3mask,h4mask
;print,'  begin RONMAP creation: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
;
;;====================================================================
;;Masks outliers with especially high or low standard deviation
;;Before moving forward with the analysis, we clean the data set. We look for outliers, with threshold set at 4 sigma
;
;;Prepare the output 2x2 window 
;; mask outliers with especially high or low stddev
;  wset,0
;  !p.multi = [0,2,2]
;  
;;
;;We now look at the four “quadrants”, or sectors of the array; their indexes have been previously labeled h1mask…h4mask. 
;;We plot the histogram of the standard deviations; the mean of the histogram is the readout noise, in electrons. 
;;The histogauss function produces a plot, fits a Gaussian profile and returns the parameters of the fit.
;; 
;  sig_data=RON
;  histogauss,sig_data[h1mask],aa1   ; HISTOGAUSS is in the IDLASTRO library at http://idlastro.gsfc.nasa.gov
;
;;Let us remember the values of aa
;;      aa = coefficients of the Gaussian fit: Height, mean, sigma
;;               aa[0]= the height of the Gaussian
;;               aa[1]= the mean
;;               aa[2]= the standard deviation
;;               aa[3]= the half-width of the 95% conf. interval of the standard
;;                     mean
;;               aa[4]= 1/(N-1)*total( (y-mean)/sigma)^2 ) = a measure of 
;;                       normality
;
;;Next we isolate the pixels that are within +/-4 sigma from the average readout noise: mask1=1 if the condition is satisfied, 
;;0 otherwise. 
;;Again, we extract from h1mask the indexes of the pixels with a discrepancy <4sigma (with respect to the average of the quadrant). 
;;Mask 1 has values 1 or 0
;  mask1 = abs(sig_data[h1mask]-aa1[1]) lt (4*aa1[2])
;
;;idl returns 1 if the LT condition is satisfied, 0 otherwise. This procedure extracts from h1mask the values with  
;;discrepancy <4sigma (with respect to the full array). Mask 1 has values 1 or 0
;
;;Repeat for the other 3 sectors. There is here a “counter-reflection” around the y axis to have the pixels oriented normally; 
;;this because h2mask and h4mask have been originally defined flipped with respect to the horizontal (fast readout) direction
;
;  histogauss,sig_data[h2mask],aa2;_extra={charsize:0.6}
;  mask2 = rotate(abs(sig_data[h2mask]-aa2[1]) lt (4*aa2[2]),5)
;  histogauss,sig_data[h3mask],aa3
;  mask3 = abs(sig_data[h3mask]-aa3[1]) lt (4*aa3[2])
;  histogauss,sig_data[h4mask],aa4
;  mask4 = rotate(abs(sig_data[h4mask]-aa4[1]) lt (4*aa4[2]),5)
;;...again the “counter-reflection”
;
;;THE RESULTS CAN BE EXTRACTED HERE
;  RON_AVERAGES=[[aa1],[aa2],[aa3],[aa4]]
;;The variable mask contains now the indexes of all “good” pixels
;  mask = [mask1,mask2,mask3,mask4]
; 
;END



;DIR_OUT = '~/FUNCTIONAL/JWST/NIRCAM/Studies/2013/2013.NIRCAM.Noise/Results/'
;filebasename = '487'
;;
;;set detector parameters
;SET_DETECTOR_PARAMS
;;
;;get ready to find files
; DATA_PATH= '/Volumes/witserv'
; DATA_PATH = DATA_PATH+'/witserv/data14/nrc/nrc_cryo2' 
; OBS_ID = 'NRC_2D6-DARK'
; DETECTOR = '486' 
; FIND_FILES,DATA_PATH,OBS_ID,DETECTOR,FILES
PRO PRINT_RESULTS,fileinfo,darkcube,ron, dark_averages,ron_averages
;
;
   writefits,fileinfo[0]+'/'+fileinfo[1]+'/'+fileinfo[2]+'_darklinfit.fits',darkcube
   writefits,fileinfo[0]+'/'+fileinfo[1]+'/'+fileinfo[2]+'_RON.fits',ron
   openw,1,fileinfo[0]+'/'+fileinfo[1]+'/'+fileinfo[2]+'_RONresults.txt'
   printf,1,RON_AVERAGES
   close,1
  openw,2,fileinfo[0]+'/'+fileinfo[1]+'/'+fileinfo[2]+'_DARKresults.txt'
   printf,2,dark_averages
  close,2 
END   

pro openeps,nomefile
on_error,2
print,'out to ps on the file:', nomefile
set_plot,'ps'
device,bits=8,filename=nomefile,/ENCAPSULATED
end
  
pro closeps
on_error,2
device,/close
set_plot,'x'
end

  
pro analyze_dark
COMMON SYSTEM_PARAMS,st,DIR_OUT
OBS_ID = 'NRC2D6-DARK_'
DETECTOR = '487' 
cd, DIR_OUT 

; OBS_ID = 'NRC_2D6-DARK'
; DETECTOR = '486' 
files = FILE_SEARCH(OBS_ID+'*_'+DETECTOR+'_*DARKresults.txt',count=nfiles)
;
;SORT THE FILES ACCORDING TO THEIR OBS_ID
NUMBER=STRARR(N_ELEMENTS(files))
for ij=0,N_ELEMENTS(files)-1 do begin & $
;  ij=0 & 
  file=files[ij] & $
  from = STRLEN(OBS_ID) & $
  to = STREGEX(file,'_'+DETECTOR) & $
  number[IJ] = strmid(file,from,to-from) & $
endfor
  number(WHERE(strlen(number) EQ 1)) =   '0'+number(WHERE(strlen(number) EQ 1))
  print,number 
  FILES = FILES(SORT(NUMBER))
;
;  
alldarks=fltarr(2,N_ELEMENTS(files))
for ij=0,N_ELEMENTS(files)-1 do begin & $
  READCOL, FILES[ij],Peaks,Means,STDevs,mean2sigmaerror,normality,FORMAT='(F,F,F,F,F)' & $
  alldarks[*,ij]=[means[1],STDevs[1]] & $
endfor
plot,alldarks[0,*],psym=3
oploterr,alldarks[0,*],alldarks[1,*]

end  
  

PRO FOURIER
;
;FOURIER ANALYSYS
;For the Fourier Analysis, also the 4 pixels nearest to each outlier
;For better cleanup, we remove the 4 pixels adjacent to each outlier. 
;This is simply done by shifting the bad pixel mask of +/-1 in all directions
  mask *= shift(mask,-1,0)*shift(mask,1,0)*shift(mask,0,-1)*shift(mask,0,1)

;and back to single window

  !p.multi = 0
  ;wset,5

;finally take data0 and set to 0 the values of the pixels flagged as bad, frame by frame.

  for i=0,s[3]-1 do data0[*,*,i] *= mask   ;apply the mask to filter out the bad data.

;The figure below shows the bottom-left 1Mx1M corner of the mask, with two sectors: 
tvscl,mask[0:1023,0:1023]


;
;
;Fill in for the masked pixels

;Uses the nmean data from the +/-12 pixels n time domain, to avoid artifacts in the power spectrum
;The pixels set to 0 in data0 will create problems when we do the Fourier Transforms. 
;The idea here is to fill in by using the mean of the +/- 12 pixels in the readout sequence (time domain)

for i=0,3 do begin & $; This fills in masked pixels with the mean of data from +/- 12 pix in the time domain. 
  ;(helps avoid artifacts in the power spectra)

;The index i labels the quadrant, so for quadrant 1 (i=0) we have [0:511,*]. The bad pixels have index h

    h = where(mask[i*512:(i+1)*512-1,*] eq 0,nh) & $
;
;for the current quadrant, we run through the entire ramp 
   for k=0,s[3]-1 do begin & $
      mn = fltarr(nh) & $
      nn = fltarr(nh) & $
;mn sums over the values of the +/- 12 pixels. Some of them will be zero, so to find the average we will need 
;to divide by the numbers of pixels different from 0. They are counted by nn, a counter (1 or 0) that ticks only 
;if the pixel is good (i.e. non 0). It is the same operation, but the NE command makes the second statement a logical one 
      for j=-12,12 do begin  & $
        mn += (data0[i*512:(i+1)*512-1,*,k])[h-j] & $
        nn += (data0[i*512:(i+1)*512-1,*,k])[h-j] ne 0 & $
      endfor & $
;find the average dividing by the number of good pixels: 
      mn /= (nn>1) & $
;reorganize the data substituting the average values at the positions of the bad pixels; the indexes are arranged to reproduce the [512,2048] sector 
      for j=0L,nh-1 do data0[(h[j] mod (s[1]/4))+(i*512),h[j]/(s[1]/4),k] = mn[j] & $
   endfor & $
endfor
;
;out of the loop on the quadrants we get a new time  
  print,'mark3: ',systime(1)-st
;
;
;Set the spatial average of each output to zero 
;
;We are moving closer to the Fourier analysis. At the moment we have that each pixel has been fitted 
;by a linear ramp and subtracted; this has facilitated finding the outliers. New we need to go spatial 
;across the frame to have time sequences that have zero mean. 

; set means to 0.0 for each output (4) for each frame (108)
;subtract to each sector its own average to have spatial average to zero
for j=0,s[3]-1 do begin & $
   for i=0,3 do begin  & $
      frame = data0[i*512:(i+1)*512-1,*,j] & $
      h = where(frame ne 0) & $
      frame[h] -= mean(frame[h]) & $
      data0[i*512:(i+1)*512-1,*,j] = frame & $
   endfor& $
endfor
;
;data0 now has zero average, quadrant by quadrant, for each frame. 
  print,'mark4: ',systime(1)-st


;5.5 Reform data 
;From [2048,2048,n] to [512,4,2048,n], 
; then reorder to [512,2048,n,4], 
;  then rearrange quadrants 1 and 3 of the array, so that the pixels are oriented according to the readout directions. 
;   then add padding rows
;Sectors 2 and 4 need to be flipped to have them in sync with quadrants 1 and 3. Start with sector 2:
; reform the data, add padding on the rows, and take FFT

;let’s start with the creation of the reverted indexes:
  rev = reverse(lindgen(s[1]/4))
;
;and prepare the datacube. First, reform data0 from [2048,2048,n] to [512,4,2048,n]
  data0 = reform(data0,s[1]/4,4,s[2],s[3],/over)
;
;then reform again to [512,2048,n,4]. Piece of cake with IDL:
  data0 = transpose(data0,[0,2,3,1]) 
;
;now use the reverted indexes to flip the second and fourth sector 
  data0[*,*,*,1] = data0[rev,*,*,1] 
  data0[*,*,*,3] = data0[rev,*,*,3]

;
;
;Add pad=12 pixels at the end of each row
pad=12
;A final extra step: we need to expand the x size of data0 to add 12 pixels cycles, 
;needed each time we move to the next line. There are the padding rows; the data look 
;better in the Fourier space.
;
;First modify data0 to a fatter datacube: 512->524, i.e. [524,2048,n,4]
 data0 = [data0,fltarr(pad,s[2],s[3],4)] 
;
;We need to put something into these extra pixels. 
;For each sector [0-3], frame [0-87] and row [0:2046] we insert the average of the last 12 pixels 
;of each row and the first 12 pixels of next rows. 
  for k=0,3 do for j=0,s[3]-1 do for i=0L,s[2]-2 do $
    data0[512:*,i,j,k] = total(data0[500:511,i,j,k])/24+total(data0[0:11,i+1,j,k])/24 ; fixed minor bug here (i+1)
;
;where we go up to 2046 to jump at 2047 at a valid row. 
;
;At this point we can reform again the data to have a [1M, 108,4] datacube. 
;We will work on the 1M array as a 1-d signal to be analyzed with Fourier. 

  data0 = reform(data0,(s[1]/4+pad)*s[2],s[3],4,/over)
  print,'mark5: ',systime(1)-st

;
;
;
;Fourier transform
;=========================================
;;Overwriting the data; 
;pay attention to the normalization factor SQRT(nr.of pixels=524*2028) 
;
  data0 = fft(data0,dim=1,/overwrite)*sqrt((s[1]/4+pad)*s[2])
  print,'mark6: ',systime(1)-st

;We can now plot the FFT for a e.g. the first sector of the first frame
plot,data0[*,0,0]

;
;
;5.8  Calculation of the Power and Cross correlation spectra

; Calculate power and cross power for each output (RAW DATA).
;..  if (option) then begin
 
;Calculate here the Power spectrum, one for each sector; it is calculated as a mean over the 108 frames, 
;since the total has a keyword 2, which means sum over the second dimension. These are 1-d   arrays of 
;512x2048 elements that we can plot
scipix_n = 512L
refpix_r = 0L
row = 512L+pad+(refpix_r+2)*(512/scipix_n)    ; LONG 716 = 512 + 12 + (4+2)*(512/16)
f = dindgen(2048L*row)/2048L/row ; frequency array ; DOUBLE    = Array[1466368]
f5 = f*1e5  ;set the highest frequency is at 100KHz (10ms/frame)
f0 = dindgen(2048L*512L)/2048L/512L ; frequency array for normal pixels only (ignoring all gaps and references);
f05 = f0*1e5



 q1_0 = fltarr(2048L*row)
 q2_0 = fltarr(2048L*row)
 q3_0 = fltarr(2048L*row)
 q4_0 = fltarr(2048L*row)

    q1_0 += total(abs(data0[*,*,0])^2,2)/s[3]
    q2_0 += total(abs(data0[*,*,1])^2,2)/s[3]
    q3_0 += total(abs(data0[*,*,2])^2,2)/s[3]
    q4_0 += total(abs(data0[*,*,3])^2,2)/s[3]

  print,'mark7: ',systime(1)-st

  !y.title = 'Power (e!u2!n) [Raw, single ended]'
  density,/ylog,f*1e5,q1_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 1: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q1_0)/8.031746),2)
  density,/ylog,f*1e5,q2_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 2: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q2_0)/8.031746),2)
  density,/ylog,f*1e5,q3_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 3: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q3_0)/8.031746),2)
  density,/ylog,f*1e5,q4_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 4: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q4_0)/8.031746),2)

  
  density,/ylog,f*1e5,q1_0,xran=[0,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 1: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q1_0)/8.031746),2)
  density,/ylog,f*1e5,q2_0,xran=[0,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 2: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q2_0)/8.031746),2)
  density,/ylog,f*1e5,q3_0,xran=[0,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 3: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q3_0)/8.031746),2)
  density,/ylog,f*1e5,q4_0,xran=[0,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 4: sqrt(!12<!3P!12>!3/8.031746) = '+strtrim(sqrt(mean(q4_0)/8.031746),2)

;The figure on the left shows the power spectrum of the first sector measured single ended, 
;There is a strong 1/f component that extends up to about 3KHz. On the right the same figure for the reference output, 
;quite noisy as well.

;Calculate the Cross correlation spectrum; with 4 sectors, we have 6 spectra, one for each pair of sectors. 
;Again, they are averaged over 108 ramps 
if Cross_correlation EQ 1 THEN BEGIN
c12_0 = fltarr(2048L*row)
c13_0 = fltarr(2048L*row)
c14_0 = fltarr(2048L*row)
c23_0 = fltarr(2048L*row)
c24_0 = fltarr(2048L*row)
c34_0 = fltarr(2048L*row)

    c12_0 += total(data0[*,*,1]*conj(data0[*,*,0]),2)/s[3]
    c13_0 += total(data0[*,*,2]*conj(data0[*,*,0]),2)/s[3]
    c14_0 += total(data0[*,*,3]*conj(data0[*,*,0]),2)/s[3]
    c23_0 += total(data0[*,*,2]*conj(data0[*,*,1]),2)/s[3]
    c24_0 += total(data0[*,*,3]*conj(data0[*,*,1]),2)/s[3]
    c34_0 += total(data0[*,*,3]*conj(data0[*,*,2]),2)/s[3]

  
; density,/ylog,f*1e5,c12_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 1x2'
; density,/ylog,f*1e5,c13_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 1x3'
; density,/ylog,f*1e5,c14_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 1x4'
; density,/ylog,f*1e5,c23_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 2x3'
; density,/ylog,f*1e5,c24_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 2x4'
; density,/ylog,f*1e5,c34_0,xran=[8e-7,0.5]*1e5,psym=3,yra=[1e1,1e7],charsiz=1,xsty=1,/dlog,/xlog,xminor=9,ct=-1,dsize=[2048,2048]/4,title=strmid(file[0],0,18)+'*:   Output 3x4'
;
;
;
;
;
;     5.9  Calculate the spectrum of the correlation matrix, i.e. the complex coherence matrix 

; Calculate correlation matrix
 ; if (option) then begin 
    h = where(f5 ne 0 and f5 lt 2e3,nh)
    mqq0 =[[1.0                                    ,total(c12_0[h]/sqrt(q1_0[h]*q2_0[h]))/nh,total(c13_0[h]/sqrt(q1_0[h]*q3_0[h]))/nh,total(c14_0[h]/sqrt(q1_0[h]*q4_0[h]))/nh],$
         [total(c12_0[h]/sqrt(q2_0[h]*q1_0[h]))/nh,1.0                                      ,total(c23_0[h]/sqrt(q2_0[h]*q3_0[h]))/nh,total(c24_0[h]/sqrt(q2_0[h]*q4_0[h]))/nh],$
         [total(c13_0[h]/sqrt(q3_0[h]*q1_0[h]))/nh,total(c23_0[h]/sqrt(q3_0[h]*q2_0[h]))/nh,1.0                                      ,total(c34_0[h]/sqrt(q3_0[h]*q4_0[h]))/nh],$
         [total(c14_0[h]/sqrt(q4_0[h]*q1_0[h]))/nh,total(c24_0[h]/sqrt(q4_0[h]*q2_0[h]))/nh,total(c34_0[h]/sqrt(q4_0[h]*q3_0[h]))/nh, 1.0                                      ]]
endif
         
;And the result is 
print,float(mqq0);
;      1.00000     0.545643     0.558868     0.564672
;     0.545643      1.00000     0.591492     0.594834
;     0.558868     0.591492      1.00000     0.611869
;     0.564672     0.594834     0.611869      1.00000
  print,'mark8: ',systime(1)-st

;
;
;5.23 Row Reference pixel subtraction.
  data4 = float(fft(data0,dim=1,1,/overwrite))/sqrt((s[1]/4+pad)*s[2])
  print,'mark15: ',systime(1)-st
  stop
  
  data4 = reform(data4,(s[1]/4+pad),s[2],s[3],4,/over)
  data4 = data4[0:511,*,*,*]
  data4[*,*,*,1] = data4[rev,*,*,1]
  data4[*,*,*,3] = data4[rev,*,*,3]
  data4 = transpose(data4,[0,3,1,2])
  data4 = reform(data4,s[1],s[2],s[3],/over)

  ; Row Reference pixel subtraction.
  ; Extract corrected reference pixels and correlated with data at low freq.
  rpix = [0,1,2,3] ; [0,1,2,3,2044,2045,2046,2047]
  nref = n_elements(rpix)
  even = findgen(nref/2)*2
  h1 = [lindgen(s[2]/2),lindgen(s[2]/2)+(s[2]*524)-s[2]/2]  ;0...1023: 2 bottom rows?
  h2 = [lindgen(s[2]/2)+(s[2]*524)/2,lindgen(s[2]/2)+(s[2]*524)/2-s[2]/2]   ; 2 top rows...?
  ref1a = total(data4[rpix,*,*],1)/nref   ;average of the 4 leftmost columns [2048, 108]
  ref2a = (total(data4[rpix[even],*,*],1)-total(data0[rpix[even+1],*,*],1))/nref  ;average of difference pairs column [0,2] - column [1,3] : odd/even correction
  ref1a = fft(ref1a,dim=1,/overwrite)*sqrt(s[2])   ;fft of the refpix signal
  ref2a = fft(ref2a,dim=1,/overwrite)*sqrt(s[2])   ;fft of the odd/even pair difference
  rat0ref1a = total(data4[h1,*,0]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  rat1ref1a = total(data4[h1,*,1]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  ;rat2ref1a = total(data5[h1,*,2]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  ;rat3ref1a = total(data5[h1,*,3]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  rat0ref2a = total(data4[h2,*,0]*conj(ref2a),2)/total(abs(ref2a)^2,2)
  rat1ref2a = total(data4[h2,*,1]*conj(ref2a),2)/total(abs(ref2a)^2,2)
  ;rat2ref2a = total(data5[h2,*,2]*conj(ref2a),2)/total(abs(ref2a)^2,2)
  ;rat3ref2a = total(data5[h2,*,3]*conj(ref2a),2)/total(abs(ref2a)^2,2)
;  
;  
;5.24  Subtract reference pixel corrections for low frequencies;
;
  ; Subtract reference pixel corrections for low frequencies.
for i=0,s[3]-1 do data4[h1,i,0] -= float(rat0ref1a)*ref1a[*,i]
  for i=0,s[3]-1 do data4[h1,i,1] -= float(rat1ref1a)*ref1a[*,i]
  ;for i=0,s[3]-1 do data5[h1,i,2] -= float(rat2ref1a)*ref1a[*,i]
  ;for i=0,s[3]-1 do data5[h1,i,3] -= float(rat3ref1a)*ref1a[*,i]
;  for i=0,s[3]-1 do data4[h2,i,0] -= float(rat0ref2a)*ref2a[*,i]
;  for i=0,s[3]-1 do data4[h2,i,1] -= float(rat1ref2a)*ref2a[*,i]
  ;for i=0,s[3]-1 do data5[h2,i,2] -= float(rat2ref2a)*ref2a[*,i]
  ;for i=0,s[3]-1 do data5[h2,i,3] -= float(rat3ref2a)*ref2a[*,i]
  
  rpix = [2047,2046,2045,2044] ; [0,1,2,3,2044,2045,2046,2047]
  nref = n_elements(rpix)
  even = findgen(nref/2)*2
  h1 = [lindgen(s[2]/2),lindgen(s[2]/2)+(s[2]*(512+pad))-s[2]/2]
  h2 = [lindgen(s[2]/2)+(s[2]*524)/2,lindgen(s[2]/2)+(s[2]*524)/2-s[2]/2]
  ref1a = total(data4[rpix,*,*],1)/nref
  ref2a = (total(data4[rpix[even],*,*],1)-total(data4[rpix[even+1],*,*],1))/nref
  ref1a = fft(ref1a,dim=1,/overwrite)*sqrt(s[2])
  ref2a = fft(ref2a,dim=1,/overwrite)*sqrt(s[2])
  ;rat0ref1a = total(data5[h1,*,0]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  ;rat1ref1a = total(data5[h1,*,1]*conj(ref1a),2)/total(abs(ref1a)^2,2)
;  rat2ref1a = total(data5[h1,*,2]*conj(ref1a),2)/total(abs(ref1a)^2,2)
;  rat3ref1a = total(data5[h1,*,3]*conj(ref1a),2)/total(abs(ref1a)^2,2)
  ;rat0ref2a = total(data5[h2,*,0]*conj(ref2a),2)/total(abs(ref2a)^2,2)
  ;rat1ref2a = total(data5[h2,*,1]*conj(ref2a),2)/total(abs(ref2a)^2,2)
;  rat2ref2a = total(data5[h2,*,2]*conj(ref2a),2)/total(abs(ref2a)^2,2)
;  rat3ref2a = total(data5[h2,*,3]*conj(ref2a),2)/total(abs(ref2a)^2,2)

end