pro grismspec,phafile,tab,tab1,lower=lower,upper=upper, dshift = dshift, $
     mkwave = mkwave,coeff=ccoeff,srcwid=srcwid,centroid=centroid, $
      wcoeff=wcoeff, minbkg=minbkg
;+
; NAME:
;     GRISMSPEC
; PURPOSE:
;     Extract spectrum from a .pha file created by UVOTIMGRISM into an IDL
;     structure.     
;
;     By default, performs a robust polynomial fit to the background, 
;     Has options to allow only the lower or upper background to be extracted.
;      The system variable CALDB must be defined to point to 
;     the CALDB directory    
; CALLING SEQUENCE:
;     GRISMSPEC, filename, str, [orig_str] /Lower,/Upper, dshift='
; REQUIRED INPUTS:
;     filename - scalar string giving name of the .pha file.   The
;        .pha extension is optional
; OUTPUTS: 
;     str - IDL structure with the following self-explanatory tags
;      IDL> help,/str,tab
;        W               FLOAT     Array[910]
;        DIS             FLOAT     Array[910]
;        GROSS           DOUBLE    Array[910]
;        BKG1            DOUBLE    Array[910]
;        BKG2            DOUBLE    Array[910]
;        NET             DOUBLE    Array[910]
;        FLUX            DOUBLE    Array[910]
;        ERR             DOUBLE    Array[910]
; 
;     The .DIS tag gives the pixel number of each wawvelength.    
;     GRISMSPEC does its own customized extraction using a robust polynomial
;     to fit the background
;
;     orig_str - IDL structure giving the default UVOTIMGRISM extraction 
;            spectrum
;       IDL> help,/str,tab1
;           W               FLOAT     Array[910]
;           BKG             INT       Array[910]
;           NET             INT       Array[910]
;           FLUX            FLOAT     Array[910]

; OPTIONAL KEYWORD INPUTS:
;   /Lower - Use only the lower background (default is to average upper and 
;            lower)
;   /Upper - Use only the upper background
;   dshift - manual shift in pixels of the centroid along the dispersion 
;            direction.   Useful when true centroid cannot be determined because
;            of saturation.     Note that since the wavelength is nonlinear 
;            with pixels 
; EXAMPLE:
;     Compare the default and grismspec spectrum of a file w00055502010ugu_1.pha
;
;        IDL> grismspec,'sw00055502010ugu_1',tab,tab1
;        IDL> plot,tab.w,tab.flux,/xsty       ;Grismspec spectrum
;        IDL> oplot,tab1.w,tab1.flux,/lines   ;Original extraction
;
; REVISION HISTORY: 
;     Written W. Landsman          August 2005
;     Update require for new FITS_READ which no longer automatically inherits
;     the primary header  W. Landsman      April 2007
;     Don't assume spectrum is centered    Dec 2007
;-
if N_params() LT 2 then begin
    print,'Syntax - GRISMSPEC, filename, str, [orig_str, /LOWER, /UPPER, 
    print,'                           dshift=, /centroid, srcwid=]'          
    return
endif

caldir = '$CALDB/data/swift/uvota/cpf/arf/'
coifac = !pi*5.*5./0.557/0.557


fdecomp,phafile,disk,dir,fname,ext
if ext EQ '' then ext = 'pha'
file = disk + dir + fname + '.' + ext
fits_open,file,fcb
fits_read,fcb,tab,htab,exten=2
fits_read,fcb,im,h,exten=3
fits_close,fcb

filter = strtrim(sxpar(h,'FILTER'),2)
clock = sxpar(h,'WHEELPOS')
bin = sxpar(h,'BINX')
etime = sxpar(h,'exposure')

if filter EQ 'UGRISM' then begin
    arf = caldir + 'swugu0' + string(clock,'(i3)') + '_20041120v101.arf'
   extname = 'SPECRESPUGRISM' + strtrim(clock,2)
endif else begin
    arf = caldir + 'swugv' + string(clock,'(I4.4)') + '_20041120v101.arf'
   extname = 'SPECRESPVGRISM' + strtrim(clock,2)
endelse


w = tbget(htab,tab,'LAMBDA')
f = tbget(htab,tab,'NET')
flux = tbget(htab,tab,'FLUX')
dim = size(im,/dimen)
npts_w = N_elements(w)
wn = dim[0] - npts_w
dis = findgen(npts_w) + wn - sxpar(h,'crpix1') -1
if N_elements(dshift) EQ 1 then w = dis2wave(fix(clock),dis + dshift, bin=bin)

if keyword_set(mkwave) then w = dis2wave(clock, dis,bin=bin) $ 
else if N_elements(wcoeff) GT 2 then w = poly(dis,wcoeff)
wp = wn 

if N_params() EQ 3 then begin
 fdecomp,file,disk,dir,fname
 tab1 = mrdfits(fname + '_back.pha',1,htab1,/sil)
 obkg = reverse(tab1.counts)
 tab1 = {w:w,bkg:obkg,net:f,flux:flux}
endif

hh = sxpar(fcb.hmain,'HISTORY')
tmp = gettok(hh,' ')
keyword = gettok(hh,'=')
keyword = strtrim(keyword,2) & hh = strtrim(hh,2)
g = where(keyword EQ 'sourcex')
sourcex = float(strcompress(hh[g[0]],/remove_all))
g = where(keyword EQ 'sourcey')
sourcey = strcompress(hh[g[0]],/remove_all)
g = where(keyword EQ 'infile')
infile = strcompress(hh[g[0]],/remove_all)
g = where(keyword EQ 'ang')
angle = strcompress(hh[g[0]],/remove_all)
if N_elements(srcwid) EQ 0 then begin 
   g = where(keyword EQ 'srcwid')
   srcwid = float(strcompress(hh[g[0]],/remove_all))
