# fourier.tcl --
 #    Package to manipulate Fourier series
 #
 # Author: Arjen Markus ([email protected])
 #
 # Notes on the implementation:
 # 1. The Fourier series is represented as a tagged list, where
 #    the tag is either FOURIER or FOURIER-HALF. Aftre the tag
 #    follow two values, the amplitudes for the cosine and sine
 #    components.
 # 2. The interval of interest is always [-1,1]
 # 3. The function is supposed to be periodic on an interval [-1,1]
 #    or on the interval [-2,2], indicated as FOURIER-HALF
 #
 
 # math::fourier --
 #    Namespace for the commands
 #
 namespace eval ::math::fourier {
    namespace export sineFourier toTimeseries
 
    # Possible extension: plenty!
 }
 
 # sineFourier --
 #    Construct a Fourier series based on sines only from an expression
 #    for the amplitudes (so an antisymmetric function results)
 #
 # Arguments:
 #    ncomps      Number of components
 #    var         Name of the variable in the expression (the index of
 #                the component)
 #    expression  Expression to determine the amplitudes
 #    option      Whether the interval [-1,1] covers the whole
 #                period or only half of it (-half). Optional
 # Return value:
 #    Tagged list
 #
 proc ::math::fourier::sineFourier { ncomps var expression {option {}} } {
    if { $option == "-half" } {
       set series "FOURIER-HALF"
    } else {
       set series "FOURIER"
    }
 
    #
    # TODO: properly evaluate other variables in caller's scope
    #
    lappend series 0.0 0.0
 
    for { set $var 0 } { [set $var] < $ncomps } { incr $var } {
       set ampl [expr $expression]
       lappend series 0.0 $ampl
    }
 
    return $series
 }
 
 # toTimeseries --
 #    Return the values for equally spaced points
 #
 # Arguments:
 #    series      Fourier series to use
 #    npoints     Number of points
 # Return value:
 #    List of values
 #
 proc ::math::fourier::toTimeseries { series npoints } {
    set tag [lindex $series 0]
    if { $tag == "FOURIER-HALF" } {
       set period -0.5
    } elseif { $tag == "FOURIER" } {
       set period  0.0
    } else {
       error "Series is not a Fourier series"
    }
 
    if { $npoints < 2 } {
       error "Number of points must be at least 2"
    }
 
    set ncomps [expr {([llength $series]-1)/2}]
    set delx   [expr {2.0/($npoints-1)}]
 
    set a0     [lindex $series 1]
    for { set j 0 } { $j < $npoints } { incr j } {
       lappend result $a0
    }
 
    set elem 3
    for { set i 1 } { $i < $ncomps } { incr i } {
       set acos [lindex $series $elem]
       incr elem
       set asin [lindex $series $elem]
       incr elem
 
       set period [expr {$period+1.0}]
 
       set newresult {}
       for { set j 0 } { $j < $npoints } { incr j } {
          set x      [expr {-1.0+$j*$delx}]
          set resx   [lindex $result $j]
 
          set acosx  [expr {$acos * cos(3.1415926*$period*$x)}]
          set asinx  [expr {$asin * sin(3.1415926*$period*$x)}]
          set resx   [expr {$resx + $acosx + $asinx}]
          lappend newresult $resx
       }
       set result $newresult
    }
 
    return $result
 }
 
 # plotTimeseries --
 #    Plot a timeseries
 #
 # Arguments:
 #    series      Fourier series to use
 #    npoints     Number of points
 #    colour      Which colour to use
 # Return value:
 #    List of values
 #
 proc ::math::fourier::plotTimeseries { series npoints colour } {
    set yvalues [toTimeseries $series $npoints]
 
    set x0      -1.0
    set dx      [expr {2.0/($npoints-1)}]
    set xvalues {}
    set x1      $x0
    set y1      [lindex $yvalues 0]
    for { set i 1 } { $i < $npoints } { incr i } {
       set x2 [expr {$x0+$dx*$i}]
       set y2 [lindex $yvalues $i]
 
       foreach {sx1 sy1} [scaleToCanvas $x1 $y1] {break}
       foreach {sx2 sy2} [scaleToCanvas $x2 $y2] {break}

       .c create line $sx1 $sy1 $sx2 $sy2 -fill $colour

       set x1 $x2
       set y1 $y2
    }
 }

 # scaleToCanvas --
 #    Scale points for drawing on the canvas
 #
 # Arguments:
 #    x           X coordinate
 #    y           Y coordinate
 # Return value:
 #    List of scaled x and y
 #
 proc scaleToCanvas { x y } {
    global scale
 
    set sx [expr {$scale(pxmin)+($x-$scale(xmin))*$scale(xfactor)}]
    set sy [expr {$scale(pymax)-($scale(ymax)-$y)*$scale(yfactor)}]

    return [list $sx $sy]
 }

 # setCanvasScales --
 #    Set the drawing area and the extreme coordinates for drawing
 #
 # Arguments:
 #    area        List of pixel coordinates: xmin, ymin, xmax, ymax
 #    coords      List of world coordinates: xmin, ymin, xmax, ymax
 # Return value:
 #    None
 # Side effect:
 #    Calculation of coordinate transformation
 #
 proc setCanvasScale { area coords } {
    global scale
 
    set scale(pxmin) [lindex $area   0]
    set scale(pymax) [lindex $area   1] ;# Because of inversion of y-coordinate
    set scale(xmin)  [lindex $coords 0]
    set scale(ymax)  [lindex $coords 3]

    set xmin         [lindex $coords 0]
    set xmax         [lindex $coords 2]
    set pxmax        [lindex $area   2]
    set xfactor [expr {($pxmax-$scale(pxmin))/($xmax-$xmin)}]

    set ymin         [lindex $coords 1]
    set ymax         [lindex $coords 3]
    set pymin        [lindex $area   3]
    set yfactor [expr {($scale(pymax)-$pymin)/($ymax-$ymin)}]

    set scale(xfactor) $xfactor
    set scale(yfactor) $yfactor

    return
 } 
 
 #
 # Some simple tests
 #
 catch { console show }

  if {[file tail $::argv0] == [file tail [info script]]} {
    namespace import ::math::fourier::*

    # This Fourier series converges to the function:
    # f(x) = -1/2 if x <   0 
    # f(x) =   1/2 if x >= 0 

    set f [sineFourier 20 k {1.0/($k+0.5)/3.1415926} -half]
    puts $f
    set values [toTimeseries $f 41]
    puts $values

    if { [info exists tk_version] } {
       canvas .c -width 400 -height 300
       pack   .c -fill both

       setCanvasScale {50 50 380 280} {-1.0 -1.0  1.0 1.0}
       foreach {xaxis1 yaxis1} [scaleToCanvas -1.0 -1.0] {break}
       foreach {xaxis2 yaxis2} [scaleToCanvas  1.0  1.0] {break}
       foreach {xaxis3 yaxis3} [scaleToCanvas  0.0  0.0] {break}

       .c create line $xaxis1 $yaxis3 $xaxis2 $yaxis3 -fill black
       .c create line $xaxis3 $yaxis1 $xaxis3 $yaxis2 -fill black

       foreach {ncomps colour} {
          1 blue 2 green 3 yellow 4 red 5 magenta 20 black} {

          set series [sineFourier $ncomps k {1.0/($k+0.5)/3.1415926} -half]
          plotTimeseries $series 100 $colour
       }
    }
 }