; A Dave Walton IDL program for Light Transfer Curves (LTCs). It uses either
;drift-scan fits files or pairs
;of 20x20 pixel images from CICADA .dat files and works out the signal and
;standard deviation. It does this for left/right and bottom/top ports
;and for low and high gains.

pro addsysvars
  defsysv, '!TEXTOUT', 1
  defsysv, '!TEXTUNIT', 0
end


device, retain=2
!p.multi=[0,1,2]

;First, read the data
flnam=pickfile(title='Select .dat or .fit file', /must_exist)

wNlr=2

;Labels
portstring=[[' BL', ' BR'], $
            [' TL', ' TR']]

gainstring=[' NormGain', ' HiGain']


;Extract the simple filename
if strupcase(!version.os_family) EQ 'WINDOWS' then cSeparator=92b $
  else cSeparator=47b ;Forward or back-slash
whereslash=where(byte(flnam) eq cseparator)
whereslash=whereslash(n_elements(whereslash)-1) ;Last slash
shortflnam=strmid(flnam, whereslash+1, 1000)

wDotpos=strpos(shortflnam, '.');
flext=strmid(shortflnam, wDotpos, strlen(shortflnam)-wDotpos)

if strlowcase(flext) EQ '.fit' then begin
  addsysvars
  fits_info, flnam, N_ext=wN_ext
  if wN_ext EQ 1 then wNbt=1
  temp=readfits(flnam, header, ext=0)
  short_ccd_type=strupcase(strmid(strtrim(sxpar(header, 'CCD_TYPE'), 2), 0, 5))
  case short_ccd_type of ;Pixel size, cm
    'CCD42': fPix_cm=13.5e-4
    'CCD47': fPix_cm=13.0e-4
    'CCD64': fPix_cm=16.0e-4
    'CCD74': fPix_cm=12.0e-4
    'CCD75': fPix_cm=12.0e-4
    'CCD80': fPix_cm=12.0e-4
    else: begin
            print, 'Unknown CCD type, Enter pixel size(microns): '
            read, fPix_cm
            fPix_cm=1.0e-4*fPix_cm
          end
  endcase
  fTintincr_logain=sxpar(header, 'EXPOSURE')
  s=size(temp)
  case s[2] of
    140: wNy=20
    230: wNy=10
    else: wNy=1
  endcase
  wNtint=s[2]/wNy
  wNx=s[1]
  case wN_ext of
    1: begin
         wHigainflag=0
         wNgain=1 & wNbt=1
       end
    3: begin
         case short_ccd_type of
           'CCD74': begin ;4 ports
                      wHigainflag=0
                      wNgain=1 & wNbt=2
                    end
           'CCD75': begin ;4 ports
                      wHigainflag=0
                      wNgain=1 & wNbt=2
                    end
           else:    begin ;2 ports
                      wHigainflag=1
                      wNgain=2 & wNbt=1
                    end
         endcase
       end
    7: begin
         wHigainflag=1
         wNgain=2 & wNbt=2
       end
    else: begin
            wHigainflag=0
            wNgain=1 & wNbt=1
          end
  endcase

  warrAll_data=intarr(wNx, wNy, 2*wNtint, wNlr, wNbt, wNgain)

  wI_ext=0
  for wIgain=0, wNgain-1 do begin
    for wIbt=0, wNbt-1 do begin
      for wIlr=0, wNlr-1 do begin
        temp=readfits(flnam, header, ext=wI_ext)
        warrAll_data[*, *, 2*indgen(wNtint), wIlr, wIbt, wIgain]= $
          reform(temp[*, *, 0], wNx, wNy, wNtint)
        warrAll_data[*, *, 2*indgen(wNtint)+1, wIlr, wIbt, wIgain]= $
          reform(temp[*, *, 1], wNx, wNy, wNtint)

        wI_ext=wI_ext+1 ;For next time through
      endfor
    endfor
    gainstring[wIgain]=' Gain='+ $
      strtrim(string(sxpar(header, 'AMPAGAIN'), format='(F5.2)'), 2)
  endfor
  fTintincr_higain=sxpar(header, 'EXPOSURE')

  ;Values to print to file
  case strupcase(strmid(shortflnam, 3, 1)) of
    'R': begin
           fVrd=float(sxpar(header, 'V_RD'))
           fVod=float(sxpar(header, 'V_OD'))
           parmstring= ' Vrd=' + strtrim(string(fVrd, format='(F4.1)'), 2) + $
           ', Vod=' + strtrim(string(fVod, format='(F4.1)'), 2)
         end

    'S': begin
           fVss=float(sxpar(header, 'V_SS'))
           parmstring= ' Vss=' + strtrim(string(fVss, format='(F4.1)'), 2)
         end

    else: parmstring=''
  endcase
 endif $

 else begin
  fPixcm=16.0e-4 ;CCD64
  wNtint=7
  wNx=20 & wNy=20
  fTintincr_logain=.996
  fTintincr_higain=.2097
  ;wNx x wNy window
  openr, 101, flnam
  stats=fstat(101)
  if stats.size EQ 44800L then begin
    wHigainflag=1
    wNgain=2
   endif else begin
    wHigainflag=0
    wNgain=1
  endelse
  warrAll_data=intarr(wNx, wNy, 2*wNtint, wNports, wNgain)

  readu, 101, warrAll_data
  close, 101

  if strupcase(!version.arch) EQ 'SPARC' then $
    byteorder, warrAll_data ;Swap bytes if a Sun

  if strupcase(strmid(shortflnam, 3, 1)) EQ 'R' then begin
    wVrd=16+2*fix(strmid(shortflnam, 4, 1))
    wVod=28+2*fix(strmid(shortflnam, 6, 1))
    parmstring= ' Vrd=' + strtrim(string(wVrd), 2) + ', Vod=' + $
                strtrim(string(wVod), 2)
  endif

  if strupcase(strmid(shortflnam, 3, 1)) EQ 'S' then begin
    wVss=5+fix(strmid(shortflnam, 4, 1))
    parmstring= ' Vss=' + strtrim(string(wVss), 2)
  endif
