FUNCTION multiply_curves, input1, input2, xrange=xrange, qstop=qstop

; MODIFICATION HISTORY
progver = 'v2007-May-15' ;--- (M.Weber) Written. (But based heavily on
;                             previous work.
progver = 'v2008-Oct-01' ;--- (M.Weber) Reviewed version.

  qstop = keyword_set(qstop)


; === Set x-range for product =========================


  x_lows   = [min(input1[*,0]), min(input2[*,0])]
  x_highs  = [max(input1[*,0]), max(input2[*,0])]

  if not keyword_set(xrange) then begin

    x_lo_out = min(x_lows)
    x_lo_in  = max(x_lows)
    x_hi_in  = min(x_highs)
    x_hi_out = max(x_highs)

    ss1 = where((input1[*,0] ge x_lo_in) and (input1[*,0] le x_hi_in), $
                count1)
    ss2 = where((input2[*,0] ge x_lo_in) and (input2[*,0] le x_hi_in), $
                count2)

    if (count1 ge count2) then xrange = input1[ss1,0] $
                          else xrange = input2[ss2,0]
  endif else begin

    ss1 = where((input1[*,0] ge min(xrange)) and (input1[*,0] le $
                max(xrange)), count1)
    ss2 = where((input2[*,0] ge min(xrange)) and (input2[*,0] le $
                max(xrange)), count2)
    if ((count1 eq 0) or (count2 eq 0)) then begin
      print, 'ERROR: The inputs do not both overlap the ' $
             + 'given x-range. --- Aborting.'
      return, 0
    endif
  endelse

  if (input1[ss1[0],0] gt x_lows[0])  then begin
    ss1 = [ss1[0]-1,ss1]
    if (input1[ss1[count1],0] lt x_highs[0]) then ss1 = [ss1,ss1[count1]+1]
  endif else begin
    if (input1[ss1[count1-1],0] lt x_highs[0]) then ss1 = [ss1,ss1[count1-1]+1]
  endelse
  if (input2[ss2[0],0] gt x_lows[1])  then begin
    ss2 = [ss2[0]-1,ss2]
    if (input2[ss2[count2],0] lt x_highs[1]) then ss2 = [ss2,ss2[count2]+1]
  endif else begin
    if (input2[ss2[count2-1],0] lt x_highs[1]) then ss2 = [ss2,ss2[count2-1]+1]
  endelse


; === Interpolate both curves onto the x-range ========

  input1b = interpol(input1[ss1,1], input1[ss1,0], xrange)
  ss = where((xrange lt x_lows[0]) or (xrange gt x_highs[0]), count)
  if (count gt 0) then input1b[ss] = 0.

  input2b = interpol(input2[ss2,1], input2[ss2,0], xrange)
  ss = where((xrange lt x_lows[1]) or (xrange gt x_highs[1]), count)
  if (count gt 0) then input2b[ss] = 0.


; === Multiply curves =================================

if qstop then stop

  output = input1b * input2b


; === Finish ==========================================

  return, output


END


; =========================================================================