# 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.2: improved sorting of polygons, september 2002 
# 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 zdepth [lindex [rotateY $angle $point1] 3] 
		set normal [normalToPlane $point1 $point2 $point4] 

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

    # rotateY -- 
    # Rotate around the y axis 
    # 
    # Arguments: 
    # angle Angle over which to rotate 
    # point Point to rotate 
    # 
    # Result: 
    # New coordinates 
    # 
    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] 
    } 

    # 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, 
    # 0, if they intersect 
    # 
    # Note: 
    # The third element of the list that constitutes a polygon (index 2!) 
    # is supposed to be a vector normal to the polygon. For reasons of 
    # speed. 
    # 
    proc compareView2 {polygon1 polygon2} { return 1 } 

    proc compareView {polygon1 polygon2} { 
	variable angle 
	variable viewpoint 

	# 
	# Determine the relevant sign 
	# 
	foreach {dummy nx ny nz} [lindex $polygon1 2] {break} 
	foreach {dummy vx vy vz} $viewpoint {break} 
	set vn [expr {$nx*$vx+$ny*$vy+$nz*$vz}] 
	if { $vn < 0 } { 
	    set nx [expr -$nx] 
	    set ny [expr -$ny] 
	    set nz [expr -$nz] 
	} 

	# 
	# For each point of the second polygon, does the inproduct have 
	# the same sign? If the number of identical signs is 4, then 
	# polygon lies either on the same side as the viewpoint or 
	# on the opposite side, but we have reached a conclusion. 
	# 
	set count 0 
	foreach point [lrange $polygon2 3 end] { 
	    foreach {dummy px py pz} $point {break} 
	    set pn [expr {$nx*$px+$ny*$py+$nz*$pz}] 
	    if { $pn > 0.0 } { 
		incr count 
	    } else { 
		incr count -1 
	    } 
	} 

	if { $count == 4 } { 
	    return -1 
	} elseif { $count == -4 } { 
	    return 1 
	} else { 
	    set rc [compareView $polygon2 $polygon1] 
	    expr {-$rc} 
	} 
    } 

    # 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)]] 

	# 
	# Sort the polygons first 
	# Note: 
	# The comparison is too simple, but for now it should work 
	# 
	set plane_polygons [lrange $polygons 1 end] 

	set sorted_polygons [lsort -command compareView $plane_polygons] 

	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 [rotateY $angle $point]]] 
	    } 
	    .cnv create polygon $coords -fill $colour -outline black 
	} 
    } 

    # 
    # Initialise the variables 
    # 
    variable angle [expr {0.25*3.1415926}] 
    variable scale 100.0 
    variable yoffset -1.5 
    variable zoffset 2.0 

} ;# End of namespace 

# 
# Run the program 
# 

canvas .cnv -width 300 -height 300 -background white 
pack .cnv -fill both 
set torus [::Solids3D::constructTorus 16 16] 
::Solids3D::display $torus 