endif else srcwid=float(srcwid)
g = where(keyword EQ 'bkgwid1')
bkgwid1 = fix(strcompress(hh[g[0]],/remove_all))
g = where(keyword EQ 'bkgwid2')
bkgwid2 = fix(strcompress(hh[g[0]],/remove_all))
g = where(keyword EQ 'bkgoff1')
bkgoff1 = fix(strcompress(hh[g[0]],/remove_all))
g = where(keyword EQ 'bkgoff2')
bkgoff2 = fix(strcompress(hh[g[0]],/remove_all))

object = sxpar(h,'OBJECT') 
dateo = sxpar(h,'DATE-OBS')
expt = string(sxpar(h,'EXPOSURE'),f='(f6.1)') 
imx = im[wp:*,*]
imz = smooth(im,17)
imz = imz[wp:*,*]	   

wid2 = (srcwid-1)/2
cen = bkgwid1 + bkgoff1 + wid2
if keyword_set(centroid) then ggrism,phafile,ccoeff
if N_elements(ccoeff) GT 0 then begin 
    
     yy = round( poly(wp+findgen(npts_w), ccoeff)) > 0
     gross = fltarr(npts_w)
  
     y1 = (yy - wid2)>0
     y2 = (yy + wid2) < (dim[1]-1)
     fac = (2.*wid2+1)/(y2-y1+1)
     
     for i=0,npts_w-1 do $ 
           gross[i] = fac[i]*total(imx[i,y1[i]:y2[i]])
         
endif else begin 
     gross = total(imx[*,cen-wid2:cen+wid2],2)
     yy = replicate(cen,npts_w)
endelse  
 
    coicnt = fltarr(npts_w)
     xcoi1 = (wp + findgen(npts_w) -7) > 0
     xcoi2 = (xcoi1 + 14) < (dim[0]-1)
     ycoi1 = (yy -7) > 0
     ycoi2 = (ycoi1 + 14) < (dim[1]-1)
     pixarea = (xcoi2-xcoi1+1)*(ycoi2-ycoi1+1)
     for i=0,npts_w-1 do $
           coicnt[i] = total(im[xcoi1[i]:xcoi2[i], ycoi1[i]:ycoi2[i]])
     coicnt = coicnt*25*!pi/(pixarea*0.557*0.557)    ;Scale to a 5" circle	   
     coicnt = coicnt/etime     ;Scale to cps
 
 
b2= (yy- wid2 -bkgoff1-1)
b1 = b2-(bkgwid1-1) > 0
npts = b2-b1+1

bkg1 = yy*0.
for i=0,Npts_w-1 do bkg1[i] = total(imx[i,b1[i]:b2[i]],2)
bkg1 = bkg1*bkgwid1/float(npts)
xx = findgen(N_elements(w))
coeff = robust_poly_fit(xx,bkg1,5,yfit,sig)
bkg1 = yfit

b1 = yy + bkgoff2 +wid2 +1 < (dim[1]-1)
b2 = b1 + (bkgwid2-1) < (dim[1]-1)
npts = b2-b1+1
bkg2 = yy*0.
for i=0,Npts_w-1 do bkg2[i] = total(imx[i,b1[i]:b2[i]],2)
bkg2 = bkg2*bkgwid2/float(npts)
coeff = robust_poly_fit(xx,bkg2,5,yfit)
bkg2 = yfit


if keyword_set(lower) then begin 
   bkg = bkg1 
   bwid = bkgwid1
endif   else if keyword_set(upper) then begin
   bkg = bkg2 
   bwid = bkgwid2
endif else begin   
 if keyword_set(minbkg) then begin 
      bkg = bkg1 < bkg2
      bwid = bkgwid1
      endif else begin 
bkg = (bkg1 + bkg2)
bwid = bkgwid1 + bkgwid2
      endelse 
endelse

coinc_corr,coicnt, 0, sat, ratio=sratio
coibkg = bkg*coifac/bwid/etime
coinc_corr,coibkg, 0, sat, ratio=bratio
net = gross- bkg*srcwid/bwid
coicnt = gross*sratio- bkg*bratio*srcwid/bwid
ratio = coicnt/net
net = gross - bkg*srcwid/bwid
err = sqrt( gross )

fits_open,arf,fcb
;exten = where(fcb.extname EQ extname,Ng)
fits_close,fcb
exten = [1]
;if Ng EQ 0 then message,'Unrecognized extension ' + extname
tab = mrdfits(arf,exten[0],htab,/sil)
if tag_exist(tab,'wave_max') then $
  warf = (tab.wave_max + tab.wave_min)/2. else $
  warf = (tab.wave_lo + tab.wave_hi)/2.
 wdiff = w[1:*] - w
 wdiff = [wdiff,wdiff[N_elements(wdiff)-1] ]

quadterp,warf,tab.specresp,w,farea


flux = (net/farea/wdiff/etime)*1.9862e-8/w
err = (err/farea/wdiff/etime)*1.9862e-8/w

tab = {w:w,dis:dis,gross:gross,bkg1:bkg1,bkg2:bkg2,net:net,flux:flux,err:err, $
        y:yy,coicnt:ratio}
	return
end