; $Id: mean_slope_70modules.pro,v 0.1 2003/01/31 15:00:00 dkelly Exp $
;
;+
; NAME:
;	MEAN_SLOPE_70MODULES
;
; PURPOSE:
;       This routine calculates the mean slope from a set of MIPS 70um
;       images.  The mean can be calculated over the full set of DCEs or
;       over a user-specified range of DCEs.  The mips_sloper
;	portion of the DAT must be run to provide input to this
;	procedure.  Any DCEs in which a pixel was affected by a cosmic
;	ray are not used in computing the mean for that pixel.
;
; CATEGORY:
;	MIPS utilities.
;
; CALLING SEQUENCE:
;	MEAN_SLOPE_70Modules, Filebase
;
; INPUTS:
;	Filebase: Base filename of a mips_sloper .red.fits file.  
;
; KEYWORD PARAMETERS:
;       NO_CR:  The default is to average only the DCEs in which the pixel
;               in question was not struck by a cosmic ray.
;               If you want to include all DCEs without regard to cosmic ray
;               hits, set the keyword NO_CR.
;   N_DCE_REJ:  If you want to exclude the first x DCEs from the average,
;               set N_DCE_REJ = x.
;     MAX_DCE:  If you want to exclude all DCEs above number y from the
;               average, set MAX_DCE = y.
;
; OUTPUTS:
;	This routine outputs numbers for the array mean over the set of 
;       DCEs, for the mean over each side of the array, and means for
;       each of the 8 modules.
;
; OPTIONAL OUTPUTS:
;
; EXAMPLE:   If you are averaging a set of 10 DCEs, you want to exclude
;            the first 4 and the last DCE from the average, and you want 
;            to include pixels that were hit by cosmic rays:
;
;            mean_slope_70modules, n_dce_rej = 4, max_dce=9, /no_cr, 'filebase'
;
; MODIFICATION HISTORY:
;       Based on mean_slope.pro, written in Jan 2003.
; 	Written by:	Doug Kelly (June 25, 2003)
;			
;-
; ----------------------------------------------------------------------

@kplot
@comp_std

PRO mean_slope_70modules,filebase,no_cr=no_cr,n_dce_rej=n_dce_rej,max_dce=max_dce

; get the reduced data

fits_open,filebase+'.red.fits',fcb

fits_read,fcb,trash,header,exten_no=0  ; get the main header for the file

if (keyword_set(max_dce)) then begin
  end_reject = fcb.nextend - max_dce
  fcb.nextend = min([max_dce,fcb.nextend])
endif else begin
  end_reject = 0
endelse

if (not keyword_set(n_dce_rej)) then n_dce_rej = 0
n_dce_begin = 1 + n_dce_rej
n_exp = fcb.nextend - n_dce_rej

for i = n_dce_begin, fcb.nextend do begin
    fits_read,fcb,cube,exten_no=i
    if (i EQ n_dce_begin) then begin
        size_cube = size(cube)
        meta_cube = fltarr(size_cube[1],size_cube[2],size_cube[3],n_exp)
    endif
    meta_cube[*,*,*,i-1-n_dce_rej] = cube
endfor
fits_close,fcb

x_npts = size_cube[1]
y_npts = size_cube[2]

exp_vals = findgen(n_exp) + 1
ave_counts = fltarr(n_exp)

; This block calculates and then plots the array mean for each DCE 
; in the set.  It is useful for seeing if there are drifts due to
; responsivity changes or latents.

for i = 0,(n_exp-1) do begin
    ave_counts[i] = total(meta_cube[*,*,0,i],/NAN)
    index = where(finite(meta_cube[*,*,0,i]),n_pixels)
    ave_counts[i] = ave_counts[i] / n_pixels
    if (keyword_set(no_cr)) then begin
        n_pixels = x_npts*y_npts
    endif
endfor

kplot,exp_vals,ave_counts,psym=100, $;,yrange=[1.9e4,2e4],xrange=[1,60], $
  kplot_type='ii',xtitle='exposure',ytitle='counts [DN/sec]',color=1,background=255
wait,2
;kplot,[0],[0],xrange=[1,n_exp],yrange=[0,20000],/nodata

;scrn2gif,filebase+'_tot.gif'

count_image = fltarr(x_npts,y_npts)
rdnoise_image = fltarr(x_npts,y_npts)
noise_image = fltarr(x_npts,y_npts)

ave_npts_used = 0.0
n_shot_noise = 0
ave_shot_noise = 0.0
;cos_val = 64
for i = 0,(x_npts-1) do begin
  for j = 0,(y_npts-1) do begin 

      indxs = where(((meta_cube[i,j,2,*] / 64) mod 2) LT 1,n_indxs)
      if (n_indxs lt 1) then $
          print,'Pixel ',i,j,' :  all DCEs rejected'
      if (keyword_set(no_cr)) then begin
          indxs = indgen(n_exp)
          n_indxs = n_exp
      endif
      ave_npts_used = ave_npts_used + n_indxs