endelse


;Put two images side-by-side and tvscl, to look for structure
for wIbt=0, wNbt-1 do begin
for wIlr=0, wNlr-1 do begin
  for wIgain=0, wNgain-1 do begin
    im2=intarr(wNx+wNx+1, wNy)
    im2[0:(wNx-1), *]=warrAll_data[*, *, 0, wIlr, wIbt, wIgain]
    im2[(wNx+1):(wNx+wNx), *]=warrAll_data[*, *, 1, wIlr, wIbt, wIgain]
    im2[wNx,*]=min(warrAll_data[*, *, 0:1, wIlr, wIbt, wIgain])
    im3=rebin(im2, 10*(wNx+wNx+1), 10*wNy, /sample)
    window, wIgain, title=portstring[wIlr, wIbt] + gainstring[wIgain]
    tvscl, im3
    print, portstring[wIlr, wIbt] + gainstring[wIgain] + 'Column means'
    print, total(warrAll_data[*, *, 0, wIlr, wIbt, wIgain], 2)/wNy
    print
  endfor
  print, 'Enter number of columns to eliminate:'
  read, wNcols
  if wNcols GT 0 then begin
    wCols=intarr(wNcols)
    print, 'Enter columns to eliminate:'
    read, wCols
    print, 'Columns=:', wCols
  endif
endfor
endfor


;Sort out
;Declare arrays, window size etc. Window chosen to avoid (probably
;insignificant) bright and dark defects.
wXmin=0 & wXmax=wNx-1
wXrange=wXmax-wXmin+1 & fXrange=float(wXrange)
wYmin=0 & wYmax=wNy-1
wYrange=wYmax-wYmin+1 & fYrange=float(wYrange)
dSigmasquare=dblarr(wNtint)
dSigplusoff=dblarr(wNtint, wNlr, wNbt, wNgain)
fTint=fltarr(wNtint, wNgain)
fTint[*, 0]=fTintincr_logain*findgen(wNtint)
if wHigainflag NE 0 then fTint[*, 1]=fTintincr_higain*findgen(wNtint)
farrData0=fltarr(wXrange, wYrange)
farrData1=fltarr(wXrange, wYrange)

