1# -*-mode: tcl; fill-column: 75; tab-width: 8; coding: iso-latin-1-unix -*-
2#
3#       $Id: NPlot3d.tcl,v 1.9 2009-11-16 22:41:47 villate Exp $
4#
5###### NPlot3d.tcl ######
6############################################################
7# Netmath       Copyright (C) 1998 William F. Schelter     #
8# For distribution under GNU public License.  See COPYING. #
9############################################################
10
11# source plotting.tcl ; source nplot3d.tcl ; catch { destroy .plot3d} ;  plot3d -zfun "" -data $sample -xradius 10 -yradius 10
12# newidea:
13# { plot3d
14#  { gridequal {minx maxx} {miny maxy}
15#   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
16#  { grid {x0 x1  xm} {y0 y1 yn } miny maxy}
17#   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
18#  { xyzgrid {{x00 y00 z00 x01 y01 z01 .. x0  }{x0 x1  xm} {y0 y1 yn } miny maxy}
19#   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
20# tclMesh(2*[0,0,0,0,0;1,1,1,1,1]-1,2*[0,1,1,0,0;0,1,1,0,0]-1,2*[0,0,1,1,0;0,0,1,1,0]-1)
21
22#     { gridequal {
23
24# z00 z01 .. all belong to x=minx and y = miny,.... up y=maxy in n+1 steps
25#{ grid {minx maxx} {miny maxy}
26#  {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
27# }
28# where a mesh(1) {z00 z01 z11 z10} above
29
30
31
32# { mesh {{{x00 y00 z00 } { x01 y01 z01} { x02 y02 z02}  ..}{{x10 y10 z10} {x11 y11 z11} ......} ..}}
33# mesh(1) = P00 P01 P11 P10
34
35set sample { variable_grid { 0 1 2 } { 3 4 5} { {21 111 2} {3 4 5 } {6 7 8 }}}
36set sample { variable_grid { 0 1 2 } { 3 4 5} { {0 1 2} {3 4 5 } {6 7 8 }}}
37set sample { matrix_mesh {{0 1} { 2 3 } {4 5 }}  {{0 1} { 2 3 } {4 5 }}  {{0 1} { 2 3 } {4 5 }} }
38set sample { matrix_mesh {{0 1 2} {0 1 2 } {0 1 2 }} {{3 4 5} {3 4 5} {3 4 5}} { {0 1 2} {3 4 5 } {6 7 8 }}}
39set sample1 { variable_grid  { 1 2 3 4 5 6 7 8 9 10 }
40    { 1 2 3 }
41    {  { 0 0 0 0 0 0 0 0 0 0 }
42	{ 0 0.68404 1.28558 1.73205 1.96962 1.96962 1.73205 1.28558 0.68404 2.44921e-16 }
43	{ 0 1.36808 2.57115 3.4641 3.93923 3.93923 3.4641 2.57115 1.36808 4.89843e-16 }
44    }
45}
46
47set sample { matrix_mesh  {  { 0 0 0 0 0 }
48    { 1 1 1 1 1 }
49}  {  { 0 1 1 0 0 }
50    { 0 1 1 0 0 }
51}  {  { 0 0 1 1 0 }
52    { 0 0 1 1 0 }
53}
54}
55
56proc  fixupZ { } {
57    uplevel 1 {
58	if { [catch { expr $z + 0 } ] } {
59	    set z nam
60	}
61    }
62}
63
64proc vectorLength { v } {
65    expr { sqrt(1.0 * [lindex $v 0]*[lindex $v 0] + [lindex $v 1]*[lindex $v 1] + [lindex $v 2]*[lindex $v 2]) }
66}
67
68proc normalizeToLengthOne { v } {
69    set norm [expr { sqrt(1.0 * [lindex $v 0]*[lindex $v 0] + [lindex $v 1]*[lindex $v 1] + [lindex $v 2]*[lindex $v 2]) }]
70    if { $norm != 0.0 } {
71	return [list [expr { [lindex $v 0] / $norm  } ] \
72		    [expr { [lindex $v 1] / $norm  } ] \
73		    [expr { [lindex $v 2] / $norm  } ] ]
74
75    } else {
76	return "1.0 0.0 0.0 "
77    }
78}
79
80
81
82proc vectorCross { x1 x2 }  {
83    list \
84	[expr { [lindex $x1 1]*[lindex $x2 2]- [lindex $x2 1]*[lindex $x1 2]}] \
85	[expr { [lindex $x1 2]*[lindex $x2 0]- [lindex $x2 2]*[lindex $x1 0] } ] \
86	[expr { [lindex $x1 0]*[lindex $x2 1]- [lindex $x2 0]*[lindex $x1 1] }]
87}
88
89proc linspace { a b n } {
90    if { $n < 2 } { error [M [mc "from %s to %s requires at least 2 points"] "$a" "$b" ] }
91    set del [expr {($b - $a)*1.0/($n -1)  }]
92    for { set i 0 } { $i < $n } { incr i } {
93	lappend ans [expr {$a + $del * $i}]
94    }
95    return $ans
96}
97
98
99proc addOnePlot3d { win data } {
100    upvar #0 plot3dMeshes$win meshes
101    #puts " adding meshes = plot3dMeshes$win"
102    #puts "data=$data"
103    linkLocal $win points zmax zmin zcenter zradius rotationcenter xradius yradius xmin xmax ymin ymax lmesh
104    makeLocal $win flatten
105    oset $win colorfun plot3dcolorFun
106    oset $win cmap c1
107    global plot3dOptions
108    catch { unset  meshes }
109    set points ""
110
111    # check whether zradius is a number or "auto"
112    set dotruncate [expr ![catch {expr {$zradius + 1} }]]
113    if { $dotruncate } {
114		set zzmax [expr {$zcenter + $zradius}]
115		set zzmin [expr {$zcenter - $zradius}]
116    }
117
118    set k [llength $points]
119    set type [lindex $data 0]
120    # in general the data should be a list of plots..
121    if { [lsearch {grid mesh variable_grid matrix_mesh }  $type ]>=0 } {
122	set alldata [list $data]
123    } else {set alldata $data}
124    set lmesh {}
125    set first_mesh 0
126    foreach data $alldata {
127	set type [lindex $data 0]
128	if { "$type" == "grid" } {
129	    desetq "xmin xmax" [lindex $data 1]
130	    desetq "ymin ymax" [lindex $data 2]
131	    set pts [lindex $data 3]
132
133	    set ncols [llength $pts]
134	    set nrows  [llength [lindex $pts 0]]
135	    set data [list variable_grid [linspace $xmin $xmax $ncols] \
136			  [linspace $ymin $ymax $nrows] \
137			  $pts ]
138	    set type "variable_grid"
139	}
140	if { "$type" == "variable_grid" } {
141	    set first_mesh [llength $lmesh]
142	    desetq "xrow yrow zmat" [lrange $data 1 end]
143	    # puts "xrow=$xrow,yrow=$yrow,zmat=$zmat"
144	    set nx [expr {[llength $xrow] -1}]
145	    set ny [expr {[llength $yrow] -1}]
146	    #puts "nx=$nx,ny=$ny"
147	    #	set xmin [lindex $xrow 0]
148	    #	set xmax [lindex $xrow $nx]
149	    #	set ymin [lindex $yrow 0]
150	    #	set ymax [lindex $yrow $ny]
151	    desetq "xmin xmax" [minMax $xrow ""]
152	    desetq "ymin ymax" [minMax $yrow ""]
153	    desetq "zmin zmax" [matrixMinMax [list $zmat]]
154	    if { $dotruncate } {
155		set zmax  $zzmax
156		set zmin  $zzmin
157	    }
158	    #	puts "and now"
159	    #	dshow nx xmin xmax ymin ymax zmin zmax
160
161	    for {set j 0} { $j <= $ny } { incr j} {
162		set y [lindex $yrow $j]
163		set row [lindex $zmat $j]
164		for {set i 0} { $i <= $nx } { incr i} {
165		    set x [lindex $xrow $i]
166		    set z [lindex $row $i]
167		    #puts "x=$x,y=$y,z=$z, at ($i,$j)"
168		    fixupZ
169		    if { $j != $ny && $i != $nx } {
170			lappend lmesh [list $k [expr { $k+3 }] \
171					   [expr { $k+3+($nx+1)*3 }] \
172					   [expr { $k+($nx+1)*3 }]]
173		    }
174		    incr k 3
175		    lappend points $x $y $z
176		}
177	    }
178	    setupPlot3dColors $win $first_mesh
179	} elseif { "$type" == "matrix_mesh" } {
180	    set first_mesh [llength $lmesh]
181	    desetq "xmat ymat zmat" [lrange $data 1 end]
182	    foreach v {x y z} {
183		desetq "${v}min ${v}max" [matrixMinMax [list [set ${v}mat]]]
184	    }
185	    if { $dotruncate } {
186		set zmax  $zzmax
187		set zmin  $zzmin
188	    }
189	    set nj [expr {[llength [lindex $xmat 0]] -1 }]
190	    set ni [expr {[llength $xmat ] -1 }]
191	    set i -1
192	    set k [llength $points]
193	    foreach rowx $xmat rowy $ymat rowz $zmat {
194		set j -1
195		incr i
196		if { [llength $rowx] != [llength $rowy] } {
197		    error [concat [mc "mismatch"] "rowx:$rowx,rowy:$rowy"]
198		}
199		if { [llength $rowx] != [llength $rowz] } {
200		    error [concat [mc "mismatch"] "rowx:$rowx,rowz:$rowz"]
201		}
202		foreach x $rowx y $rowy z $rowz {
203		    incr j
204		    if { $j != $nj && $i != $ni } {
205			#puts "tes=($i,$j) $x, $y, $z"
206			lappend lmesh [ list \
207					    $k [expr { $k+3 } ] [expr { $k + 3  + ($nj+1)*3}] \
208					    [expr { $k+($nj+1)*3 }] ]
209		    }
210		    incr k 3
211		    lappend points $x $y $z
212		}
213	    }
214	    setupPlot3dColors $win $first_mesh
215	} elseif { 0 && "$type" == "mesh" } {
216	    set first_mesh [llength $lmesh]
217	    # walk thru compute the xmin, xmax, ymin , ymax...
218	    # and then go thru setting up the mesh array..
219	    # and maybe setting up the color map for these meshes..
220	    #
221	    # { mesh {{{x00 y00 z00 } { x01 y01 z01} { x02 y02 z02}  ..}{{x10 y10 z10} {x11 y11 z11} ......} ..}}
222	    # mesh(1) = P00 P01 P11 P10
223	    set mdata [lindex $data 1]
224	    set nx [llength $mdata]
225	    set ny [llength [lindex $mdata 0]]
226
227	    for {set i 0} { $i <= $nx } { incr i} {
228		set pts [lindex $mdata $i]
229		set j 0
230		foreach { x y z} $pts {
231		    fixupZ $z
232		    if { $j != $ny && $i != $nx } {
233			lappend lmesh [list
234				       $k [expr { $k+3 }] [expr { $k+3+($ny+1)*3 }] \
235					   [expr { $k+($ny+1)*3 }] ]
236		    }
237		}
238		incr k 3
239		lappend points $x $y $z
240		incr j
241	    }
242	    setupPlot3dColors $win $first_mesh
243	} elseif { "[assq $type $plot3dOptions notthere]" != "notthere" } {
244	    oset $win $type [lindex $data 1]
245	    if { $type == "zradius" } {
246		# check whether zradius is a number or "auto"
247		set dotruncate [expr ![catch {expr {$zradius + 1} }]]
248		if { $dotruncate } {
249		    set zzmax [expr {$zcenter + $zradius}]
250		    set zzmin [expr {$zcenter - $zradius}]
251		}
252	    }
253	}
254	if { $first_mesh != [llength $lmesh] } {
255	    # set up the min/max values for the complete plot
256	    foreach v { x y z } {
257		if { [info exists ${v}gmin] } {
258		    if { [set ${v}min] < [set ${v}gmin] } {
259			set ${v}gmin [set ${v}min]
260		    }
261		} else {
262		    set ${v}gmin [set ${v}min]
263		}
264		if { [info exists ${v}gmax] } {
265		    if { [set ${v}max] > [set ${v}gmax] } {
266			set ${v}gmax [set ${v}max]
267		    }
268		} else {
269		    set ${v}gmax [set ${v}max]
270		}
271	    }
272	}
273    }
274    foreach v { x y z } {
275	set ${v}min [set ${v}gmin]
276	set ${v}max [set ${v}gmax]
277	set a [set ${v}min]
278	set b  [set ${v}max]
279	if { $a == $b } {
280	    set ${v}min [expr {$a -1}]
281	    set ${v}max [expr {$a +1}]
282	}
283	set ${v}radius [expr {($b - $a)/2.0}]
284	set ${v}center [expr {($b + $a)/2.0}]
285    }
286    set rotationcenter "[expr {.5*($xmax + $xmin)}] [expr {.5*($ymax + $ymin)}]   [expr {.5*($zmax + $zmin)}] "
287
288    #puts "meshes data=[array get meshes]"
289    #global plot3dMeshes.plot3d
290    #puts "array names plot3dMeshes.plot3d = [array names plot3dMeshes.plot3d]"
291}
292
293proc vectorDiff { x1 x2 } {
294    list [expr { [lindex $x1 0] - [lindex $x2 0] }] \
295	[expr { [lindex $x1 1] - [lindex $x2 1] }] \
296	[expr { [lindex $x1 2] - [lindex $x2 2] }]
297}
298
299
300proc oneCircle { old2 old1 pt radius nsides { angle 0 } } {
301    set dt  [expr {  3.14159265358979323*2.0/($nsides-1.0) + $angle }]
302    for  { set i 0 } { $i < $nsides } { incr i } {
303	set t [expr {$dt*$i }]
304	lappend ans [expr { $radius*([lindex $old2 0]*cos($t) + [lindex $old1 0] * sin($t)) + [lindex $pt 0] } ] \
305	    [expr { $radius*([lindex $old2 1]*cos($t) + [lindex $old1 1] * sin($t)) + [lindex $pt 1] } ] \
306	    [expr { $radius*([lindex $old2 2]*cos($t) + [lindex $old1 2] * sin($t)) + [lindex $pt 2] } ]
307    }
308    return $ans
309}
310
311proc curve3d { xfun yfun zfun trange } {
312    foreach u { x y z} {
313	set res [parseConvert [set ${u}fun] -variables t]
314	proc _${u}fun { t } [list expr [lindex [lindex $res 0] 0]]
315    }
316}
317
318proc tubeFromCurveData { pts nsides radius } {
319    set n [llength $pts] ;
320    set closed [ expr { [vectorLength [vectorDiff [lindex $pts 0] [lindex $pts end]]] < .02} ]
321    if { $closed } {
322	set f1 [expr {$n -2}]
323	set f2 1
324    } else {
325	set f1 0
326	set f2 1
327    }
328    set delta [vectorDiff [lindex $pts $f2] [lindex $pts $f1]]
329    if { [lindex $delta 0] == 0 && [lindex $delta 1] == 0 && [lindex $delta 2] == 0 } { set delta "0 0 1.0" }
330    set old ".6543654 0.0765456443 0.2965433"
331    set old1 [normalizeToLengthOne [vectorCross $delta $old]]
332    set n1 $old1
333    set n2 [normalizeToLengthOne [vectorCross $delta $old1]]
334    set first1 $n1 ; set first2 $n2
335
336    lappend ans [oneCircle $n2   old1 [lindex $pts 0]]
337    for { set j 1 } { $j < $n -1 } { incr j } {
338	set delta [vectorDiff [lindex $pts $j] [lindex $pts [expr {$j+1}]]]
339	if { [lindex $delta 0] == 0 && [lindex $delta 1] == 0 && [lindex $delta 2] == 0 } { set delta $old
340	}
341	set old $delta
342	set old1 [normalizeToLengthOne [vectorCross $delta $n1]]
343	set old2 [normalizeToLengthOne [vectorCross $delta $n2]]
344	set n2 $old1
345	set n1 $old2
346	lappend ans [oneCircle $n2 $n1 [lindex $pts $j] $radius $nsides]
347    }
348    if { $closed } {
349	set f2 1 ; set f1 [expr {$n -2}] ; set f3 0
350    } else {
351	set f1 [expr {$n -2}] ; set f2 [expr {$n-1}] ; set f3 $f2
352    }
353
354    set delta [vectorDiff [lindex $pts $f2] [lindex $pts $f1]]
355    if { [lindex $delta 0] == 0 && [lindex $delta 1] == 0 && \
356	     [lindex $delta 2] == 0 } { set delta $old }
357    set old1 [normalizeToLengthOne [vectorCross delta $n1]]
358    set old2 [normalizeToLengthOne [vectorCross $n2 $delta]]
359    set n2 $old1 ; set n1 $old2
360    if { $closed } {
361	set angle [vangle $first1 $n1]
362	set n1 $first1 ; st n2 $first2;
363    }
364    lappend ans [oneCircle $n2 $n1 [lindex $pts $f3] $radius $nsides $angle]
365    return $ans
366}
367
368
369#
370#-----------------------------------------------------------------
371#
372# vangle --  angle between two unit vectors
373#
374#  Results: an angle
375#
376#  Side Effects: none.
377#
378#----------------------------------------------------------------
379#
380proc vangle { x1 x2 } {
381    set dot [expr { [lindex $x1 0]*[lindex $x2 0] +\
382			[lindex $x1 1]*[lindex $x2 1] +\
383			[lindex $x1 2]*[lindex $x2 2]} ]
384    if { $dot >= 1 } { return 0.0 }
385    if { $dot <= -1.0 } { return 3.141592653589 }
386    return [expr { acos($dot) } ]
387}
388
389## endsource nplot3d.tcl
390