# solids3d.tcl --
#
#    Package for displaying 3D solid bodies
#    (sample Workbench module)
#
# Notes:
#    This package is a quick hack to get started only
#
# Version information:
#    version 0.3: added mapping algorithm for general transformations,
#                 september 2002
#    version 0.2: improved sorting of polygons, september 2002
#                 (with help from Donal Fellows)
#    version 0.1: initial implementation, september 2002

package provide Solids3D 0.1

namespace eval ::Solids3D {

    namespace export func deriv

    variable display_options

# normalToPlane --
#    Return the normal vector to a plane given as three points
#
# Arguments:
#    point1       First point
#    point2       Second point
#    point3       Third point
#
# Result:
#    Vector {x,y,z}, a normal vector to the plane
#
proc normalToPlane {point1 point2 point3} {
    foreach {dummy px1 py1 pz1} $point1 {break}
    foreach {dummy px2 py2 pz2} $point2 {break}
    foreach {dummy px3 py3 pz3} $point3 {break}
    set vx1 [expr {$px2-$px1}]
    set vy1 [expr {$py2-$py1}]
    set vz1 [expr {$pz2-$pz1}]
    set vx2 [expr {$px3-$px1}]
    set vy2 [expr {$py3-$py1}]
    set vz2 [expr {$pz3-$pz1}]

    set nx  [expr {$vy1*$vz2-$vz1*$vy2}]
    set ny  [expr {$vz1*$vx2-$vx1*$vz2}]
    set nz  [expr {$vx1*$vy2-$vy1*$vx2}]

    return [list VECTOR $nx $ny $nz]
}

# pointOnTorus --
#    Return the coordinates of a point on a torus, as given by
#    two parameters (angles)
#
# Arguments:
#    phi          Angle with respect to x-axis
#    theta        Angle with respect to "inner" axis of torus
#
# Result:
#    Point {x,y,z}, the coordinates of the point
#
proc pointOnTorus {phi theta} {
    set cosphi [expr {cos($phi)}]
    set sinphi [expr {sin($phi)}]
    set costh4 [expr {0.25*cos($theta)}]
    set sinth4 [expr {0.25*sin($theta)}]
    set x [expr {$cosphi*(1.0+$costh4)}]
    set y [expr {$sinphi*(1.0+$costh4)}]
    set z [expr {1.0+$sinth4}]

    return [list POINT $x $y $z]
}

# constructTorus --
#    Return a list of polygons, together forming an approximation to a
#    torus
#
# Arguments:
#    nophi        Number of steps along main perimeter
#    notheta      Number of steps along secondary perimeter
#
# Result:
#    List of polygons
#
proc constructTorus {nophi notheta} {
    variable angle

    set pi     3.1415926
    set dphi   [expr {2.0*$pi/double($nophi)}]
    set dtheta [expr {2.0*$pi/double($notheta)}]

    set polygons {POLYGON-LIST}

    for { set iphi 0 } { $iphi < $nophi } { incr iphi } {
      set phi1 [expr {$dphi*double($iphi)}]
      set phi2 [expr {$dphi*double($iphi+1)}]

      for { set itheta 0 } { $itheta < $notheta } { incr itheta } {
         set theta1 [expr {$dtheta*double($itheta)}]
         set theta2 [expr {$dtheta*double($itheta+1)}]

         set point1 [pointOnTorus $phi1 $theta1]
         set point2 [pointOnTorus $phi1 $theta2]
         set point3 [pointOnTorus $phi2 $theta2]
         set point4 [pointOnTorus $phi2 $theta1]

         set red    [expr {int(125*(1.0+cos($phi1)))}]
         set green  [expr {int(125*(1.0+sin($phi1)))}]
         set blue   [expr {int(125*(1.0+sin($theta1)))}]
         set colour [format "#%2.2X%2.2X%2.2X" $red $green $blue]
         set normal [normalToPlane $point1 $point2 $point4]

         lappend polygons \
            [list POLYGON $colour $normal $point1 $point2 $point3 $point4]
       }
    }
    return $polygons
}

# constructCube --
#    Return a list of polygons, together forming a cube
#
# Arguments:
#    None
#
# Result:
#    List of polygons
#
proc constructCube {} {

    set v  {VECTOR  1  1  1}
    set p1 {POINT  -1 -1 -1}
    set p2 {POINT   1 -1 -1}
    set p3 {POINT   1  1 -1}
    set p4 {POINT  -1  1 -1}
    set p5 {POINT  -1 -1  1}
    set p6 {POINT   1 -1  1}
    set p7 {POINT   1  1  1}
    set p8 {POINT  -1  1  1}
    set f1 [list POLYGON blue   $v $p1 $p2 $p3 $p4]
    set f2 [list POLYGON green  $v $p1 $p2 $p6 $p5]
    set f3 [list POLYGON yellow $v $p2 $p3 $p7 $p6]
    set f4 [list POLYGON purple $v $p3 $p4 $p8 $p7]
    set f5 [list POLYGON orange $v $p1 $p4 $p8 $p5]
    set f6 [list POLYGON red    $v $p5 $p6 $p7 $p8]

    list POLYGON-LIST $f1 $f2 $f3 $f4 $f5 $f6
}

# constructSpiral --
#    Return a list of polygons, together forming a spiral
#
# Arguments:
#    nosegments    Number of segments
#
# Result:
#    List of polygons
#
proc constructSpiral {nosegments} {

    set result {POLYGON-LIST}

    set dz   [expr {2.0/double($nosegments)}]
    set dphi [expr {4.0*3.1415926/double($nosegments)}]
    set r1   0.45
    set r2   0.95
    set v  {VECTOR  1  1  1}
    for {set i 0} {$i < $nosegments} {incr i} {
       set z1   [expr {-1.0+$dz*$i}]
       set z2   [expr {$z1+$dz}]
       set phi1 [expr {$dphi*$i}]
       set phi2 [expr {$phi1+$dphi}]
       set x1   [expr {$r1*cos($phi1)}]
       set x2   [expr {$r2*cos($phi1)}]
       set y1   [expr {$r1*sin($phi1)}]
       set y2   [expr {$r2*sin($phi1)}]
       set x3   [expr {$r1*cos($phi2)}]
       set x4   [expr {$r2*cos($phi2)}]
       set y3   [expr {$r1*sin($phi2)}]
       set y4   [expr {$r2*sin($phi2)}]

       set p1   [list POINT $x1 $y1 $z1]
       set p2   [list POINT $x2 $y2 $z1]
       set p3   [list POINT $x3 $y3 $z2]
       set p4   [list POINT $x4 $y4 $z2]

       set red    [expr {int(125*(1.0+cos($phi1)))}]
       set green  [expr {int(125*(1.0+sin($phi1)))}]
       set blue   255
       set colour [format "#%2.2X%2.2X%2.2X" $red $green $blue]

       lappend result [list POLYGON $colour $v $p1 $p2 $p4 $p3]
    }
    return $result
}

# rotateX,Y,Z --
#    Rotate around the x,y or z axis
#
# Arguments:
#    angle        Angle over which to rotate
#    point        Point to rotate
#
# Result:
#    New coordinates
#
proc rotateX {angle point} {
    foreach {dummy x y z} $point break
    set xr $x
    set yr [expr {$y*cos($angle)-$z*sin($angle)}]
    set zr [expr {$y*sin($angle)+$z*cos($angle)}]

    return [list POINT $xr $yr $zr]
}
proc rotateY {angle point} {
    foreach {dummy x y z} $point break
    set xr [expr {$x*cos($angle)-$z*sin($angle)}]
    set yr $y
    set zr [expr {$x*sin($angle)+$z*cos($angle)}]

    return [list POINT $xr $yr $zr]
}
proc rotateZ {angle point} {
    foreach {dummy x y z} $point break
    set xr [expr {$x*cos($angle)-$y*sin($angle)}]
    set yr [expr {$x*sin($angle)+$y*cos($angle)}]
    set zr $z

    return [list POINT $xr $yr $zr]
}

# projectOnYZ --
#    Project a point on the YZ-plane (and scale as we do this)
#
# Arguments:
#    point        Point to rotate
#
# Result:
#    XY coordinates (for the screen)
#
proc projectOnYZ {point} {
    variable scale
    variable yoffset
    variable zoffset
    foreach {dummy x y z} $point break
    set xp [expr {int($scale*($y-$yoffset))}]
    set yp [expr {int($scale*($zoffset-$z))}]

    return [list $xp $yp]
}

# compareView --
#    Compare the relative position of two polygons w.r.t. the viewpoint
#
# Arguments:
#    polygon1     First polygon
#    polygon2     Second polygon
#
# Result:
#    1, -1 if the first polygon lies before or after the second,
#
# Note:
#    The comparison is made by looking at the x-coordinate. As the
#    projection is in a plane perpendicular to the x-axis, this is
#    all we need to do.
#
proc compareView {polygon1 polygon2} {
    variable angle
    variable viewpoint

    #
    # Determine the smallest z-coordinate for each point
    #
    set xmin1 {}
    foreach point [lrange $polygon1 3 end] {
       foreach {dummy px py pz} $point {break}
       if { $xmin1 == {} || $xmin1 > $px } {
          set xmin1 $px
       }
    }

    set xmin2 {}
    foreach point [lrange $polygon2 3 end] {
       foreach {dummy px py pz} $point {break}
       if { $xmin2 == {} || $xmin2 > $px } {
          set xmin2 $px
       }
    }

    if { $xmin1 < $xmin2 } {
       return 1
    } else {
       return -1
    }
}

# applyTorsion --
#    Transformation modelling torsion around Y-axis
#
# Arguments:
#    point        Point to "tord"
#
# Result:
#    New point
#
proc applyTorsion {point} {
    variable torsion_angle

    foreach {type x y z} $point break
    set angle [expr {$torsion_angle*$y}]
    set xr [expr {$x*cos($angle)-$z*sin($angle)}]
    set yr $y
    set zr [expr {$x*sin($angle)+$z*cos($angle)}]

    return [list $type $xr $yr $zr]
}

# rotateCube --
#    Rotate a cube over three axes
#
# Arguments:
#    point        Point to rotate
#
# Result:
#    New point
#
proc rotateCube {point} {
    variable angle

    set anglex [expr {1.23*$angle}]
    set angley [expr {$angle+0.2}]
    set anglez [expr {-0.81*$angle+0.44}]

    rotateX $anglex [rotateY $angley [rotateZ $anglez $point]]
}

# scaleZ --
#    Scale a point in the Z-axis
#
# Arguments:
#    point        Point to "scale"
#
# Result:
#    New point
#
proc scaleZ {point} {
    variable scale_factor

    foreach {type x y z} $point break
    set zr [expr {$z*$scale_factor}]

    return [list $type $x $y $zr]
}

# map --
#    Map a list of polygons via a transformation of single points and
#    vectors
#
# Arguments:
#    mapPoint     Transformation of a single point or vector
#    polygons     List of 3D polygons
#
# Result:
#    List of transformed polygons
#
proc map {mapPoint polygons} {
    #
    # Check the type of the second argument (it should be a list whose
    # first element gives information about the type):
    # - For polygon-lists, call [map] on the individual polygons
    #   and assemble
    # - For polygons, call [mapPoint] on the individual points
    #   and assemble
    #
    switch -- [lindex $polygons 0] {
    "POLYGON-LIST" {
       set result [lindex $polygons 0]
       foreach polygon [lrange $polygons 1 end] {
          lappend result [map $mapPoint $polygon]
       }
    }
    "POLYGON" {
       # $polygons is actually a single polygon now:
       # element 0 - the keyword
       # element 1 - the colour
       # element 2 - the normal vector (should be re-evaluated,
       #             but right now it is not used)
       set result [lrange $polygons 0 2]
       foreach point [lrange $polygons 3 end] {
          lappend result [$mapPoint $point]
       }
    }
    default {
       set result $polygons
    }
    }

    return $result
}

# display --
#    Quick and dirty implementation to display a set of polygons
#
# Arguments:
#    polygons     List of 3D polygons
#
# Result:
#    None
#
# Side effect:
#    Display of polygons in the canvas
#
proc display {polygons} {
    variable angle
    variable viewpoint

    set viewpoint [list POINT [expr cos($angle)] 0.0 [expr sin($angle)]]

    #
    # Rotate the polygons first, so that the new coordinates can
    # be used directly for sorting and plotting
    #
    set plane_polygons   [lrange $polygons 1 end]
    set rotated_polygons {}

    foreach polygon $plane_polygons {
       set points [lrange $polygon 3 end]
       set coords [lrange $polygon 0 2]

       foreach point $points {
          lappend coords [rotateY $angle $point]
       }
       lappend rotated_polygons $coords
    }

    #
    # Now sort the polygons w.r.t. the z-coordinate
    #
    set sorted_polygons [lsort -command compareView $rotated_polygons]

    .cnv delete all
    foreach polygon $sorted_polygons {
       set colour [lindex $polygon 1]
       set points [lrange $polygon 3 end]
       set coords {}
       foreach point $points {
          set coords [concat $coords [projectOnYZ $point]]
       }
      .cnv create polygon $coords -fill $colour -outline black
    }
}

#
# Initialise the variables
#
variable angle   [expr {0.25*3.1415926}]
variable scale    70.0
variable yoffset  -2.0
variable zoffset   2.0
variable torsion_angle 0.0

} ;# End of namespace

