interfaces/xmaxima/Tkmaxima/NPlot3d.tcl

   1 ############################################################
   2 # NPlot3d.tc                                              #
   3 # Copyright (C) 1998 William F. Schelter                   #
   4 # For distribution under GNU public License.  See COPYING. #
   5 #                                                          #
   6 #     Time-stamp: "2024-04-03 19:04:26 villate"            #
   7 ############################################################
   8
   9 # source plotting.tcl ; source nplot3d.tcl ; catch { destroy .plot3d} ;  plot3d -zfun "" -data $sample -xradius 10 -yradius 10
  10 # newidea:
  11 # { plot3d
  12 #  { gridequal {minx maxx} {miny maxy}
  13 #   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
  14 #  { grid {x0 x1  xm} {y0 y1 yn } miny maxy}
  15 #   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
  16 #  { xyzgrid {{x00 y00 z00 x01 y01 z01 .. x0  }{x0 x1  xm} {y0 y1 yn } miny maxy}
  17 #   {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
  18 # 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)
  19
  20 #     { gridequal {
  21
  22 # z00 z01 .. all belong to x=minx and y = miny,.... up y=maxy in n+1 steps
  23 #{ grid {minx maxx} {miny maxy}
  24 #  {{z00 z01 z02 .. z0n } { z10 z11 z12 .. z1n} {..  } ...}
  25 # }
  26 # where a mesh(1) {z00 z01 z11 z10} above
  27
  28
  29
  30 # { mesh {{{x00 y00 z00 } { x01 y01 z01} { x02 y02 z02}  ..}{{x10 y10 z10} {x11 y11 z11} ......} ..}}
  31 # mesh(1) = P00 P01 P11 P10
  32
  33 set sample { variable_grid { 0 1 2 } { 3 4 5} { {21 111 2} {3 4 5 } {6 7 8 }}}
  34 set sample { variable_grid { 0 1 2 } { 3 4 5} { {0 1 2} {3 4 5 } {6 7 8 }}}
  35 set sample { matrix_mesh {{0 1} { 2 3 } {4 5 }}  {{0 1} { 2 3 } {4 5 }}  {{0 1} { 2 3 } {4 5 }} }
  36 set 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 }}}
  37 set sample1 { variable_grid  { 1 2 3 4 5 6 7 8 9 10 }
  38     { 1 2 3 }
  39     {  { 0 0 0 0 0 0 0 0 0 0 }
  40         { 0 0.68404 1.28558 1.73205 1.96962 1.96962 1.73205 1.28558 0.68404 2.44921e-16 }
  41         { 0 1.36808 2.57115 3.4641 3.93923 3.93923 3.4641 2.57115 1.36808 4.89843e-16 }
  42     }
  43 }
  44
  45 set sample { matrix_mesh  {  { 0 0 0 0 0 }
  46     { 1 1 1 1 1 }
  47 }  {  { 0 1 1 0 0 }
  48     { 0 1 1 0 0 }
  49 }  {  { 0 0 1 1 0 }
  50     { 0 0 1 1 0 }
  51 }
  52 }
  53
  54 proc  fixupZ { } {
  55     uplevel 1 {
  56         if { [catch { expr $z + 0 } ] } {
  57             set z nam
  58         }
  59     }
  60 }
  61
  62 proc vectorLength { v } {
  63     expr { sqrt(1.0 * [lindex $v 0]*[lindex $v 0] + [lindex $v 1]*[lindex $v 1] + [lindex $v 2]*[lindex $v 2]) }
  64 }
  65
  66 proc normalizeToLengthOne { v } {
  67     set norm [expr { sqrt(1.0 * [lindex $v 0]*[lindex $v 0] + [lindex $v 1]*[lindex $v 1] + [lindex $v 2]*[lindex $v 2]) }]
  68     if { $norm != 0.0 } {
  69         return [list [expr { [lindex $v 0] / $norm  } ] \
  70                     [expr { [lindex $v 1] / $norm  } ] \
  71                     [expr { [lindex $v 2] / $norm  } ] ]
  72
  73     } else {
  74         return "1.0 0.0 0.0 "
  75     }
  76 }
  77
  78
  79
  80 proc vectorCross { x1 x2 }  {
  81     list \
  82         [expr { [lindex $x1 1]*[lindex $x2 2]- [lindex $x2 1]*[lindex $x1 2]}] \
  83         [expr { [lindex $x1 2]*[lindex $x2 0]- [lindex $x2 2]*[lindex $x1 0] } ] \
  84         [expr { [lindex $x1 0]*[lindex $x2 1]- [lindex $x2 0]*[lindex $x1 1] }]
  85 }
  86
  87 proc linspace { a b n } {
  88     if { $n < 2 } { error [mc "from %s to %s requires at least 2 points" $a $b ] }
  89     set del [expr {($b - $a)*1.0/($n -1)  }]
  90     for { set i 0 } { $i < $n } { incr i } {
  91         lappend ans [expr {$a + $del * $i}]
  92     }
  93     return $ans
  94 }
  95
  96
  97 proc addOnePlot3d { win data } {
  98     upvar #0 plot3dMeshes$win meshes
  99     #puts " adding meshes = plot3dMeshes$win"
 100     #puts "data=$data"
 101     linkLocal $win points zmax zmin zcenter zradius rotationcenter xradius yradius xmin xmax ymin ymax lmesh
 102     makeLocal $win flatten
 103     oset $win colorfun plot3dcolorFun
 104     oset $win cmap c1
 105     global plot3dOptions
 106     catch { unset  meshes }
 107     set points ""
 108
 109     # check whether zradius is a number or "auto"
 110     set dotruncate [expr ![catch {expr {$zradius + 1} }]]
 111     if { $dotruncate } {
 112                 set zzmax [expr {$zcenter + $zradius}]
 113                 set zzmin [expr {$zcenter - $zradius}]
 114     }
 115
 116     set k [llength $points]
 117     set type [lindex $data 0]
 118     # in general the data should be a list of plots..
 119     if { [lsearch {grid mesh variable_grid matrix_mesh }  $type ]>=0 } {
 120         set alldata [list $data]
 121     } else {set alldata $data}
 122     set lmesh {}
 123     set first_mesh 0
 124     foreach data $alldata {
 125         set type [lindex $data 0]
 126         if { "$type" == "grid" } {
 127             desetq "xmin xmax" [lindex $data 1]
 128             desetq "ymin ymax" [lindex $data 2]
 129             set pts [lindex $data 3]
 130
 131             set ncols [llength $pts]
 132             set nrows  [llength [lindex $pts 0]]
 133             set data [list variable_grid [linspace $xmin $xmax $ncols] \
 134                           [linspace $ymin $ymax $nrows] \
 135                           $pts ]
 136             set type "variable_grid"
 137         }
 138         if { "$type" == "variable_grid" } {
 139             set first_mesh [llength $lmesh]
 140             desetq "xrow yrow zmat" [lrange $data 1 end]
 141             # puts "xrow=$xrow,yrow=$yrow,zmat=$zmat"
 142             set nx [expr {[llength $xrow] -1}]
 143             set ny [expr {[llength $yrow] -1}]
 144             #puts "nx=$nx,ny=$ny"
 145             #   set xmin [lindex $xrow 0]
 146             #   set xmax [lindex $xrow $nx]
 147             #   set ymin [lindex $yrow 0]
 148             #   set ymax [lindex $yrow $ny]
 149             desetq "xmin xmax" [minMax $xrow ""]
 150             desetq "ymin ymax" [minMax $yrow ""]
 151             desetq "zmin zmax" [matrixMinMax [list $zmat]]
 152             if { $dotruncate } {
 153                 set zmax  $zzmax
 154                 set zmin  $zzmin
 155             }
 156             #   puts "and now"
 157             #   dshow nx xmin xmax ymin ymax zmin zmax
 158
 159             for {set j 0} { $j <= $ny } { incr j} {
 160                 set y [lindex $yrow $j]
 161                 set row [lindex $zmat $j]
 162                 for {set i 0} { $i <= $nx } { incr i} {
 163                     set x [lindex $xrow $i]
 164                     set z [lindex $row $i]
 165                     #puts "x=$x,y=$y,z=$z, at ($i,$j)"
 166                     fixupZ
 167                     if { $j != $ny && $i != $nx } {
 168                         lappend lmesh [list $k [expr { $k+3 }] \
 169                                            [expr { $k+3+($nx+1)*3 }] \
 170                                            [expr { $k+($nx+1)*3 }]]
 171                     }
 172                     incr k 3
 173                     lappend points $x $y $z
 174                 }
 175             }
 176             setupPlot3dColors $win $first_mesh
 177         } elseif { "$type" == "matrix_mesh" } {
 178             set first_mesh [llength $lmesh]
 179             desetq "xmat ymat zmat" [lrange $data 1 end]
 180             foreach v {x y z} {
 181                 desetq "${v}min ${v}max" [matrixMinMax [list [set ${v}mat]]]
 182             }
 183             if { $dotruncate } {
 184                 set zmax  $zzmax
 185                 set zmin  $zzmin
 186             }
 187             set nj [expr {[llength [lindex $xmat 0]] -1 }]
 188             set ni [expr {[llength $xmat ] -1 }]
 189             set i -1
 190             set k [llength $points]
 191             foreach rowx $xmat rowy $ymat rowz $zmat {
 192                 set j -1
 193                 incr i
 194                 if { [llength $rowx] != [llength $rowy] } {
 195                     error [concat [mc "mismatch"] "rowx:$rowx,rowy:$rowy"]
 196                 }
 197                 if { [llength $rowx] != [llength $rowz] } {
 198                     error [concat [mc "mismatch"] "rowx:$rowx,rowz:$rowz"]
 199                 }
 200                 foreach x $rowx y $rowy z $rowz {
 201                     incr j
 202                     if { $j != $nj && $i != $ni } {
 203                         #puts "tes=($i,$j) $x, $y, $z"
 204                         lappend lmesh [ list \
 205                                             $k [expr { $k+3 } ] [expr { $k + 3  + ($nj+1)*3}] \
 206                                             [expr { $k+($nj+1)*3 }] ]
 207                     }
 208                     incr k 3
 209                     lappend points $x $y $z
 210                 }
 211             }
 212             setupPlot3dColors $win $first_mesh
 213         } elseif { 0 && "$type" == "mesh" } {
 214             set first_mesh [llength $lmesh]
 215             # walk thru compute the xmin, xmax, ymin , ymax...
 216             # and then go thru setting up the mesh array..
 217             # and maybe setting up the color map for these meshes..
 218             #
 219             # { mesh {{{x00 y00 z00 } { x01 y01 z01} { x02 y02 z02}  ..}{{x10 y10 z10} {x11 y11 z11} ......} ..}}
 220             # mesh(1) = P00 P01 P11 P10
 221             set mdata [lindex $data 1]
 222             set nx [llength $mdata]
 223             set ny [llength [lindex $mdata 0]]
 224
 225             for {set i 0} { $i <= $nx } { incr i} {
 226                 set pts [lindex $mdata $i]
 227                 set j 0
 228                 foreach { x y z} $pts {
 229                     fixupZ $z
 230                     if { $j != $ny && $i != $nx } {
 231                         lappend lmesh [list
 232                                        $k [expr { $k+3 }] [expr { $k+3+($ny+1)*3 }] \
 233                                            [expr { $k+($ny+1)*3 }] ]
 234                     }
 235                 }
 236                 incr k 3
 237                 lappend points $x $y $z
 238                 incr j
 239             }
 240             setupPlot3dColors $win $first_mesh
 241         } elseif { "[assq $type $plot3dOptions notthere]" != "notthere" } {
 242             oset $win $type [lindex $data 1]
 243             if { $type == "zradius" } {
 244                 # check whether zradius is a number or "auto"
 245                 set dotruncate [expr ![catch {expr {$zradius + 1} }]]
 246                 if { $dotruncate } {
 247                     set zzmax [expr {$zcenter + $zradius}]
 248                     set zzmin [expr {$zcenter - $zradius}]
 249                 }
 250             }
 251         }
 252         if { $first_mesh != [llength $lmesh] } {
 253             # set up the min/max values for the complete plot
 254             foreach v { x y z } {
 255                 if { [info exists ${v}gmin] } {
 256                     if { [set ${v}min] < [set ${v}gmin] } {
 257                         set ${v}gmin [set ${v}min]
 258                     }
 259                 } else {
 260                     set ${v}gmin [set ${v}min]
 261                 }
 262                 if { [info exists ${v}gmax] } {
 263                     if { [set ${v}max] > [set ${v}gmax] } {
 264                         set ${v}gmax [set ${v}max]
 265                     }
 266                 } else {
 267                     set ${v}gmax [set ${v}max]
 268                 }
 269             }
 270         }
 271     }
 272     foreach v { x y z } {
 273         set ${v}min [set ${v}gmin]
 274         set ${v}max [set ${v}gmax]
 275         set a [set ${v}min]
 276         set b  [set ${v}max]
 277         if { $a == $b } {
 278             set ${v}min [expr {$a -1}]
 279             set ${v}max [expr {$a +1}]
 280         }
 281         set ${v}radius [expr {($b - $a)/2.0}]
 282         set ${v}center [expr {($b + $a)/2.0}]
 283     }
 284     set rotationcenter [list [list [expr {.5*($xmax+$xmin)}]] \
 285                             [list [expr {.5*($ymax+$ymin)}]] \
 286                             [list [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
 293 proc 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
 300 proc 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
 311 proc 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
 318 proc 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 #
 380 proc 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