;Stats arrays etc.
fEperdn=fltarr(wNlr, wNbt, wNgain)
fYfit_sigplusoff=fltarr(wNtint, wNlr, wNbt, wNgain, /nozero)
fSigma_sigplusoff=fltarr(2, wNlr, wNbt, wNgain) ;Mean+SD, l+r, b+t, NormGain+HiGain
fSigma_sigmasquare=fltarr(2, wNlr, wNbt, wNgain) ;Mean+SD, l+r, b+t, NormGain+HiGain
fCurrent_dn=fltarr(wNlr, wNbt, wNgain) & fCurrent_e=fltarr(wNlr, wNbt, wNgain)
fCM=fltarr(2) ;C+M
fRead_noise_sq_dn2=fltarr(wNlr, wNbt, wNgain)
fRead_noise_dn2=fltarr(wNlr, wNbt, wNgain)
fADCoff=fltarr(wNlr, wNbt, wNgain)
fYfit_sigmasquare=fltarr(wNtint, wNlr, wNbt, wNgain)
dSigmasquare=dblarr(wNtint) & dSigplusoff=dblarr(wNtint, wNlr, wNbt, wNgain)
dSigplusresiduals=dblarr(wNtint, wNlr, wNbt, wNgain)
dSig=dblarr(wNtint, wNlr, wNbt, wNgain)
dSigmasquare=dblarr(wNtint, wNlr, wNbt, wNgain)
dSigmaresiduals=dblarr(wNtint, wNlr, wNbt, wNgain)
faZero=replicate(0.0, wNtint)

