; $Id: vrst_160module_meanslope.pro,v 0.1 2003/01/31 15:00:00 dkelly Exp $
;
;+
; NAME:
;	VRST_160MODULE_MEANSLOPE
;
; PURPOSE:
;       This routine works from the output of mips_sloper and calculates 
;       the mean slope for each module in each extension.  The mean values
;       are written out to file.  
;
; CATEGORY:
;	MIPS utilities.
;
; CALLING SEQUENCE:
;	vrst_160module_meanslope, Filebase
;
; INPUTS:
;	Filebase: Base filename of a mips_sloper .red.fits file.  
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;	This routine outputs mean slopes for each module in each extension
;       of the .red.fits file.  The output filename is filebase.meanslopes 
;
; OPTIONAL OUTPUTS:
;
; MODIFICATION HISTORY:
;       Based on mean_slope.pro, written in Jan 2003.
; 	Written by:	Doug Kelly (October 31, 2003)
;			
;-
; ----------------------------------------------------------------------

PRO vrst_160module_meanslope,filebase

; 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

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

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

ave_counts = fltarr(fcb.nextend,4,3)

; 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,(fcb.nextend-1) do begin
 for j = 0, 3 do begin
   mincol = 5*j
   maxcol = 5*j + 4
  for k = 0, 1 do begin
     ave_counts[i,j,k] = total(meta_cube[mincol:maxcol,2*k,0,i],/NAN)
     index = where(finite(meta_cube[mincol:maxcol,2*k,0,i]),n_pixels)
     ave_counts[i,j,k] = ave_counts[i,j,k] / n_pixels
  endfor
     ave_counts[i,j,2] = (ave_counts[i,j,0] + ave_counts[i,j,1])/2.0
     if (j EQ 2) then ave_counts[i,j,2] = ave_counts[i,j,1]
 endfor
endfor

; perform linear regression fits to Vrst data on each readout and module
opt_vrst = fltarr(4,3)
module_vrst = fltarr(4)
igood = indgen(11)*2 + 3
even_counts = ave_counts[igood,*,*]
bias_voltage = [5,4,3,2,1,0,-1,-2,-3,-4,-5]

for j = 0, 3 do begin
  for k = 0, 2 do begin
    if (j EQ 2) AND (k EQ 0) then begin
     opt_vrst[j,k] = 0.0
    endif else begin
     ifit = where(finite(even_counts[*,j,k]),n_good)
     if (n_good GT 0) then begin
       X = even_counts[ifit,j,k]
       Y = bias_voltage[ifit]
       measure_errors = abs(1.0/even_counts[ifit,j,k])
       result = POLY_FIT(even_counts[ifit,j,k],bias_voltage[ifit],1, $
         MEASURE_ERRORS=measure_errors[ifit])
       opt_vrst[j,k] = result[0]
     endif else begin
       opt_vrst[j,k] = -300
       print,j,k,opt_vrst[j,k]
     endelse
    endelse
  endfor
  module_vrst[3-j] = total(opt_vrst[j,0:1])/2.0
  if (j EQ 2) then module_vrst[3-j] = opt_vrst[j,1]
endfor

; setup output file for readout report

openw,unit1,filebase+'.meanslopes',/get_lun
printf,unit1,''
print,''
printf,unit1,filebase
print,filebase

; do statistics by module

printf,unit1,''
;print,''
printf,unit1,'          ***array statistics (counts/sec)***'
;print,'          ***array statistics (counts/sec)***'
printf,unit1,'      DCE   module Readout1  Readout2'
;print,'      DCE   module Readout1  Readout2'
for j = 0, 3 do begin
  for i = 0,(fcb.nextend -1) do begin
    printf,unit1,i,4-j,ave_counts[i,j,0],ave_counts[i,j,1], $
      format="(2x,I6,2x,I6,2x,F8.1,2x,F8.1)"
  endfor
  print,'     Module  Opt.Vrst'
  print,j+1,module_vrst[j]
endfor

printf,unit1,''
printf,unit1,'     Module  Readout  Opt.Vrst'
for j = 0, 3 do begin
  for k = 0, 2 do begin
     printf,unit1,4-j,k,opt_vrst[j,k]
  endfor
endfor

printf,unit1,''
printf,unit1,'Vrst averaged over readouts in module:'
printf,unit1,'   1       2       3       4
printf,unit1,module_vrst,$
 format="(F5.2,3x,F5.2,3x,F5.2,3x,F5.2)"

free_lun,unit1

; plot mean slopes 
!p.background=255
!p.color=0
loadct,39
; configure to put four plots on one page, 2 columns x 2 rows
!p.multi = [0,2,2,0,1]

psland, fname=filebase+'.mean.ps', wide=10.5, high=7.5

for j = 0, 3 do begin
  MaxSlope = Max(ave_counts[*,j,*], min=MinSlope)
  if (j EQ 2) then MaxSlope = Max(ave_counts[*,j,1],min=MinSlope)
  plot, findgen(fcb.nextend), ave_counts[*,j,0], yrange=[(MinSlope*1.3),(MaxSlope*1.3)], $
     xtitle='DCE', ytitle='Mean Slope (DN/s)', $
     title='160um Vrst Optimization', $
     subtitle = (j+1), charsize=1.3, color=0,/nodata
  plots, [0,fcb.nextend], [0,0], linestyle=2, color=0
  plots, [13,13], [(MinSlope*1.3),(MaxSlope*1.3)], linestyle=2, color=0
  for k = 0, 2 do begin
    oplot, findgen(fcb.nextend), ave_counts[*,j,k], color=(48*k+10)
;    xyouts, 15, (MinSlope*(0.5 + 0.1*k)),(k+1), color=(48*k+10)
  endfor
  xyouts, 15, (MinSlope*0.5),"R.O.1", color=(10)
  xyouts, 15, (MinSlope*0.65),"R.O.2", color=(58)
  xyouts, 15, (MinSlope*0.8),"Module", color=(106)
  if (j EQ 0) then begin
    xyouts, 3, MaxSlope*1.15, filebase, color=0
  endif
endfor

device,/close
set_plot, 'x'

END 

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