pro gs_cont,xval,cont,order
;+
;  Routine to interactively define the continuum
;
;
;  Input:
;          xval     array of x-values at which continuum will be calculated
;          order    order of polynomial to fit (optional - prompted for)
;
;  Output:
;          cont     output array - the calculated continuum at the x-values
;                                  specified by xval
;
;  Use:
;       Plot a spectrum and then..
;       IDL> gs_cont,xval, cont [,order]
;
;
;  CDP   Feb 92
;-

xc = fltarr(100)
yc = fltarr(100)

print,'Use cursor to locate continuum points...'
print,'Input same point twice to end.'
xlast = 0.0
ylast = 0.0
xpts = 0
for i=0,99 do begin
   cursor,x1,y1,3,/data   
   oplot,[x1,x1],[y1,y1],psym=6
   if xlast eq x1 and ylast eq y1 then i=99
   xlast = x1
   ylast = y1
   xc(i) = x1
   yc(i) = y1
   xpts = xpts + 1
endfor

xpts = xpts - 1
xc = xc(0:xpts-1)
yc = yc(0:xpts-1)

ans = 'N'
cont = fltarr(n_elements(xval))


repeat begin

   if n_params() lt 3 then begin
      order = 0
      while (order lt 1 or order gt 4) do read,$
                          'What order of polynomial to fit? (max 4)> ',order
   endif

   coeff = poly_fit(xc,yc,order,dum)
   tc = fltarr(5)
   tc(*) = 0.0
   tc(0:order) = coeff
   
   for i=0,n_elements(xval)-1 do begin
      x = xval(i)
      cont(i) = x*(x*(x*(x*tc(4)+tc(3))+tc(2))+tc(1))+tc(0) 
   endfor

   oplot,xval,cont
   read,'Try again? >',ans

endrep until strupcase(ans) eq 'N'

end