PRO normalise, creg ;; convert flux spectra to normalized spectra


; COMMON fileblock
COMMON fitsblock
COMMON procblock

;; Input
;
;  creg = [w1l, wlu, w2l, w2u, ....]
;         pairs of wavelength defining regions for fitting the continuum


;; evaluate the median spectrum

specm=fltarr(nwave)
FOR iw=0,nwave-1 DO BEGIN 
specm(iw) = median(spec(iw,0:nspec-1))
ENDFOR


cont = fltarr(nwave,nspec)

;; define continuum regions

ncont = n_elements(creg)
creg_loc = where(wave ge creg(0) and wave le creg(1))

FOR ic = 0,ncont-1,2 DO $
  creg_loc =[creg_loc,where(wave ge creg(ic) and wave le creg(ic+1))]


;; create 2d contiuum fit


cforder = 3                   ;; try a cubic fit for continuum

FOR is = 0,nspec-1 DO BEGIN

;; fit 1d polynomial to continuum

   poly = poly_fit(wave(creg_loc)-wave(0), spec(creg_loc,is), cforder)

;; evaluate the continuum polynomial

   cont(*,is) = poly(cforder)
   FOR ip = cforder-1,0,-1 DO BEGIN
      cont(*,is) = cont(*,is)*(wave-wave(0)) + poly(ip)
   ENDFOR

ENDFOR


;; divide 2d spectrum by 2d continuum

spec = spec / cont


;; find continuum for median spectrum


cont = fltarr(nwave)

   poly = poly_fit(wave(creg_loc)-wave(0), specm(creg_loc), cforder)
   cont = poly(cforder)
   FOR ip = cforder-1,0,-1 DO BEGIN
      cont = cont*(wave-wave(0)) + poly(ip)
   ENDFOR


window,8,title="normalise.pro"
p_save = !p.multi
!p.multi=[0,0,2,0,0]
plot,wave,specm
oplot,wave,cont

  ;; normalise median spectrum

specm = specm/cont

plot,wave,specm
!p.multi=p_save


END