# calculus.test --
#    Test cases for the Calculus package
#
package require tcltest
source [file join [pwd] [file dirname [info script]] calculus.tcl]

puts "math::calculus [package present math::calculus]"


package require log
log::lvSuppress notice


namespace eval ::math::calculus::test {

namespace import ::tcltest::test
namespace import ::tcltest::cleanupTests
namespace import ::math::calculus::*

#
# Simple test functions - exact result predictable!
#
proc const_func { x } {
   return 1
}
proc linear_func { x } {
   return $x
}
proc downward_linear { x } {
   return [expr {100.0-$x}]
}

#
# Test the Integral proc
#
test "Integral-1.0" "Integral of constant function" {
   integral 0 100 100 const_func
} 100.0

test "Integral-1.1" "Integral of linear function" {
   integral 0 100 100 linear_func
} 5000.0

test "Integral-1.2" "Integral of downward linear function" {
   integral 0 100 100 downward_linear
} 5000.0

test "Integral-1.3" "Integral of expression" {
   integralExpr 0 100 100 {100.0-$x}
} 5000.0


proc const_func2d { x y } {
   return 1
}
proc linear_func2d { x y } {
   return $x
}

test "Integral2D-1.0" "Integral of constant 2D function" {
   integral2D { 0 100 10 } { 0 50 1 } const_func2d
} 5000.0
test "Integral2D-1.1" "Integral of constant 2D function (different step)" {
   integral2D { 0 100 1 } { 0 50 1 } const_func2d
} 5000.0
test "Integral2D-1.2" "Integral of linear 2D function" {
   integral2D { 0 100 10 } { 0 50 1 } linear_func2d
} 250000.0

#
# Note:
# Test cases for integral3D are missing!
#

#
# Test cases: yet to be brought into the tcltest form!
#

# xvec should one long!
proc const_func { t xvec } { return 1.0 }

# xvec should be two long!
proc dampened_oscillator { t xvec } {
   set x  [lindex $xvec 0]
   set x1 [lindex $xvec 1]
   return [list $x1 [expr {-$x1-$x}]]
}

foreach method {eulerStep heunStep rungeKuttaStep} {
   log::log notice "Method: $method"

   set xvec   0.0
   set t      0.0
   set tstep  1.0
   for { set i 0 } { $i < 10 } { incr i } {
      set result [$method $t $tstep $xvec const_func]
      log::log notice "Result ($t): $result"
      set t      [expr {$t+$tstep}]
      set xvec   $result
   }

   set xvec   { 1.0 0.0 }
   set t      0.0
   set tstep  0.1
   for { set i 0 } { $i < 20 } { incr i } {
      set result [$method $t $tstep $xvec dampened_oscillator]
      log::log notice "Result ($t): $result"
      set t      [expr {$t+$tstep}]
      set xvec   $result
   }
}

#
# Boundary value problems:
#
proc coeffs { x } { return {1.0 0.0 0.0} }
proc forces { x } { return 0.0 }

log::log notice [boundaryValueSecondOrder coeffs forces {0.0 1.0} {100.0 0.0} 10]
log::log notice [boundaryValueSecondOrder coeffs forces {0.0 0.0} {100.0 1.0} 10]

#
# Determining the root of an equation
# use simple functions
#
proc func  { x } { expr {$x*$x-1.0} }
proc deriv { x } { expr {2.0*$x} }

test "NewtonRaphson-1.0" "Result should be 1" {
   set result [newtonRaphson func deriv 2.0]
   if { abs($result-1.0) < 0.0001 } {
      set answer 1
   }
} 1
test "NewtonRaphson-1.1" "Result should be -1" {
   set result [newtonRaphson func deriv -0.5]
   if { abs($result+1.0) < 0.0001 } {
      set answer 1
   }
} 1

proc func2  { x } { expr {$x*exp($x)-1.0} }
proc deriv2 { x } { expr {exp($x)+$x*exp($x)} }

test "NewtonRaphson-2.1" "Result should be nearly 0.56714" {
   set result [newtonRaphson func2 deriv2 2.0]
   if { abs($result-0.56714) < 0.0001 } {
      set answer 1
   }
} 1

test "NewtonRaphson-2.1" "Result should be nearly 0.56714" {
   set result [newtonRaphson func2 deriv2 -0.5]
   if { abs($result-0.56714) < 0.0001 } {
      set answer 1
   }
} 1

cleanupTests
}

namespace delete ::math::calculus::test