;+
; kastbias.pro
;
; Subtract overscan region from Kast images using either one or
; two amplifiers.
;
; Usage:
;   kastbias, InFileName, OutFileName, legendre=legendre
;
; Author:  Elinor L. Gates, Lick Observatory, September 2008
;         Version 2.1 E. Gates Mar 2012 
;
; Version: 2.1
;   Extracts data regions and overscan regions from Kast FITS data.
;   Averages or medians the overscan region(s).
;   Optionally fits Legendre polynomial to overscan region(s).  
;   Subtracts overscan fit from data.  
;   Writes overscan subtracted data to a new FITS file. 
;
; Inputs:
;   InFileName - Input file name
;   OutFileName - Output file name
;   Fit - [optional] no or legendre fitting to overscan. 
;        Default is no fitting. 
;   Median - [optional] Median or average the overscan columns before
;            fitting.  Default is averaging.
;         
;
; Example:
;   kastbias,'r1000.fits','r1000new.fits'
;   kastbias,'r100.fits','r100new.fits',/legendre
;   kastbias,'r10.fits','r10os.fits',/median
;-

pro kastbias,infile,outfile,legendre=legendre,median=median

; read in FITS file and header
im=readfits(infile,h)

; read necessary parameters from FITS header
xsize=sxpar(h,'NAXIS1') ; number of columns (including overscan)
ysize=sxpar(h,'NAXIS2') ; number of rows
xorig=sxpar(h,'CRVAL1') ; readout origin on unbinned chip
yorig=sxpar(h,'CRVAL2') ; readout origin on unbinned chip
cdelt1=sxpar(h,'CDELT1') ; column change per pixel (i.e. binning and direction)
cdelt2=sxpar(h,'CDELT2') ; row change per pixel (i.e. binning and direction)
rover=sxpar(h,'ROVER') ; number of overscan rows for each amplifier
cover=sxpar(h,'COVER') ; number of overscan columns for each amplifier
inxsize=sxpar(h,'DNAXIS1') ; columns in sensor
inysize=sxpar(h,'DNAXIS2') ; rows in sensor
ampsx=sxpar(h,'AMPSCOL') ; amplifiers per column
ampsy=sxpar(h,'AMPSROW') ; amplifiers per row

; determine number and sizes of overscan and data regions
namps=ampsx*ampsy
if (rover gt 0) then begin 
 over=rover 
 print,'Program does not yet handle row overscans'
 exit
endif else over=cover
if (over eq 0) then begin
    print,'No overscan region. Exitting.'
    exit                   
endif    

; Single Amplifier Mode
if (namps eq 1) then begin

  ; extract overscan regions (each over scan region is cover pixels wide)
  biassec=float(im[xsize-cover:xsize-1,*])

  ; average or median overscan regions along the readout axis
  n=(size(biassec))[2]
  overscan=fltarr(n)
  if keyword_set(median) then $
    for i=0,n-1 do overscan[i]=median(biassec[*,i]) else $
    for i=0,n-1 do overscan[i]=total(biassec[*,i])/float(cover)

  ; do fit requested
  x=findgen(n)
 if keyword_set(legendre) then begin
  ; fit third order legendre polynomial to overscan region
   ofit=svdfit(x,overscan,4,yfit=overscanfit,/legendre,/double)
endif else overscanfit=overscan
 
  ; extract data regions and subtract overscan fits
  imnew=float(im[0:xsize-1-cover,*])
  m=(size(imnew))[1]
  for i=0,m-1 do imnew[i,*]=imnew[i,*]-overscanfit 

endif

; Two Amplifier Mode
if (namps eq 2) then begin

  ;extract overscan regions (each over scan region is cover pixels wide)
  biasseca=float(im[xsize-cover*2:xsize-1-cover,*])
  biassecb=float(im[xsize-cover:xsize-1,*])

  ;average or median overscan regions along the readout axis
  na=(size(biasseca))[2] 
  nb=(size(biassecb))[2]
  oa=fltarr(na)
  ob=fltarr(nb)
  if keyword_set(median) then begin
    for i=0,na-1 do oa[i]=median(biasseca[*,i])
    for i=0,nb-1 do ob[i]=median(biassecb[*,i])
  endif else begin
    for i=0,na-1 do oa[i]=total(biasseca[*,i])/float(cover)
    for i=0,nb-1 do ob[i]=total(biassecb[*,i])/float(cover)
  endelse

 if keyword_set(legendre) then begin
  ;fit third order legendre polynomial to overscan regions
  xa=findgen(na)
  xb=findgen(nb)
  oafit=svdfit(xa,oa,4,yfit=oaf,/double,/legendre)
  obfit=svdfit(xb,ob,4,yfit=obf,/double,/legendre)
endif else begin 
  oaf=oa
  obf=ob
endelse

  ;extract data regions and subtract overscan fits
  hsize=abs(inxsize/cdelt1)
  ;calc x origin of readout in binned units if cdelt1 negative or positive 
  if (cdelt1 lt 0) then xorig=(xorig-(xsize-2*cover)*abs(cdelt1))/abs(cdelt1) $
    else xorig=xorig/cdelt1
  x0=xorig+xsize-1-cover*2

  if (x0 lt hsize/2) then begin ; all data on Left amplifier
   imnew=float(im[0:xsize-1-cover*2,*])
   m=(size(imnew))[1]
   for i=0,m-1 do imnew[i,*]=imnew[i,*]-oaf 
  endif

  if (xorig ge hsize/2) then begin ; all data on Right amplifier
   imnew=float(im[0:xsize-1-cover*2,*])
   m=(size(imnew))[1]
   for i=0,m-1 do imnew[i,*]=imnew[i,*]-obf 
  endif

  if ( (xorig lt hsize/2) and (x0 gt hsize/2) ) then begin ; data on both amps
   x1=hsize/2-xorig
   ima=float(im[0:x1-1,*])
   imb=float(im[x1:xsize-1-cover*2,*])
   ma=(size(ima))[1]
   mb=(size(imb))[1]

   for i=0,ma-1 do ima[i,*]=ima[i,*]-oaf
   for i=0,mb-1 do imb[i,*]=imb[i,*]-obf

   ;rejoin the two subimages into a single image
   imnew=fltarr(xsize-cover*2,ysize)
   imnew[0:x1-1,*]=ima
   imnew[x1:*,*]=imb
  endif 

endif

if (namps gt 2) then begin
    print,'More than two amplifiers not yet supported. Exiting.'
endif

; update header information
bzero=sxpar(h,'BZERO')
bzero=float(bzero)
if bzero ne 0 then begin
  sxaddpar,h,'O_BZERO',bzero,'Original Data is unsigned Integer'
  sxaddpar,h,'BZERO',0.0,/PDU
endif

bscale=sxpar(h,'BSCALE')
bscale=float(bscale)
if bscale ne 1.0 then begin
  sxaddpar,h,'O_BSCALE',bscale,'Original Data is scaled'
  sxaddpar,h,'BSCALE',1.0,/PDU
endif

sxaddpar,h,'HISTORY',"= 'OVERSCAN SUBTRACTED'"

; write overscan subtracted image to file
writefits,outfile,imnew,h

end