;Process the data
for wGain=0, wNgain-1 do begin
for wIbt=0, wNbt-1 do begin
for wIlr=0, wNlr-1 do begin
  for wI=0, wNtint-1 do begin
    farrData0= $
      float(warrAll_data[wXmin:wXmax, wYmin:wYmax, 2*wI, wIlr, wIbt, wGain])
    farrData1= $
      float(warrAll_data[wXmin:wXmax, wYmin:wYmax, 2*wI+1, wIlr, wIbt, wGain])
    ;Adjust each line to eliminate line-to-line variation. Note that this
    ;eliminates one degree of freedom from the data, hence the dSD_factor
    ;to correct the standard deviation
    dSD_factor=double(wXrange)/(wXrange-1)
    for wY=0, wYrange-1 do begin
      fLine_total0=total(farrData0[*, wY])
      fLine_total1=total(farrData1[*, wY])
      fHalfdiff=(fLine_total1-fLine_total0)/(2.0*fXrange)
      farrData0[*, wY]=farrData0[*, wY]+fHalfdiff
      farrData1[*, wY]=farrData1[*, wY]-fHalfdiff
    endfor

    dDiff=double(farrData1-farrData0)
    dSum=double(farrData0+farrData1)

    ;window, 1
    ;tvscl, rebin(dDiff, 200, 200, /sample)

    ;Calculate LTC stats
    dDiffsquare=dDiff*dDiff
    dSigmasquare[wI, wIlr, wIbt, wGain]= $
      total(dDiffsquare)/(2.0*fXrange*fYrange)*dSD_factor
    dSigplusoff[wI, wIlr, wIbt, wGain]= $
      total(dSum)/(2.0*fXrange*fYrange)
  endfor


  ;Plot ADC v. Tint. Then plot residuals
  ;Do line fit. Skip 1st point as it's not necessarily on the line.
  ;Also, skip any zero points
  whereynzero=where(dSigplusoff[*, wIlr, wIbt, wGain] GT 0.0d)
  wMinindex=fix(1L > whereynzero[0]) 
  fSigma_temp=0.0
  if min(dSigplusoff[wMinindex:(wNtint-1), wIlr, wIbt, wGain]) EQ $
     max(dSigplusoff[wMinindex:(wNtint-1), wIlr, wIbt, wGain]) then begin
    fCM=float([dSigplusoff[wMinindex, wIlr, wIbt, wGain], 0.0])
    fYfit_sigplusoff=replicate(float(dSigplusoff[wMinindex, wIlr, wIbt, wGain]), $
                               wNtint-1)
    fSigma_temp=[0.0, 0.0]
   endif else if (wNtint-wMinindex) EQ 2 then begin
     fCM[1]=float((dSigplusoff[wNtint-1, wIlr, wIbt, wGain]- $
                   dSigplusoff[wMinindex, wIlr, wIbt, wGain])/ $
                  (fTint[wNtint-1, wGain]-fTint[wMinindex, wGain]))
     fCM[0]=float(dSigplusoff[wNtint-1, wIlr, wIbt, wGain]- $
                  fCM[1]*fTint[wNtint-1, wGain])
     fYfit_sigplusoff= $
       float(dSigplusoff[wMinindex:(wNtint-1), wIlr, wIbt, wGain])
   endif else $
    fCM= svdfit(fTint[wMinindex:(wNtint-1), wGain], $
                dSigplusoff[wMinindex:(wNtint-1), wIlr, wIbt, wGain], $
                2, yfit=fYfit_sigplusoff, variance=fSigma_temp)

  fSigma_sigplusoff[*, wIlr, wIbt, wGain]=fSigma_temp

  if whereynzero[0] NE 0L then $ ;If some zeroes then
    fAdcoff[wIlr, wIbt, wGain]=fCM[0] $
   else fAdcoff[wIlr, wIbt, wGain]=dSigplusoff[0, wIlr, wIbt, wGain]
  fCurrent_dn[wIlr, wIbt, wGain]=fCM(1)

  ;Add eliminated points back in
  fYfit_sigplusoff=[fCM[1]*fTint[*, wGain]+fCM[0]]

  window, 2*wIlr+4*wIbt+8*wGain
  ;Plot it this way so ADCoff doesn't dominate
  plot, fTint[*, wGain], fYfit_sigplusoff, /nodata, $
          charsize=1.4, xtitle='T_Int', ytitle='ADC o/p(DN)', $
          title=flnam+portstring[wIlr, wIbt]+gainstring[wGain]

  oplot, fTint[*, wGain], dSigplusoff[*, wIlr, wIbt, wGain], psym=-1

  dSigplusresiduals[*, wIlr, wIbt, wGain]= $
                         dSigplusoff[*, wIlr, wIbt, wGain]-fYfit_sigplusoff
  ;Set pt.0 and other zeroes to zero temporarily
  dSigplusresidualstemp=dSigplusresiduals[*, wIlr, wIbt, wGain]
  dSigplusresidualstemp[0:(wMinindex-1)]=0.0d
  plot, fTint[*, wGain], dSigplusresidualstemp, /nodata, $
          charsize=1.4, xtitle='T_Int', ytitle='Residuals(DN)'
  oplot, fTint[*, wGain], dSigplusresiduals[*, wIlr, wIbt, wGain], psym=-1


  ;Plot Sigma^2 v. signal(DN). Then plot residuals
  window, 1 + 2*wIlr + 4*wIbt + 8*wGain
  dSig[*, wIlr, wIbt, wGain]=dSigplusoff[*, wIlr, wIbt, wGain] $
                               -fAdcoff[wIlr, wIbt, wGain]
  plot, dSig[*, wIlr, wIbt, wGain], dSigmasquare[*, wIlr, wIbt, wGain], $
          psym=-1, charsize=1.4, xtitle='Signal(DN)', ytitle='Sigma^2(DN^2)', $
          title=flnam+portstring[wIlr, wIbt]+gainstring[wGain]

  ;Line fit
  fyfit_temp=0.0
  if min(dSig[*, wIlr, wIbt, wGain]) EQ $
     max(dSig[*, wIlr, wIbt, wGain]) then begin
    fYfit_temp=faZero
    fSigma_temp=[0.0d, 0.0d]
   endif else begin
    fCM=svdfit(dSig[*, wIlr, wIbt, wGain], $
               dSigmasquare[*, wIlr, wIbt, wGain], 2, $
                    yfit=fYfit_temp, $
                    variance=fSigma_temp)
  endelse
  fSigma_temp=sqrt(fSigma_temp)

  fYfit_sigmasquare[*, wIlr, wIbt, wGain]=fYfit_temp
  fSigma_sigmasquare[*, wIlr, wIbt, wGain]=fSigma_temp

  ;fRead_noise_sq_dn2(wPort, wGain)=fCM(0) ;Fit version of read noise
  fRead_noise_sq_dn2[wIlr, wIbt, wGain]= $
                      dSigmasquare[0, wIlr, wIbt, wGain] ;Data pt. version
  if fCM[1] EQ 0.0 then fEperdn[wIlr, wIbt, wGain]=0.0 $
   else fEperdn[wIlr, wIbt, wGain]=1.0/fCM[1]

  oplot, dSig[*, wIlr, wIbt, wGain], fYfit_sigmasquare[*, wIlr, wIbt, wGain]

  dSigmaresiduals[*, wIlr, wIbt, wGain]= $
         dSigmasquare[*, wIlr, wIbt, wGain]-fYfit_sigmasquare[*, wIlr, wIbt, wGain]
  plot, dSig[*, wIlr, wIbt, wGain], dSigmaresiduals[*, wIlr, wIbt, wGain], $
          psym=-1, charsize=1.4, xtitle='Signal(DN)', ytitle='Residuals(DN^2)'
