# ::math::integrate --
 #
 # calculate the area under a curve defined by a
 # set of (x,y) data pairs.  the x data must increase
 # monotonically throughout the data set for the
 # calculation to be meaningful, therefore the
 # monotonic condition is tested, and an error is thrown
 # if the x value is found to be decreasing.
 #
 # Arguments:
 #    (x,y) pairs in list format
 # 
 # Results:
 #    Area and error bound
 #    (calculation is "Simpson's rule")
 #
 # at least 5 data pairs are required, and if the number
 # of data pairs is even, a padding value of (x0,0) will
 # be added.
 #
 # 04-13-2000 -- by: Phil Ehrens at The LIGO Lab
 
 namespace eval math {} 

 proc ::math::integrate { xy_pairs } {
      
      set length [ llength $xy_pairs ]
      
      if { $length < 10 } {
         return -code error "at least 5 x,y pairs must be given"
      }   
      
      ;## are we dealing with x,y pairs?
      if { [ expr $length % 2 ] } {
         return -code error "unmatched xy pair in input"
      }
      
      ;## are there an even number of pairs?  Augment.
      if { ! [ expr $length % 4 ] } {
         set xy_pairs [ concat [ lindex $xy_pairs 0 ] 0 $xy_pairs ]
      }
      set x0   [ lindex $xy_pairs 0     ]
      set x1   [ lindex $xy_pairs 2     ]
      set xn   [ lindex $xy_pairs end-1 ]
      set xnminus1 [ lindex $xy_pairs end-3 ]
     
      if { $x1 < $x0 } {
         return -code error "monotonicity broken by x1"
      }
 
      if { $xn < $xnminus1 } {
         return -code error "monotonicity broken by xn"
      }   
      
      ;## handle the assymetrical elements 0, n, and n-1.
      set sum [ expr {[ lindex $xy_pairs 1 ] + [ lindex $xy_pairs end ]} ]
      set sum [ expr {$sum + (4*[ lindex $xy_pairs end-2 ])} ]
 
      set data [ lrange $xy_pairs 2 end-4 ]
      
      set xmax $x1
      set i 1
      foreach {x1 y1 x2 y2} $data {
         incr i
         if { $x1 < $xmax } {
            return -code error "monotonicity broken by x$i"
         }
         set xmax $x1
         incr i
         if { $x2 < $xmax } {
            return -code error "monotonicity broken by x$i"
         }
         set xmax $x2
         set sum [ expr {$sum + (4*$y1) + (2*$y2)} ]
      }   
      
      if { $xmax > $xnminus1 } {
         return -code error "monotonicity broken by xn-1"
      }   
     
      set h [ expr { ( $xn - $x0 ) / $i } ]
      set area [ expr { ( $h / 3.0 ) * $sum } ]
      set err_bound  [ expr { ( ( $xn - $x0 ) / 180.0 ) * pow($h,4) * $xn } ]  
      return [ list $area $err_bound ]
 }