;+
; ccd2bias.pro
;
;  Subtract overscan region from CCD-2 images in one-amp mode
;
;  Usage:
;    ccd2bias, InFileName, OutFileName
;
;  Author:  Elinor L. Gates, Lick Observatory, December 2006
;
;  Version: 2
;    Adapted from IRAF procedure pfcambias.cl by Mike Bolte.
;    Extracts data regions and overscan regions from CCD-2 FITS data.  
;    Fits Legendre polynomial to overscan region. Subtracts overscan 
;    fit from data.  Writes overscan subtracted data to new FITS file.
;
;  Example:
;    ccd2bias,'d100.fits','d100_new.fits'
;-

pro ccd2bias,infile,outfile

;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
cdelt1=sxpar(h,'CDELT1') ; column change per pixel (i.e. binning and direction)
cover=sxpar(h,'COVER') ; number of overscan columns for each amplifier

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

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

;fit third order legendre polynomial to overscan regions
x=findgen(n)
ofit=svdfit(x,overscan,2,yfit=overscanfit,/legendre)

;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 

;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