endfor
endfor
endfor


for wGain=0, wNgain-1 do begin
for wIbt=0, wNbt-1 do begin
print
print
print, flnam+gainstring[wGain]+portstring[0, wIbt]+portstring[1, wIbt]
print, '                          Left                     Right'
print, format='(A, 2(G10.3, A, G10.3))', $
       'For T_Int=0, ADC   =', fAdcoff[0, wIbt, wGain], ' +/-', $
                            fSigma_sigplusoff[0, 0, wIbt, wGain], $
                            fAdcoff[1, wIbt, wGain], ' +/-', $
                            fSigma_sigplusoff[0, 1, wIbt, wGain]
print, format='(A, 2(G10.3, A, G10.3))', 'Responsivity(e-/DN)=', $
                       fEperdn[0, wIbt, wGain], ' +/-', $
                       fSigma_sigplusoff[1, 0, wIbt, wGain], $
                       fEperdn[1, wIbt, wGain], ' +/-', $
                       fSigma_sigplusoff[1, 1, wIbt, wGain]
print

fCurrent_e=fCurrent_dn*fEperdn
print, format='(A, 2G24.3)', 'Signal (DN/s)      =', fCurrent_dn[*, wIbt, wGain]
print, format='(A, 2G24.3)', 'Signal (e-/s)      =', fCurrent_e[*, wIbt, wGain]
print, format='(A, 2G24.3)', 'Signal (nA cm-2)   =', fCurrent_e[*, wIbt, wGain] $
                                   *1.602e-19*1e9/(fPix_cm*fPix_cm)

print
print, format='(A, 2(G10.3, A, G10.3))', 'Rd Noise Sq(DN^2)  =', $
           fRead_noise_sq_dn2[0, wIbt, wGain], ' +/-', $
           fSigma_sigmasquare[0, 0, wIbt, wGain], $
           fRead_noise_sq_dn2[1, wIbt, wGain], ' +/-', $
           fSigma_sigmasquare[0, 1, wIbt, wGain]


fRead_noise_sq_e2=fRead_noise_sq_dn2*fEperdn*fEperdn
print, format='(A, 2(G10.3, A, G10.3))', 'Rd Noise Sq(e-^2)  =', $
 fRead_noise_sq_e2[0, wIbt, wGain], ' +/-', $
 fSigma_sigmasquare[0, 0, wIbt, wGain]*fEperdn[0, wIbt, wGain]*fEperdn[0, wIbt, wGain],$
 fRead_noise_sq_e2[1, wIbt, wGain], ' +/-', $
 fSigma_sigmasquare[0, 1, wIbt, wGain]*fEperdn[1, wIbt, wGain]*fEperdn[1, wIbt, wGain]

