PRO ucles_merge, nw, no, wb,wi, fn,fs, wm,fm, thresh


  ; Procedure to merge orders from a  2d echelle frame
  ; into a single spectrum.

  ; Input
  ; nw            n wavelength / order
  ; no            n orders
  ; wb   nw x no  wavelengths
  ; wi        no  number of valid wavelengths in ech order
  ; fs   nw x no  resampled rectification function
  ; fn   nw x no  resampled normalized spectrum

  ; thresh        cutoff low fluxes at order ends

  ; Output
  ; wm            merged wavelength grid
  ; fm            merged fluxes


  PRINT,'UCLES_MERGE: Merging echelle orders...'
      colors = GetColor(/Load)

  no = n_elements(wb)/nw

  ;; do not use data where local mean is small
  ;; fs is now the weighting function to be used in
  ;; the order merging process. 

  FOR j = 0,no-1 DO BEGIN
   sf = fs(0:wi(j),j)
   bar = MAX ( sf ) 
   msk = WHERE ( sf LT thresh * bar, cnt )
   IF ( cnt GT 0 ) THEN fs(msk,j) = 0.
  ENDFOR

  window,1
    shade_surf,fs,zrange=[0,max(fs)]
  window,2


  ; Find overlaps and merge

  wm = reform(wb(0:wi(no-1),no-1))
  fm = reform(fn(0:wi(no-1),no-1))


    !y.range=[0,1.5]
    FOR j = no-2,0,-1 DO BEGIN

    im = n_elements(wm) - 1
      plot,wm,fm,xrange=[wb(0,j+1),wb(wi(j),j)]
    ;    print,wm(0),wm(im),im
    ;   print,wb(0,j),wb(wi(j),j),j


    ilm = 0L
    WHILE ( wm(ilm) LT wb(0,j)+0.001 AND ilm LT im ) DO ilm = ilm + 1
    ium = im

    ilj = 0
    iuj = im - ilm
;    WHILE ( wb(iuj,j) GT wm(im)-0.001 ) DO iuj = iuj - 1

    k = j+1
    ilk = wi(k) - (ium-ilm)
    iuk = wi(k)

    IF (ium-ilm NE iuj-ilj OR ium-ilm NE iuk-ilk) THEN BEGIN
    print,'overlap error',ilj,iuj-ilj,ilm,ium-ilm,iuk,iuk-ilk
    print,wm(0),wm(ilm),wm(ium)
    print,0.,   wb(ilk,k),wb(iuk,k),wb(wi(k),k)
    print,0.,   wb(ilj,j),wb(iuj,j),wb(wi(j),j)
    STOP
    ENDIF

  ; check that overlap actually exsists before plotting it
    IF ( ilj NE 0 OR iuj NE 0 ) THEN BEGIN

    wdge = ( 0.5 + 0.5 * cos ( !pi * findgen(ium-ilm+1) / (ium-ilm) ) ) ^ 2


    oplot,wb(ilj+1:iuj,j),fn(ilj:iuj,j),color=colors.red
    oplot,wb(ilj+1:iuj,j),fs(ilj:iuj,j),color=colors.blue
    oplot,wb(ilj+1:iuj,j),(1.-wdge),color=colors.green
    oplot,wb(iuj+1:wi(j),j),fn(iuj+1:wi(j),j),color=colors.yellow
       ;    read,i
    wait,0.3

    fm(ilm:ium) = ( fm(ilm:ium) * fs(ilk:iuk,k)^2 * wdge  $
                  + fn(ilj:iuj,j) * fs(ilj:iuj,j)^2 * (1.-wdge) ) $
                  / ( fs(ilk:iuk,k)^2  * wdge + fs(ilj:iuj,j)^2 * (1.-wdge) )

    ENDIF

    wm = [wm , wb(iuj+1:wi(j),j)]
    fm = [fm , fn(iuj+1:wi(j),j)]


  ENDFOR
      
  END