; $Id: mips246_darkcal.pro,v 0.1 2003/01/31 15:00:00 dkelly Exp $
;
;+
; NAME:
;	MIPS246_DARKCAL
;
; PURPOSE:
;       This routine was written to process stim repeatability data
;       from IOC task mips-246.  It works from the output of mips_caler.
;       It ignores stim flash DCEs and a user-specified number of DCEs
;       after the stim and then calculates the mean brightness of 70um
;       or 160um array data on a module-by-module basis.  The results
;       are written out to file.  The full-array averages are plotted
;       to screen, providing a visual guide for deciding how many DCEs
;       to skip after a stim flash.
;
; CATEGORY:
;	MIPS utilities.
;
; CALLING SEQUENCE:
;	MIPS246_DARKCAL, filebase, ymin=ymin, ymax=ymax, skip=skip
;
; INPUTS:
;       Filebase: Base filename of a mips_caler .cal.fits file.
;	
;
; KEYWORD PARAMETERS:
;       YMIN: minimum plot limit for the y-axis
;       YMAX: maximum plot limit for the y-axis
;       SKIP: number of DCEs to skip after a stim in calculating means
;
; OUTPUTS:
;	This routine outputs numbers for the array mean, for the mean
;       over each half of the array, and also means for each of the
;       individual modules.
;
; OPTIONAL OUTPUTS:
;
; EXAMPLE:      
;               mips246_darkcal, 'mips_2001-267T15.23.05.888_0_A70',
;               ymin=5000, ymax=15000, skip=2
;
; MODIFICATION HISTORY:
; 	Written by:	Doug Kelly (June 29, 2003)
;       Based on:  mips246_stimcal.pro   and   mean_slope_70Modules.pro
;			
;-
; ----------------------------------------------------------------------

@kplot

PRO mips246_darkcal,filebase,ymin=ymin,ymax=ymax,skip=skip

; open a file into which dark level results will be written
openw,unit1,filebase+'.darklevels',/get_lun

if (not keyword_set(ymin)) then ymin = 0.0
if (not keyword_set(ymax)) then ymax = 1.0
if (not keyword_set(skip)) then skip = 0
printf,unit1,'Num DCEs skipped after each stim flash = ' + strtrim(string(skip),2)

calcycle = 6   ; DCEs/cycle -- be sure to set to  6 for mips-246
numstims = 47   ; be sure to set to 47 for mips-246
firstStimOffset = 1    ; set to  1 for mips-246

num_darks = (numstims-1)*(calcycle-1-skip)+firstStimOffset
ave_counts = fltarr(num_darks)
exp_vals = findgen(num_darks) + 1
kplot,exp_vals,ave_counts,psym=0,/nodata, $
  yrange=[ymin,ymax],xrange=[0,num_darks+3], $
  kplot_type='ii',xtitle='Dark DCE number',ytitle='counts [DN/sec]', $
  color=1,background=255

; get the reduced data

fits_open,filebase+'.cal.fits',fcb
; note: fcb.nextend=122 for this data set
printf,unit1,'Filename = ',filebase+'.cal.fits'
printf,unit1,'    module    counts/sec   read noise   max slope'
if (firstStimOffset > 0) then begin
    for j = 1,(firstStimOffset) do begin
      fits_read,fcb,cube,header,exten_no=j  
      if (j EQ 1) then begin
          size_cube = size(cube)
          meta_cube = fltarr(size_cube[1],size_cube[2],size_cube[3],num_darks)
          x_npts = size_cube[1]
          y_npts = size_cube[2]
      endif
      meta_cube[*,*,*,j-1] = cube
    endfor
endif
for j = 0,(numstims-2) do begin
  for k = (1+skip),(calcycle-1) do begin
    fits_read,fcb,cube,header,exten_no=calcycle*j+1+firstStimOffset+k  
    if (j EQ 0) AND (k EQ (1+skip)) AND (firstStimOffset EQ 0) then begin
        size_cube = size(cube)
        meta_cube = fltarr(size_cube[1],size_cube[2],size_cube[3],num_darks)
        x_npts = size_cube[1]
        y_npts = size_cube[2]
    endif
    meta_cube[*,*,*,((calcycle-1-skip)*j+firstStimOffset+(k-skip-1))] = cube
  endfor
endfor
fits_close,fcb
; now that all of the data have been read in for this particular file,
; calculate and plot the array mean for each DCE in the set.

; Now go ahead and calculate means and standard deviations
if (y_npts LT 32) then begin
  for j = 0,(num_darks-1) do begin
      ave_counts[j] = total(meta_cube[*,0,0,j],/NAN)
      index = where(finite(meta_cube[*,0,0,j]),n_pix0)
      ave_counts[j] = ave_counts[j] + total(meta_cube[*,2,0,j],/NAN)
      index = where(finite(meta_cube[*,2,0,j]),n_pix2)
      n_pixels = n_pix0 + n_pix2
      ave_counts[j] = ave_counts[j] / n_pixels
;      print, 'pixel ', j, ', npixels =', n_pixels
  endfor
endif else begin
  for j = 0,(num_darks-1) do begin
      ave_counts[j] = total(meta_cube[*,*,0,j],/NAN)
      index = where(finite(meta_cube[*,*,0,j]),n_pixels)
      ave_counts[j] = ave_counts[j] / n_pixels
  endfor
endelse
printf,unit1,'Mean flux =',mean(ave_counts),'+/-',stddev(ave_counts)
kplot,exp_vals,ave_counts,psym=100, $
  yrange=[ymin,ymax],xrange=[0,num_darks+5], kplot_type='ii', $
  xtitle='Dark DCE number',ytitle='counts [DN/sec]',color=1,background=255
;ans = ''
;read,'Continue: ',ans
; wait,2

free_lun,unit1

END 

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