print
print, format='(A, 2G20.3)', 'Read Noise(DN)     =', $
       sqrt(fRead_noise_sq_dn2[*, wIbt, wGain])
print, format='(A, 2G20.3)', 'Read Noise(e-)     =', $
       sqrt(fRead_noise_sq_e2[*, wIbt, wGain])

endfor
endfor


;******************* Write to text file ***********************
wDotpos=strpos(flnam, '.');
txtnam=strmid(flnam, 0, wDotpos) + '.txt'
openw, 101, txtnam

for wGain=0, wNgain-1 do begin
for wIbt=0, wNbt-1 do begin
printf, 101
printf, 101
printf, 101, flnam+gainstring[wGain]
printf, 101, '                          Left                     Right'
printf, 101, format='(A, 2(G10.3, A, G10.3))', $
       'For T_Int=0, ADC   =', fAdcoff[0, wIbt, wGain], ' +/-', $
                            fSigma_sigplusoff[0, 0, wIbt, wGain], $
                            fAdcoff[1, wIbt, wGain], ' +/-', $
                            fSigma_sigplusoff[0, 1, wIbt, wGain]
printf, 101, format='(A, 2(G10.3, A, G10.3))', 'Responsivity(e-/DN)=', $
                       fEperdn[0, wIbt, wGain], ' +/-', $
                       fSigma_sigplusoff[1, 0, wIbt, wGain], $
                       fEperdn[1, wIbt, wGain], ' +/-', $
                       fSigma_sigplusoff[1, 1, wIbt, wGain]
printf, 101

fCurrent_e=fCurrent_dn*fEperdn
printf, 101, format='(A, 2G24.3)', 'Signal (DN/s)      =', fCurrent_dn[*, wIbt, wGain]
printf, 101, format='(A, 2G24.3)', 'Signal (e-/s)      =', fCurrent_e[*, wIbt, wGain]
printf, 101, format='(A, 2G24.3)', 'Signal (nA cm-2)   =', fCurrent_e[*, wIbt, wGain] $
                                   *1.602e-19*1e9/(16.0e-4*16.0e-4)

printf, 101
printf, 101, format='(A, 2(G10.3, A, G10.3))', 'Rd Noise Sq(DN^2)  =', $
           fRead_noise_sq_dn2[0, wIbt, wGain], ' +/-', $
           fSigma_sigmasquare[0, 0, wIbt, wGain], $
           fRead_noise_sq_dn2[1, wIbt, wGain], ' +/-', $
           fSigma_sigmasquare[0, 1, wIbt, wGain]


fRead_noise_sq_e2=fRead_noise_sq_dn2*fEperdn*fEperdn
printf, 101, format='(A, 2(G10.3, A, G10.3))', 'Rd Noise Sq(e-^2)  =', $
 fRead_noise_sq_e2[0, wIbt, wGain], ' +/-', $
 fSigma_sigmasquare[0, 0, wIbt, wGain]*fEperdn[0, wIbt, wGain]*fEperdn[0, wIbt, wGain],$
 fRead_noise_sq_e2[1, wIbt, wGain], ' +/-', $
 fSigma_sigmasquare[0, 1, wIbt, wGain]*fEperdn[1, wIbt, wGain]*fEperdn[1, wIbt, wGain]

printf, 101
printf, 101, format='(A, 2G20.3)', 'Read Noise(DN)     =', $
       sqrt(fRead_noise_sq_dn2[*, wIbt, wGain])
printf, 101, format='(A, 2G20.3)', 'Read Noise(e-)     =', $
       sqrt(fRead_noise_sq_e2[*, wIbt, wGain])

endfor
endfor
close, 101


;******************* Write to text file for test report *******************
txlnam=strmid(flnam, 0, wDotpos) + '.txl'
openw, 101, txlnam
printf, 101, flnam, parmstring