;       print,i,j,n_indxs

; compute average counts [total electrons] for each pixel

      if (n_indxs GT 0) then begin
          counts = total(meta_cube[i,j,0,indxs],/NAN)/float(n_indxs)
      endif else begin
          counts = 0.0
      endelse
      count_image[i,j] = counts
  endfor
endfor 
ave_npts_used = ave_npts_used/n_pixels

; output counts array

fxhmake,header,count_image
fxwrite,filebase+'.cnt.fits',header,count_image

; setup output file for readout report

openw,unit1,filebase+'.cnt.report',/get_lun
printf,unit1,filebase
printf,unit1,'Num DCEs rejected at start = ' + strtrim(string(n_dce_rej),2)
printf,unit1,'Num DCEs rejected at end   = ' + strtrim(string(end_reject),2)
printf,unit1,'counts/sec = electrons/sec'
printf,unit1,"average # exposures used = " +  $
  strtrim(string(ave_npts_used,format="(F5.2)"),2) + " per pixel"

print,'Num DCEs rejected at start = ' + strtrim(string(n_dce_rej),2)
print,'Num DCEs rejected at end   = ' + strtrim(string(end_reject),2)
print,'counts/sec = electrons/sec'
print,"average # exposures used = " + $
  strtrim(string(ave_npts_used,format="(F5.2)"),2) + " per pixel"

; do statistics by module, array-side, and full array

ar_counts_std = fltarr(11)
ar_counts = fltarr(11)
max_slope = fltarr(11)

for i = 0, 7 do begin
   mincol = 4*i
   maxcol = 4*i + 3
   ar_counts_std[i] = comp_std(count_image[mincol:maxcol,*],ave=ave_val)
   ar_counts[i] = ave_val
   max_slope[i] = max(count_image[mincol:maxcol,*], /NAN)
endfor
for i = 0, 1 do begin
   mincol = 16*i
   maxcol = 16*i + 15
   ar_counts_std[8+i] = comp_std(count_image[mincol:maxcol,*],ave=ave_val)
   ar_counts[8+i] = ave_val
   max_slope[8+i] = max(count_image[mincol:maxcol,*], /NAN)
endfor
ar_counts_std[10] = comp_std(count_image[*,*],ave=ave_val)
ar_counts[10] = ave_val
max_slope[10] = max(count_image, /NAN)

printf,unit1,''
printf,unit1,'***array statistics***'
printf,unit1,'    module    counts/sec   read noise   max slope'
for i = 0, 7 do begin
  printf,unit1,i+1,ar_counts[i],ar_counts_std[i], max_slope[i], $
    format="(4x,I6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
endfor
printf,unit1,"side A",ar_counts[8],ar_counts_std[8], max_slope[8], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
printf,unit1,"side B",ar_counts[9],ar_counts_std[9], max_slope[9], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
printf,unit1,"Array",ar_counts[10],ar_counts_std[10], max_slope[10], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"

print,'***array statistics***'
print,'     module     counts/sec   read noise   max slope'
for i = 0, 7 do begin
  print,i+1,ar_counts[i],ar_counts_std[i], max_slope[i], $
    format="(4x,I6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
endfor
print,"side A",ar_counts[8],ar_counts_std[8], max_slope[8], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
print,"side B",ar_counts[9],ar_counts_std[9], max_slope[9], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"
print,"Array",ar_counts[10],ar_counts_std[10], max_slope[10], $
  format="(4x,A6,4x,F8.1,' +/- ',F8.1,5x,F8.1)"

free_lun,unit1

; compute histogram

kxrange = krange(count_image[*,*])
min_histo = kxrange[0]
max_histo = kxrange[1]
npts_histo = 100
binsize = (max_histo - min_histo)/npts_histo
npts_histo = (max_histo - min_histo)/binsize
val_bins = min_histo + (findgen(npts_histo) + 0.5)*binsize
histo = histogram(count_image[*,*],binsize=binsize,min=min_histo,max=max_histo)

kxrange = krange(val_bins)
kxrange[0] = 0.0
kyrange = krange(histo)
kyrange[0] = 0.0
kplot,val_bins,histo,xrange=kxrange,yrange=kyrange, $
  xtitle='counts [e!U-!N/sec]',ytitle='# pixels', $
  position=[0.15,0.15,0.9,0.9],title=filebase,color=1,background=255

;scrn2gif,filebase+'_cnt.gif'

END 

; ----------------------------------------------------------------------