;+
; pfcambias.pro
;
;  Subtract overscan region from PFCam images in two-amp mode
;
;  Usage:
;    Pfcambias, InFileName, OutFileName
;
;  Author:  Elinor L. Gates, Lick Observatory, April 2000
;
;  Version: 3
;    Adapted from IRAF procedure pfcambias.cl by Mike Bolte.
;    Extracts data regions and overscan regions from PFCam FITS data.  
;    Fits Legendre polynomial to overscan region. Subtracts overscan 
;    fit from data.  Writes overscan subtracted data to new FITS file.
;    Updated to work with Fairchild 4096x4096 array and UCAM controller
;    2006.
;
;  Example:
;    pfcambias,'d1000.ccd','d1000new.fits'
;-

pro pfcambias,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)
biasseca=float(im[xsize-cover*2:xsize-1-cover,*])
biassecb=float(im[xsize-cover:xsize-1,*])

;average overscan regions along the readout axis
na=(size(biasseca))[2]
nb=(size(biassecb))[2]
oa=fltarr(na)
ob=fltarr(nb)
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)

;fit third order legendre polynomial to overscan regions
xa=findgen(na)
xb=findgen(nb)
oafit=svdfit(xa,oa,2,yfit=oaf,/legendre)
obfit=svdfit(xb,ob,2,yfit=obf,/legendre)

;extract data regions and subtract overscan fits
hsize=abs(4096/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 

;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