for wIbt=0, wNbt-1 do begin
for wIlr=0, wNlr-1 do begin
for wGain=0, wNgain-1 do begin
printf, 101, portstring[wIlr, wIbt] + '   ' + gainstring[wGain]
printf, 101, $
  '       T_Int      Signal    Residual     Noise^2    Residual of Noise^2'

for wI=0, wNtint-1 do begin
  printf, 101, format='(5G12.4)', fTint[wI, wGain], $
            dSig[wI, wIlr, wIbt, wGain], dSigplusresiduals[wI, wIlr, wIbt, wGain], $
            dSigmasquare[wI, wIlr, wIbt, wGain], dSigmaresiduals[wI, wIlr, wIbt, wGain]
endfor

  printf, 101, format='(A, G10.3, A, G10.3)', $
       'For T_Int=0, ADC   =', fAdcoff[wIlr, wIbt, wGain], ' +/-', $
                            fSigma_sigplusoff[0, wIlr, wIbt, wGain]
  printf, 101, format='(A, G10.3, A, G10.3)', 'Responsivity(e-/DN)=', $
                       fEperdn[wIlr, wIbt, wGain], ' +/-', $
                       fSigma_sigplusoff[1, wIlr, wIbt, wGain]
  printf, 101, format='(A, G10.3)', 'Read Noise(DN)     =', $
       sqrt(fRead_noise_sq_dn2[wIlr, wIbt, wGain])
  printf, 101, format='(A, 2G20.3)', 'Read Noise(e-)     =', $
       sqrt(fRead_noise_sq_e2[wIlr, wIbt, wGain])
  printf, 101
endfor
endfor
endfor

close, 101


;************* Postscript plot **********************************
set_plot, 'ps'
psnam=strmid(flnam, 0, wDotpos) + '.ps'
device, filename=psnam, yoffset=0.5, ysize=10, /inches
!p.multi=[0, 2, 4]


for wIbt=0, wNbt-1 do begin
for wIlr=0, wNlr-1 do begin
for wGain=0, wNgain-1 do begin
  plot, fTint[*, wGain], dSigplusoff[*, wIlr, wIbt, wGain], $
        psym=-1, charsize=1.4, xtitle='T_Int', ytitle='ADC o/p(DN)', $
        title=flnam+portstring[wIlr, wIbt]+gainstring[wGain]

  !p.multi[0]=!p.multi[0]-1
  plot, fTint[*, wGain], dSigplusresiduals[*, wIlr, wIbt, wGain], $
        psym=-1, charsize=1.4, xtitle='T_Int', ytitle='Residuals(DN)'

  !p.multi[0]=!p.multi[0]+2
  plot, dSig[*, wIlr, wIbt, wGain], dSigmasquare[*, wIlr, wIbt, wGain], $
        psym=-1, charsize=1.4, xtitle='Signal(DN)', ytitle='Sigma^2(DN^2)' ;, $
        ;title=gainstring[wGain]
  oplot, dSig[*, wIlr, wIbt, wGain], fYfit_sigmasquare[*, wIlr, wIbt, wGain]

  !p.multi[0]=!p.multi[0]-1
  plot, dSig[*, wIlr, wIbt, wGain], dSigmaresiduals[*, wIlr, wIbt, wGain], $
        psym=-1, charsize=1.4, xtitle='Signal(DN)', ytitle='Residuals(DN^2)'
endfor
endfor
endfor
device, /close


if strupcase(!version.os_family) EQ 'WINDOWS' then $
  set_plot, 'win' $
 else set_plot, 'x'

!p.multi=[0, 1, 1]

print, 'Hit a key to delete the windows'
dummy1=get_kbrd(1)

wdelete
wdelete
wdelete
wdelete
wdelete
wdelete
wdelete
wdelete
if wHigainflag NE 0 then begin
  wdelete
  wdelete
  wdelete
  wdelete
  wdelete
  wdelete
  wdelete
  wdelete
endif

end