package require Tk
  package require Plotchart
  console show

  # plotf --
  #     Convenience procedure to plot a function of one variable or
  # call a user defined proc with arguments: 
  #
  # Arguments:
  #     range       Range over which to plot the function (x-axis)
  #     expression  Expression by which the function is defined
  #                 (independent variable: t)
  #
  proc plotf {range expression args} {
    # intial is the initial conditions
    # nosteps - number of time steps to progress.
    # args are arguments to be interpreted; currently
    #  -title "XXX"
    # -nosteps Number of steps to use

    # Very simple for now: no parameters
    set nosteps 100
    foreach {name value} $args {
        switch -- $name {
            "-nosteps" {set nosteps $value}
            "-initial" {set x $value}
            default {}
        }
    }
    
    foreach {tmin tmax} $range {break}
    set tstep [expr {($tmax-$tmin)/double($nosteps)}]

    if { ![winfo exists .c] } {
        pack [canvas .c]
    } else {
        .c delete all
    }

    set tvalues {}
    set xvalues {}
    set xmin    {}
    set xmax    {}

    for {set i 0} {$i < $nosteps+1} {incr i} {
        set t [expr {$tmin + $i*$tstep}]
        if {[info command $expression]!={}} {
    # Note that RungeKutta expects start of time step as time
            set x [$expression [expr {$t-$tstep}] $x $tstep]
        } else {
            set x [expr $expression]
        }
        lappend tvalues $t
        lappend xvalues $x
        if {$xmin=={}} {set nvalues [llength $x]}
        foreach xvalu $x { ;# handle vector of results
            if { $xmin == {} || $xvalu < $xmin } {
                set xmin $xvalu
            }
            if { $xmax == {} || $xvalu > $xmax } {
                set xmax $xvalu
            }
        }
    }

    set tscale [::Plotchart::determineScale $tmin $tmax]
    set xscale [::Plotchart::determineScale $xmin $xmax]

    set p [::Plotchart::createXYPlot .c $tscale $xscale]
    set i 0
    foreach {name value} $args {
        switch -- $name {
            "-title" { .c create text 200 50 -text $value }
            "-colours" -
            "-colors" {
                set i 0
                foreach col $value {
                    $p dataconfig function$i -colour $col
                    incr i
                }
            }
            default {}
        }
    }
    foreach t $tvalues x $xvalues {
        set i 0
        while {$i<$nvalues} { ;# handle vector of results
            $p plot function$i $t [lindex $x $i]
            incr i
        }
    }
  }

  # now define a runge-kutta- solver
  package require math::calculus  ;# supplies RK4 algorithm
  
  proc dampened_oscillator {t xvec} { 
    # called by RungeKuttaStep. Returns list of derivatives of xvec.
        global damp
        global k
           set x  [lindex $xvec 0] ;# pos
           set speed [lindex $xvec 1] ;# speed
        # Solves the damped oscillator equation: x'' = -k*x-damp*speed
        # where x is position, and x'' the second time derivative.
        #  This equation is broken into 2 first order equations:
        # dx/dt=v
        # dv/dt= k.x-damp*speed
        #   k is the 'spring constant' of the oscillator;
        # more negative gives higher oscillation frequency.
        # ::math::calculus::rungeKuttaStep expects to receive
        # a list of derivatives at time t.
           return [list $speed [expr {$k*$x-$damp*$speed}]]
  }

  proc osc {t x tstep} { ;# return the value at the time t+dt for step tstep
    # x is a list of values (position, speed for the moving part)
     return [::math::calculus::rungeKuttaStep $t $tstep $x dampened_oscillator]
  }
    set damp -0.2
    set k -9
    set ready 1
    while {$damp<1.0} {
        plotf {0 10} osc -initial {0 1} -nosteps 100 -title "Damping $damp" -colours {red blue}
            set damp [expr {$damp+0.1}];
        after 1000 {incr ready} ;# delays next solution
        vwait ready
  }