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