; $Id: mips246_stimcal.pro,v 0.1 2003/01/31 15:00:00 dkelly Exp $
;
;+
; NAME:
;	MIPS246_STIMCAL
;
; PURPOSE:
;       This routine was written to process stim repeatability data
;       from IOC task mips-246.  It works from the output of mips_caler.
;       It calculates mean values for each of the stim DCEs for each of
;       the five biases and plots them together on one plot, showing
;       the stability of the stim values and the relationship between 
;       stim amplitude and bias voltage.  It also produces a table of
;       stim brightnesses and calculates the mean stim brightness for
;       each bias voltage and the stim repeatability.  Detailed comments
;       can be found below in the IDL code.
;
; CATEGORY:
;	MIPS utilities.
;
; CALLING SEQUENCE:
;	MIPS246_STIMCAL, ymin=ymin, ymax=ymax, skip=skip
;
; INPUTS:
;	
;
; KEYWORD PARAMETERS:
;       YMIN: minimum plot limit for the y-axis
;       YMAX: maximum plot limit for the y-axis
;       SKIP: number of DCEs to skip in calculating mean stim brightness
;
; OUTPUTS:
;	This routine outputs a single number for the array mean over the
;       set of DCEs.
;
; OPTIONAL OUTPUTS:
;
; EXAMPLE:      
;               mips246_stimcal, ymin=5000, ymax=15000, skip=10
;
; MODIFICATION HISTORY:
; 	Written by:	Doug Kelly (June 27, 2003)
;			
;-
; ----------------------------------------------------------------------

@kplot

PRO mips246_stimcal,ymin=ymin,ymax=ymax,skip=skip

; open a file into which stim brightness results will be written
openw,unit1,'mips246.stimbrightnesses',/get_lun

; produce a filelist of reduced data, *.cal.fits
filelist = findfile('*.cal.fits',count=numfiles)
print,filelist

if (not keyword_set(ymin)) then ymin = 0
if (not keyword_set(ymax)) then ymax = 20000
if (not keyword_set(skip)) then skip = 0
print,'ymin = '+ strtrim(string(ymin),2)
print,'ymax = '+ strtrim(string(ymax),2)
print,'skip = '+ strtrim(string(skip),2)

calcycle = 6    ; 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_odds = (numstims-skip-1)/2
ave_counts = fltarr(numstims-skip)
norm_counts = fltarr(num_odds)
exp_vals = findgen(numstims-skip) + 1
kplot,exp_vals,ave_counts,psym=0, $
  yrange=[ymin,ymax],xrange=[0,numstims-skip+3], $
  kplot_type='ii',xtitle='Stimflash number',ytitle='counts [DN/sec]', $
  color=1,background=255,/nodata

; get the reduced data

for i =0,(numfiles-1) do begin
  fits_open,filelist[i],fcb
; note: fcb.nextend=122 for this data set
  printf,unit1,'Filename = ',filelist[i]
  printf,unit1,'Stim Flash    ','Average Counts'
  for j = 0,(numstims-skip-1) do begin
    fits_read,fcb,cube,header,exten_no=calcycle*(j+skip)+1+firstStimOffset
    if (j EQ 0) then begin
        size_cube = size(cube)
        meta_cube = fltarr(size_cube[1],size_cube[2],size_cube[3],numstims)
        x_npts = size_cube[1]
        y_npts = size_cube[2]
    endif
    meta_cube[*,*,*,j] = cube
  endfor
  fits_close,fcb
  if (y_npts LT 32) then begin
    for j = 0,(numstims-skip-1) do begin
        ave_counts[j] = total(meta_cube[*,0,3,j],/NAN)
        index = where(finite(meta_cube[*,0,3,j]),n_pix0)
        ave_counts[j] = ave_counts[j] + total(meta_cube[*,2,3,j],/NAN)
        index = where(finite(meta_cube[*,2,3,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,(numstims-skip-1) do begin
;        ave_counts[j] = total(meta_cube[*,*,3,j],/NAN)
;        index = where(finite(meta_cube[*,*,3,j]),n_pixels)
;        ave_counts[j] = ave_counts[j] / n_pixels
        ave_counts[j] = total(meta_cube[20:31,*,3,j],/NAN)
        index = where(finite(meta_cube[20:31,*,3,j]),n_pix_123)
        ave_counts[j] = ave_counts[j] + total(meta_cube[16:19,8:31,3,j],/NAN)
        index = where(finite(meta_cube[16:19,8:31,3,j]),n_pix_4)
        n_pixels = n_pix_123 + n_pix_4
        ave_counts[j] = ave_counts[j] / n_pixels
        printf,unit1,j,ave_counts[j], format='(4x,I2,8x,F9.3)'
    endfor
    printf,unit1,'Mean flux =',mean(ave_counts),'+/-',stddev(ave_counts)
    for j = 0,(num_odds-1) do begin
       x = 2*j+1 ; ave_count index
       norm_counts[j] = 2*ave_counts[x]/(ave_counts[x-1]+ave_counts[x+1])
    endfor
    printf,unit1,'Std Dev of stim-calibrated flashes =',stddev(norm_counts)
    kplot,exp_vals,ave_counts,psym=(2*i), $
       yrange=[ymin,ymax],xrange=[0,numstims-skip+3], $
;      kplot_type='ii',xtitle='Stimflash number',ytitle='counts [DN/sec]', $
       color=1,background=255,/noerase
  endelse
endfor

free_lun,unit1


END 

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