PRO extract1d_ACS,specname,centery,numpix,outname,PYES=pyes,OPTEX=exopt

; NAME:
;       EXTRACT1D_ACS
; PURPOSE:
;       Extract a 1D spectrum from a 2D ACS slitless grism image in a given aperture. 
; EXPLANATION:
;       Program to be used together with reduced 2D ACS slitless grism spectra, in 
;       order to extract 1D spectra in a rectangular or optimal extraction aperture 
;       with any given size. Optimal extraction first fits the profile of the 
;       spectrum in the direction perpendicular to the dispersion and applies this
;       profile to perform optimal extraction. No aperture correction is applied in
;       either extraction, i.e. the result will be a lower estimate of the flux. 
;
; CALLING SEQUENCE:
;       EXTRACT1D_ACS, specname, centery, numpix, outname, [PYES], [OPTEX]
;
; INPUTS:
;       specname - Name (string) of 2D grism fits file containing the spectrum.
;       centery  - Y pixel to extract centred spectrum above and below.
;       numpix   - Size of box to extract from, in HST pixels.
;       outname  - Name (string) of output spectrum in ASCII format, with the 
;                  first column being wavelength in Angstroms, second column
;                  the extracted spectra in f_lambda (erg/s/cm2/AA), third column
;                  the error on the spectrum in f_lambda and fourth column
;                  the contamination in f_lambda.
;
; OPTIONAL INPUT KEYWORDS:
;       PYES     - Set to 1 if the extracted spectrum should be plotted.
;       OPTEX    - Set to 1 if optimal extraction is wanted.
;
; OUTPUTS:
;       Extracted spectrum is written to the 'outname' ASCII file.
;
; EXAMPLES:
;       2D grism spectrum is stored in 'UDFNICP1_spec901A.2Dstamp.fits' and we wish
;       to extract a new spectrum around central y-pixel 61, in a box of 3 pixels.
;       The result will be stored in the ASCII file 'test.tab' and it will also be 
;       plotted on the screen.
;
;       IDL> EXTRACT1D_ACS,'UDFNICP1_spec901A.2Dstamp.fits',61,3,'test.tab',PYES=yes
;
; RESTRICTIONS:
;       The extraction is performed symmetrically around the central y pixel chosen,
;       i.e. if the central y position is close to the edge of the image this will
;       contrain the extraction box which may miss flux at the edges.
;
; PROCEDURES CALLED
;       MRDFITS(), READCOL(), SXPAR(), WRTAB_EXP()
;
; REVISION HISTORY:
;       Written        K.K. Nilsson                    May 2010
on_error,2
!EXCEPT=0
angstrom = '!6!sA!r!u!9 %!6!n'

if N_params() lt 4 then begin
  print,'Syntax - EXTRACT1D, specname, centery, numpix, outname, [PYES], [OPTEX]'
  return
endif

if size(specname,/tname) NE 'STRING' then $
  message,'ERROR - Supplied SPECNAME keyword must be a scalar string'
if size(outname,/tname) NE 'STRING' then $
  message,'ERROR - Supplied OUTNAME keyword must be a scalar string'

!P.charsize=1.5
!P.charthick=2

centery=centery-1
im=mrdfits(specname,1,hdr,/silent)
wl0=sxpar(hdr,'CRVAL1')
wl0x=sxpar(hdr,'CRPIX1')
wlstep=sxpar(hdr,'CDELT1')
wl=findgen(n_elements(im(*,0)))*wlstep+wl0-(wl0x-1)*wlstep
wlind=where(wl ge 4300 and wl le 10700)
wl=wl(wlind)
majax=sxpar(hdr,'OBJSHAPA')
minax=sxpar(hdr,'OBJSHAPB')
objang=sxpar(hdr,'OBJPA')
slitwi=2.35/sqrt((cos(objang*!pi/180)^2/majax^2+sin(objang*!pi/180)^2/minax^2))

readcol,'ACS.WFC.1st.sens.7.dat',wlsens,sens,format='D,D,X',/silent
newx2=findgen(301)-150.
newres=6563./(66.0*0.11/slitwi)/40.
newgauss=double(exp(-newx2^2/(2*(newres/(2*sqrt(2.*alog(2.))))^2))/total(exp(-newx2^2/(2*(newres/(2*sqrt(2.*alog(2.))))^2))))
newflux=convol(double(sens),newgauss,/edge_truncate)
sens2=interpol(newflux,wlsens,wl)

newtab=dblarr(4,n_elements(wlind))
newtab(0,*)=wl
extnum=[1,2,4]
for o=0,2 do begin

im=mrdfits(specname,extnum(o),hdr,/silent)
sumfrom=max([0,centery-numpix/2])
sumto=min([n_elements(im(0,*))-1,centery+numpix/2])
if o eq 0 then begin

prof=dblarr(n_elements(im(0,*)))
for j=0,n_elements(prof)-1 do begin
  prof(j)=total(im(*,j))/n_elements(im(*,j))
endfor

endif


if o eq 0 or o eq 2 then begin 

onedspec=dblarr(n_elements(im(*,0)))
for n=0,n_elements(im(*,0))-1 do begin
   onedspec(n)=total(im(n,sumfrom:sumto))
endfor

if keyword_set(exopt) then begin
  fluxnorm=total(onedspec)

  onedspec=dblarr(n_elements(im(*,0)))
  for n=0,n_elements(im(*,0))-1 do begin
    onedspec(n)=total(prof(sumfrom:sumto)*im(n,sumfrom:sumto))
  endfor

  onedspec=onedspec*fluxnorm/total(onedspec)

endif

endif

if o eq 1 then begin 

onedspec=dblarr(n_elements(im(*,0)))
for n=0,n_elements(im(*,0))-1 do begin
   onedspec(n)=total(im(n,sumfrom:sumto)^2)
endfor

onedspec=sqrt(onedspec)

if keyword_set(exopt) then begin

  fluxnorm=total(onedspec)

  onedspec=dblarr(n_elements(im(*,0)))
  for n=0,n_elements(im(*,0))-1 do begin
    onedspec(n)=total(prof(sumfrom:sumto)^2*im(n,sumfrom:sumto)^2)
  endfor

onedspec=sqrt(onedspec)*fluxnorm/total(sqrt(onedspec))

endif

endif

fluxy=onedspec(wlind)/sens2(wlind)/40.

if keyword_set(pyes) and o eq 0 then begin
ymin=max([0.9*min(fluxy(where(wl ge 5500 and wl le 9500))),1E-20])
ymax=1.1*max(fluxy(where(wl ge 5500 and wl le 9500)))
  plot,wl,fluxy,/ylog,xra=[5000,10000],/xsty,xtitle='Wavelength ('+angstrom+')',$
       ytitle='Flux (erg s!U-1!N cm!U-2!N '+angstrom+'!U-1!N)',thick=3,yra=[ymin,ymax],/ysty
fluxy0=fluxy
endif
if keyword_set(pyes) and o eq 1 then begin
  oploterror,wl,fluxy0,fluxy,thick=3
endif
if keyword_set(pyes) and o eq 2 then begin
  oplot,wl,fluxy,thick=3,linesty=1
endif

newtab(o+1,*)=fluxy

endfor

wrtab_exp,newtab,outname

end