#
# Set up a simple UI
#
frame       .f
radiobutton .f.radio1 -value torus  -variable object -text "Contorted torus" -anchor w
#radiobutton .f.radio2 -value orange -variable object -text "Orange peals"    -anchor w
radiobutton .f.radio3 -value cube   -variable object -text "Colourful cube"  -anchor w
radiobutton .f.radio4 -value spring -variable object -text "Spring"          -anchor w
label       .f.label  -text " "                   -anchor w
button      .f.show   -text "Show"  -command showAnimation
pack .f.radio1 -side top -fill x
#pack .f.radio2 -side top -fill x
pack .f.radio3 -side top -fill x
pack .f.radio4 -side top -fill x
pack .f.label  -side top -fill x
pack .f.show   -side top -fill x
pack .f -side left


canvas .cnv -width 300 -height 300 -background white
pack .cnv -fill both
button      .exit   -text "Exit"  -command exit
pack .exit -side bottom

set ::object "torus"

#
# The animation commands showAnimation and next
#

set ::torus  [::Solids3D::constructTorus 16 16]
set ::cube   [::Solids3D::constructCube]
set ::spring [::Solids3D::constructSpiral 40]

proc showAnimation {} {
    global object
    global step

    set step 0
    next $step
}

proc next {step} {
    global object
    set maxstep 40

    switch -- $object {
    "torus" {
       set ::Solids3D::torsion_angle [expr {1.57*$step/20.0}]
       ::Solids3D::display [::Solids3D::map applyTorsion $::torus]
    }
    "cube" {
       set maxstep 100
       set ::Solids3D::angle [expr {0.1*$step}]
       ::Solids3D::display [::Solids3D::map rotateCube $::cube]
    }
    "orange" {
       set ::Solids3D::angle [expr {0.1*$step}]
       ::Solids3D::display [::Solids3D::map rotateCube $::cube]
    }
    "spring" {
       set ::Solids3D::angle [expr {0.03*$step}]
       set ::Solids3D::scale_factor [expr {($step+5)/30.0}]
       ::Solids3D::display \
          [::Solids3D::map rotateCube \
             [::Solids3D::map scaleZ $::spring]]
    }
    }
    if { $step < $maxstep } {
       incr step
       after 100 next $step
    }
}