~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/sandbox/xplt/gistdemo3d.py~

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
## Automatically adapted for scipy Oct 31, 2005 by 
 
2
 
 
3
# $Id: gistdemo3d.py 1592 2006-02-01 16:49:05Z pearu $
 
4
#  ---------------------------------------------------------------------
 
5
#
 
6
#  NAME:     gistdemo3d.py
 
7
#
 
8
#  PURPOSE:  Demonstrate 3D functions in pygist.
 
9
#
 
10
#  CHANGES:
 
11
#  12/06/03 mdh Change demo5_n from an array of floats to an array of ints
 
12
#               (no reason to use floats and the change avoids a deprecation
 
13
#               warning.
 
14
#  04/22/03 llc This version is the original demo5 in gist3d with the
 
15
#               tests using the non-existent PR module commented out.
 
16
#               PR is part of the outdated pyPDB.  These tests read in
 
17
#               data files in Portable Data Base (PDB) or Silo format
 
18
#               and were likely not included in the original distribution
 
19
#               due to size.
 
20
#
 
21
#  ---------------------------------------------------------------------
 
22
# Copyright (c) 1996, 1997, The Regents of the University of California.
 
23
# All rights reserved.  See Legal.htm for full text and disclaimer.
 
24
 
 
25
from plwf import *
 
26
from pl3d import *
 
27
from movie import *
 
28
from slice3 import *
 
29
from yorick import *
 
30
from gist import *
 
31
from arrayfns import *
 
32
from numpy.random import * 
 
33
 
 
34
print "Type gistdemo3d.run() or gistdemo3d.run(i), i = 1, 2, or 3."
 
35
print "Partway, plots are written to Postscript file; see talk.ps."
 
36
 
 
37
window3 (hcp = "talk.ps", dump = 1)
 
38
 
 
39
palette ("gray.gp")
 
40
 
 
41
demo5_n = 20 * ones (3)
 
42
 
 
43
making_movie = 0
 
44
 
 
45
def demo5_light (i) :
 
46
   global making_movie
 
47
   if i >= 30 : return 0
 
48
   theta = pi / 4 + (i - 1) * 2 * pi/29
 
49
   light3 (sdir = array ( [cos(theta), .25, sin(theta)], Float))
 
50
   # without an explicit call to draw3, the light3 function would
 
51
   # cause no changes until Python paused for input from the keyboard,
 
52
   # since unlike the primitive plotting functions (plg, plf, plfp, ...)
 
53
   # the fma call made by the movie function will not trigger the
 
54
   # 3D display list
 
55
   # any movie frame display function which uses the 3D drawing
 
56
   # functions in pl3d.i will need to do this
 
57
   # the !making_movie flag supresses the fma in draw3 if this function
 
58
   # is called by movie (which issues its own fma), but allows it
 
59
   # otherwise 
 
60
 
 
61
   draw3 ( not making_movie )
 
62
   return 1
 
63
 
 
64
def paws ( ) :
 
65
    i = raw_input ("Type in any string to continue; ^C to return to prompt. ")
 
66
    return
 
67
 
 
68
def run (*itest) :
 
69
   """run() or run(i)
 
70
     Run examples of use of pl3d.i, plwf.i, and slice3.i.  With
 
71
     argument I = 1, 2, or 3, run that particular demonstration.
 
72
     Read the source code to understand the details of how the
 
73
     various effects are obtained.
 
74
 
 
75
     run (1) demonstrates the various effects which can be obtained
 
76
     with the plwf (plot wire frame) function.
 
77
     run (2) demonstrates shading effects controlled by the light3
 
78
     function
 
79
     run (3) demonstrates the slice3, slice2, and pl3tree functions,
 
80
     as well as changing the orientation of the 3D object
 
81
   """
 
82
   global making_movie
 
83
 
 
84
   if len (itest) == 0 or itest [0] == 1 :
 
85
      set_draw3_ (0)
 
86
      x = span (-1, 1, 64, 64)
 
87
      y = transpose (x)
 
88
      z = (x + y) * exp (-6.*(x*x+y*y))
 
89
      limits_(square = 1)
 
90
      print "Plot wire frame: plwf (z, y, x)"
 
91
      orient3 ( )
 
92
      light3 ( )
 
93
      plwf (z, y, x)
 
94
      [xmin, xmax, ymin, ymax] = draw3(1) # not necessary interactively
 
95
      limits (xmin, xmax, ymin, ymax)
 
96
      plt("opaque wire mesh", .30, .42)
 
97
      paws ( )
 
98
 
 
99
      print 'plwf(z,y,x,shade=1,ecolor="red")'
 
100
      plwf(z,y,x,shade=1,ecolor="red")
 
101
      [xmin, xmax, ymin, ymax] = draw3(1) # not necessary interactively
 
102
      limits (xmin, xmax, ymin, ymax)
 
103
      paws()
 
104
 
 
105
      print "plwf(z,y,x,shade=1,edges=0)"
 
106
      plwf(z,y,x,shade=1,edges=0)
 
107
      [xmin, xmax, ymin, ymax] = draw3(1) # not necessary interactively
 
108
      limits (xmin, xmax, ymin, ymax)
 
109
      paws ( )
 
110
 
 
111
      print "light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))"
 
112
      light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))
 
113
      [xmin, xmax, ymin, ymax] = draw3(1)
 
114
      limits (xmin, xmax, ymin, ymax)
 
115
      paws ( )
 
116
 
 
117
      print "light3 ( diffuse=.5, specular=1., sdir=array([1,.5,1]))"
 
118
      light3 ( diffuse=.5, specular=1., sdir=array([1,.5,1]))
 
119
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
120
      limits (xmin, xmax, ymin, ymax)
 
121
      paws ( )
 
122
 
 
123
      print "light3 ( ambient=.1,diffuse=.1,specular=1.,"
 
124
      print "   sdir=array([[0,0,-1],[1,.5,1]]),spower=array([4,2]))"
 
125
      light3 ( ambient=.1,diffuse=.1,specular=1.,
 
126
             sdir=array([[0,0,-1],[1,.5,1]]),spower=array([4,2]))
 
127
      [xmin, xmax, ymin, ymax] = draw3(1)
 
128
      limits (xmin, xmax, ymin, ymax)
 
129
      paws ( )
 
130
 
 
131
   if len (itest) == 0 or itest [0] == 2 :
 
132
 
 
133
      set_draw3_ (0)
 
134
      x = span (-1, 1, 64, 64)
 
135
      y = transpose (x)
 
136
      z = (x + y) * exp (-6.*(x*x+y*y))
 
137
      print "Default lighting:  light3()"
 
138
      orient3 ( )
 
139
      light3 ( )
 
140
      plwf (z,y,x,shade=1,edges=0)
 
141
      [xmin, xmax, ymin, ymax] = draw3 (1) # not necessary interactively
 
142
      limits (xmin, xmax, ymin, ymax)
 
143
      paws( )
 
144
 
 
145
      print "light3(diffuse=.2,specular=1)"
 
146
      light3(diffuse=.2,specular=1)
 
147
      limits_(square = 1)
 
148
      [xmin, xmax, ymin, ymax] = draw3(1) # not necessary interactively
 
149
      limits (xmin, xmax, ymin, ymax)
 
150
      paws()
 
151
 
 
152
      print "movie(demo5_light, lims = [xmin, xmax, ymin, ymax])"
 
153
      print "demo5_light calls:"
 
154
      print "light3 (sdir = array ( [cos(theta), .25, sin(theta)], Float))"
 
155
      making_movie = 1
 
156
      movie(demo5_light, lims = [xmin, xmax, ymin, ymax])
 
157
      making_movie = 0
 
158
      fma()
 
159
      demo5_light(1)
 
160
      paws()
 
161
      light3()
 
162
 
 
163
   if len (itest) == 0 or itest [0] == 3 :
 
164
 
 
165
      nx = demo5_n [0]
 
166
      ny = demo5_n [1]
 
167
      nz = demo5_n [2]
 
168
      xyz = zeros ( (3, nx, ny, nz), Float)
 
169
      xyz [0] = multiply.outer ( span (-1, 1, nx), ones ( (ny, nz), Float))
 
170
      xyz [1] = multiply.outer ( ones (nx, Float),
 
171
         multiply.outer ( span (-1, 1, ny), ones (nz, Float)))
 
172
      xyz [2] = multiply.outer ( ones ( (nx, ny), Float), span (-1, 1, nz))
 
173
      r = sqrt (xyz [0] ** 2 + xyz [1] **2 + xyz [2] **2)
 
174
      theta = arccos (xyz [2] / r)
 
175
      phi = arctan2 (xyz [1] , xyz [0] + logical_not (r))
 
176
      y32 = sin (theta) ** 2 * cos (theta) * cos (2 * phi)
 
177
      m3 = mesh3 (xyz, funcs = [r * (1. + y32)])
 
178
      del r, theta, phi, xyz, y32
 
179
 
 
180
      print "   test uses " + `(nx - 1) * (ny - 1) * (nz - 1)` + " cells"
 
181
      elapsed = [0., 0., 0.]
 
182
      elapsed = timer_ (elapsed)
 
183
      elapsed0 = elapsed
 
184
 
 
185
      [nv, xyzv, dum] = slice3 (m3, 1, None, None, value = .50)
 
186
          # (inner isosurface)
 
187
      [nw, xyzw, dum] = slice3 (m3, 1, None, None, value = 1.)
 
188
          # (outer isosurface)
 
189
      pxy = plane3 ( array ([0, 0, 1], Float ), zeros (3, Float))
 
190
      pyz = plane3 ( array ([1, 0, 0], Float ), zeros (3, Float))
 
191
      [np, xyzp, vp] = slice3 (m3, pyz, None, None, 1)
 
192
          # (pseudo-colored slice)
 
193
      [np, xyzp, vp] = slice2 (pxy, np, xyzp, vp)
 
194
          # (cut slice in half)
 
195
      [nv, xyzv, d1, nvb, xyzvb, d2] = \
 
196
          slice2x (pxy, nv, xyzv, None)
 
197
      [nv, xyzv, d1] = \
 
198
          slice2 (- pyz, nv, xyzv, None)
 
199
          # (...halve one of those halves)
 
200
      [nw, xyzw, d1, nwb, xyzwb, d2] = \
 
201
          slice2x ( pxy , nw, xyzw, None)
 
202
          # (split outer in halves)
 
203
      [nw, xyzw, d1] = \
 
204
          slice2 (- pyz, nw, xyzw, None)
 
205
 
 
206
      elapsed = timer_ (elapsed)
 
207
      timer_print ("slicing time", elapsed - elapsed0)
 
208
 
 
209
      fma ()
 
210
      print 'Generate palette for pl3tree:  split_palette ("earth.gp")'
 
211
      split_palette ("earth.gp")
 
212
      print "gnomon -- turn on gnomon"
 
213
      gnomon (1)
 
214
 
 
215
      print "pl3tree with 1 slicing plane, 2 isosurfaces"
 
216
      clear3 ()
 
217
      # Make sure we don't draw till ready
 
218
      set_draw3_ (0)
 
219
      pl3tree (np, xyzp, vp, pyz)
 
220
      pl3tree (nvb, xyzvb)
 
221
      pl3tree (nwb, xyzwb)
 
222
      pl3tree (nv, xyzv)
 
223
      pl3tree (nw, xyzw)
 
224
      orient3 ()
 
225
      light3 (diffuse = .2, specular = 1)
 
226
      limits ()
 
227
      limits (square=1)
 
228
      demo5_light (1)
 
229
      paws ()
 
230
      hcp ()
 
231
 
 
232
      print "spin3 animated rotation, use rot3 or orient3 for one frame"
 
233
      # don't want limits to autoscale during animation
 
234
      lims = limits ( )
 
235
      spin3 ()
 
236
      limits ( ) # back to autoscaling
 
237
      demo5_light (1)
 
238
      paws ()
 
239
 
 
240
      light3 ()
 
241
      gnomon (0)
 
242
      limits (square = 1)
 
243
      palette ("gray.gp")
 
244
 
 
245
#  .. PR is not available.  Comment out this section.
 
246
#   if len (itest) == 0 or itest [0] == 4 :
 
247
#      from PR import *
 
248
#      f = PR ('./bills_plot')
 
249
#      n_nodes = f.NumNodes
 
250
#      n_z = f.NodesOnZones
 
251
#      x = f.XNodeCoords
 
252
#      y = f.YNodeCoords
 
253
#      z = f.ZNodeCoords
 
254
#      c = f.ZNodeVelocity
 
255
#      n_zones = f.NumZones
 
256
#      # Put vertices in right order for Gist
 
257
#      n_z = transpose (
 
258
#         take (transpose (n_z), array ( [0, 4, 3, 7, 1, 5, 2, 6])))
 
259
#      m3 = mesh3 (x, y, z, funcs = [c], verts = n_z ) # [0:10])
 
260
#      [nv, xyzv, cv] = slice3 (m3, 1, None, None, 1, value = .9 * max (c) )
 
261
#      pyz = plane3 ( array ([1, 0, 0], Float ), zeros (3, Float))
 
262
#      pxz = plane3 ( array ([0, 1, 0], Float ), zeros (3, Float))
 
263
#
 
264
#      # draw a colored plane first
 
265
#      fma ()
 
266
#      clear3 ()
 
267
#      # Make sure we don't draw till ready
 
268
#      set_draw3_ (0)
 
269
#      [np, xyzp, vp] = slice3 (m3, pyz, None, None, 1)
 
270
#      pl3tree (np, xyzp, vp, pyz, split = 0)
 
271
#      palette ("rainbow.gp")
 
272
#      orient3 ()
 
273
#      demo5_light (1)
 
274
#      paws ()
 
275
#      
 
276
#
 
277
#     [nv, xyzv, d1] = \
 
278
#         slice2 (- pyz, nv, xyzv, None)
 
279
#      [nw, xyzw, cw] = slice3 (m3, 1, None, None, 1, value = .9 * min (c) )
 
280
#     [nw, xyzw, d1] = \
 
281
#         slice2 (- pyz, nw, xyzw, None)
 
282
#      [nvi, xyzvi, cvi] = slice3 (m3, 1, None, None, 1, value = .5 * min (c) )
 
283
#      [nvi, xyzvi, cvi] = \
 
284
#          slice2 (- pyz, nvi, xyzvi, cvi)
 
285
#      [nvj, xyzvj, cvj] = slice3 (m3, 1, None, None, 1, value = .5 * max (c) )
 
286
#      [nvj, xyzvj, cvj] = \
 
287
#          slice2 (- pyz, nvj, xyzvj, cvj)
 
288
#
 
289
#      fma ()
 
290
#      print "gnomon -- turn on gnomon"
 
291
#      gnomon (1)
 
292
#      clear3 ()
 
293
#      # Make sure we don't draw till ready
 
294
#      set_draw3_ (0)
 
295
#      pl3tree (nv, xyzv) # , cv)
 
296
#      pl3tree (nw, xyzw) # , cw)
 
297
#      pl3tree (nvi, xyzvi) # , cvi)
 
298
#      pl3tree (nvj, xyzvj) # , cvi)
 
299
#      orient3 ()
 
300
#      light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
 
301
#      limits (square=1)
 
302
#      palette ("gray.gp")
 
303
#      demo5_light (1)
 
304
#      paws ()
 
305
#      
 
306
#      print "spin3 animated rotation, use rot3 or orient3 for one frame"
 
307
#      # don't want limits to autoscale during animation
 
308
#      spin3 ()
 
309
#      limits ( ) # back to autoscaling
 
310
#      demo5_light (1)
 
311
#      paws ()
 
312
#
 
313
#      light3 ()
 
314
#      gnomon (0)
 
315
#      palette ("gray.gp")
 
316
#
 
317
#      draw3 ( 1 )
 
318
#      paws ()
 
319
#      clear3 ()
 
320
#      del nv, xyzv, cv, nw, xyzw, cw, nvi, xyzvi, cvi, nvj, xyzvj, cvj
 
321
#      # Make sure we don't draw till ready
 
322
#      set_draw3_ (0)
 
323
#      for i in range (8) :
 
324
#         [nv, xyzv, cv] = slice3 (m3, 1, None, None, 1, value = .9 * min (c) +
 
325
#             i * (.9 * max (c) - .9 * min (c)) / 8.)
 
326
#         [nv, xyzv, d1] = \
 
327
#             slice2 (pxz, nv, xyzv, None)
 
328
#         pl3tree (nv, xyzv)
 
329
#      orient3 ()
 
330
#      light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
 
331
#      limits (square=1)
 
332
#      palette ("heat.gp")
 
333
#      demo5_light (1)
 
334
#      paws ()
 
335
#      spin3 ()
 
336
#      limits ( ) # back to autoscaling
 
337
#      demo5_light (1)
 
338
#      paws ()
 
339
#      demo5_light (1)
 
340
#      paws ()
 
341
#
 
342
#   if len (itest) == 0 or itest [0] == 5 :
 
343
#      # Try bert's data
 
344
#      from PR import PR
 
345
#      f = PR ('./berts_plot')
 
346
#      nums = array ( [63, 63, 49], Int)
 
347
#      dxs = array ( [2.5, 2.5, 10.], Float )
 
348
#      x0s = array ( [-80., -80., 0.0], Float )
 
349
#      c = f.c
 
350
#
 
351
#      m3 = mesh3 (nums, dxs, x0s, funcs = [transpose (c)])
 
352
#      [nv, xyzv, dum] = slice3 (m3, 1, None, None, value = 6.5)
 
353
#      fma ()
 
354
#      clear3 ()
 
355
#      print "gnomon -- turn on gnomon"
 
356
#      gnomon (1)
 
357
#      # Make sure we don't draw till ready
 
358
#      set_draw3_ (0)
 
359
#      palette ("rainbow.gp")
 
360
#      pl3tree (nv, xyzv)
 
361
#      orient3 ()
 
362
#      light3 (diffuse = .2, specular = 1)
 
363
#      limits (square=1)
 
364
#      demo5_light (1)
 
365
#      paws ()
 
366
#      spin3 ()
 
367
#      demo5_light (1)
 
368
#      paws ()
 
369
#   if len (itest) == 0 or itest [0] == 6 :
 
370
#      # Try Bill's irregular mesh
 
371
#      from PR import PR
 
372
#      f = PR ("ball.s0001")
 
373
#      ZLss = f.ZLstruct_shapesize
 
374
#      ZLsc = f.ZLstruct_shapecnt
 
375
#      ZLsn = f.ZLstruct_nodelist
 
376
#      x = f.sap_mesh_coord0
 
377
#      y = f.sap_mesh_coord1
 
378
#      z = f.sap_mesh_coord2
 
379
#      c = f.W_vel_data
 
380
#      # Now we need to convert this information to avs-style data
 
381
#      istart = 0 # beginning index into ZLstruct_nodelist
 
382
#      NodeError = "NodeError"
 
383
#      ntet = 0
 
384
#      nhex = 0
 
385
#      npyr = 0
 
386
#      nprism = 0
 
387
#      nz_tet = []
 
388
#      nz_hex = []
 
389
#      nz_pyr = []
 
390
#      nz_prism = []
 
391
#      for i in range (4) :
 
392
#         if ZLss [i] == 4 : # TETRAHEDRON
 
393
#            nz_tet = reshape (ZLsn [istart: istart + ZLss [i] * ZLsc [i]],
 
394
#                     (ZLsc [i], ZLss [i]))
 
395
#            ntet = ZLsc [i]
 
396
#            istart = istart + ZLss [i] * ZLsc [i]
 
397
#         elif ZLss[i] == 5 : # PYRAMID
 
398
#            nz_pyr = reshape (ZLsn [istart: istart + ZLss [i] * ZLsc [i]],
 
399
#                     (ZLsc [i], ZLss [i]))
 
400
#            npyr = ZLsc [i]
 
401
#            # Now reorder the points (bill has the apex last instead of first)
 
402
#            nz_pyr = transpose (
 
403
#               take (transpose (nz_pyr), array ( [4, 0, 1, 2, 3])))
 
404
#            istart = istart + ZLss [i] * ZLsc [i]
 
405
#         elif ZLss[i] == 6 : # PRISM
 
406
#            nz_prism = reshape (ZLsn [istart: istart + ZLss [i] * ZLsc [i]],
 
407
#                     (ZLsc [i], ZLss [i]))
 
408
#            nprism = ZLsc [i]
 
409
#            # now reorder the points (bill goes around a square face 
 
410
#            # instead of traversing the opposite sides in the same direction.
 
411
#            nz_prism = transpose (
 
412
#               take (transpose (nz_prism), array ( [0, 1, 3, 2, 4, 5])))
 
413
#            istart = istart + ZLss [i] * ZLsc [i]
 
414
#         elif ZLss[i] == 8 : # HEXAHEDRON
 
415
#            nz_hex = reshape (ZLsn [istart: istart + ZLss [i] * ZLsc [i]],
 
416
#                     (ZLsc [i], ZLss [i]))
 
417
#            # now reorder the points (bill goes around a square face 
 
418
#            # instead of traversing the opposite sides in the same direction.
 
419
#            nz_hex = transpose (
 
420
#               take (transpose (nz_hex), array ( [0, 1, 3, 2, 4, 5, 7, 6])))
 
421
#            nhex = ZLsc [i]
 
422
#            istart = istart + ZLss [i] * ZLsc [i]
 
423
#         else :
 
424
#            raise NodeError, `ZLss[i]` + "is an incorrect number of nodes."
 
425
 
 
426
#      m3 = mesh3 (x, y, z, funcs = [c], verts = [nz_tet, nz_pyr, nz_prism,
 
427
#         nz_hex])
 
428
#      [nv, xyzv, cv] = slice3 (m3, 1, None, None, 1, value = .9 * max (c) )
 
429
#      pyz = plane3 ( array ([1, 0, 0], Float ), zeros (3, Float))
 
430
#      pxz = plane3 ( array ([0, 1, 0], Float ), zeros (3, Float))
 
431
#
 
432
#      # draw a colored plane first
 
433
#      fma ()
 
434
#      clear3 ()
 
435
#      # Make sure we don't draw till ready
 
436
#      set_draw3_ (0)
 
437
#      [np, xyzp, vp] = slice3 (m3, pyz, None, None, 1)
 
438
#      pl3tree (np, xyzp, vp, pyz, split = 0)
 
439
#      palette ("rainbow.gp")
 
440
#      orient3 ()
 
441
#      limits (square=1)
 
442
#      demo5_light (1)
 
443
#      paws ()
 
444
#
 
445
#      [nw, xyzw, cw] = slice3 (m3, 1, None, None, 1, value = .9 * min (c) )
 
446
#      [nvi, xyzvi, cvi] = slice3 (m3, 1, None, None, 1, value = .1 * min (c) )
 
447
#      [nvi, xyzvi, cvi] = \
 
448
#          slice2 (- pyz, nvi, xyzvi, cvi)
 
449
#      [nvj, xyzvj, cvj] = slice3 (m3, 1, None, None, 1, value = .1 * max (c) )
 
450
#      [nvj, xyzvj, cvj] = \
 
451
#          slice2 (- pyz, nvj, xyzvj, cvj)
 
452
#      [nvii, xyzvii, cvii] = slice3 (m3, 1, None, None, 1,
 
453
#         value = 1.e-12 * min (c) )
 
454
#      [nvii, xyzvii, cvii] = \
 
455
#          slice2 (- pyz, nvii, xyzvii, cvii)
 
456
#      [nvjj, xyzvjj, cvjj] = slice3 (m3, 1, None, None, 1,
 
457
#         value = 1.e-12 * max (c) )
 
458
#      [nvjj, xyzvjj, cvjj] = \
 
459
#          slice2 (- pyz, nvjj, xyzvjj, cvjj)
 
460
#
 
461
#      fma ()
 
462
#      print "gnomon -- turn on gnomon"
 
463
#      gnomon (1)
 
464
#      clear3 ()
 
465
#      # Make sure we don't draw till ready
 
466
#      set_draw3_ (0)
 
467
#      pl3tree (nv, xyzv) # , cv)
 
468
#      pl3tree (nw, xyzw) # , cw)
 
469
#      pl3tree (nvi, xyzvi) # , cvi)
 
470
#      pl3tree (nvj, xyzvj) # , cvj)
 
471
#      pl3tree (nvii, xyzvii) # , cvii)
 
472
#      pl3tree (nvjj, xyzvjj) # , cvjj)
 
473
#      orient3 ()
 
474
#      light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
 
475
#      limits (square=1)
 
476
#      palette ("gray.gp")
 
477
#      demo5_light (1)
 
478
#      paws ()
 
479
#      palette ("heat.gp")
 
480
#      paws ()
 
481
 
 
482
   if len (itest) == 0 or itest [0] == 7 :
 
483
 
 
484
      print "Test plwf on the sombrero function"
 
485
      # compute sombrero function
 
486
      x = arange (-20, 21, dtype = Float)
 
487
      y = arange (-20, 21, dtype = Float)
 
488
      z = zeros ( (41, 41), Float)
 
489
      r = sqrt (add.outer ( x ** 2, y **2)) + 1e-6
 
490
      z = sin (r) / r
 
491
      fma ()
 
492
      clear3 ()
 
493
      gnomon (0)
 
494
      # Make sure we don't draw till ready
 
495
      set_draw3_ (0)
 
496
      palette ("rainbow.gp")
 
497
      limits (square=1)
 
498
      orient3 ()
 
499
      light3 ()
 
500
      plwf (z, fill = z, ecolor = "black")
 
501
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
502
      limits (xmin, xmax, ymin, ymax)
 
503
      paws ()
 
504
 
 
505
      print "Try smooth contours, log mode:"
 
506
      print 'plzcont (nv, xyzv, contours = 20, scale = "normal")'
 
507
      [nv, xyzv, dum] = slice3mesh (x, y, z)
 
508
      zmult = max (max (abs (x)), max (abs (y)))
 
509
      plzcont (nv, xyzv, contours = 20, scale = "normal")
 
510
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
511
      limits (xmin, xmax, ymin, ymax)
 
512
      paws ()
 
513
 
 
514
      print 'plzcont (nv, xyzv, contours = 20, scale = "lin", edges=1)'
 
515
      plzcont (nv, xyzv, contours = 20, scale = "lin", edges=1)
 
516
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
517
      limits (xmin, xmax, ymin, ymax)
 
518
      paws ()
 
519
 
 
520
      print 'plwf (z, fill = z, shade = 1, ecolor = "black")'
 
521
      plwf (z, fill = z, shade = 1, ecolor = "black")
 
522
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
523
      limits (xmin, xmax, ymin, ymax)
 
524
      paws ()
 
525
 
 
526
      print "plwf (z, fill = z, shade = 1, edges = 0)"
 
527
      plwf (z, fill = z, shade = 1, edges = 0)
 
528
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
529
      limits (xmin, xmax, ymin, ymax)
 
530
      paws ()
 
531
 
 
532
      print "light3(diffuse=.2,specular=1)"
 
533
      print "demo5_light calls:"
 
534
      print "light3 (sdir = array ( [cos(theta), .25, sin(theta)], Float))"
 
535
      light3(diffuse=.2,specular=1)
 
536
      making_movie = 1
 
537
      movie(demo5_light, lims = [xmin, xmax, ymin, ymax])
 
538
      making_movie = 0
 
539
      fma()
 
540
      demo5_light(1)
 
541
      paws ()
 
542
 
 
543
      print "plwf (z, fill = None, shade = 1, edges = 0)"
 
544
      plwf (z, fill = None, shade = 1, edges = 0)
 
545
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
546
      palette("gray.gp")
 
547
      limits (xmin, xmax, ymin, ymax)
 
548
      paws ()
 
549
 
 
550
 
 
551
   if len (itest) == 0 or itest [0] == 8 :
 
552
 
 
553
      print "Test pl3surf on the sombrero function"
 
554
      # compute sombrero function
 
555
      nc1 = 100
 
556
      nv1 = nc1 + 1
 
557
      br = - (nc1 / 2)
 
558
      tr = nc1 / 2 + 1
 
559
      x = arange (br, tr, dtype = Float) * 40. / nc1
 
560
      y = arange (br, tr, dtype = Float) * 40. / nc1
 
561
      z = zeros ( (nv1, nv1), Float)
 
562
      r = sqrt (add.outer ( x ** 2, y **2)) + 1e-6
 
563
      z = sin (r) / r
 
564
      # In order to use pl3surf, we need to construct a mesh
 
565
      # using mesh3. The way I am going to do that is to define
 
566
      # a function on the 3d mesh so that the sombrero function
 
567
      # is its 0-isosurface.
 
568
      z0 = min (ravel (z))
 
569
      z0 = z0 - .05 * abs (z0)
 
570
      maxz = max (ravel (z))
 
571
      maxz = maxz + .05 * abs (maxz)
 
572
      zmult = max (max (abs (x)), max (abs (y)))
 
573
      dz = (maxz - z0)
 
574
      nxnynz = array ( [nc1, nc1, 1], Int)
 
575
      dxdydz = array ( [1.0, 1.0, zmult*dz], Float )
 
576
      x0y0z0 = array ( [float (br), float (br), z0*zmult], Float )
 
577
      meshf = zeros ( (nv1, nv1, 2), Float )
 
578
      meshf [:, :, 0] = zmult*z - (x0y0z0 [2])
 
579
      meshf [:, :, 1] = zmult*z - (x0y0z0 [2] + dxdydz [2])
 
580
 
 
581
      m3 = mesh3 (nxnynz, dxdydz, x0y0z0, funcs = [meshf])
 
582
      fma ()
 
583
      # Make sure we don't draw till ready
 
584
      set_draw3_ (0)
 
585
      pldefault(edges=0)
 
586
      [nv, xyzv, col] = slice3 (m3, 1, None, None, value = 0.)
 
587
      orient3 ()
 
588
      pl3surf (nv, xyzv)
 
589
      lim = draw3 (1)
 
590
      limits (lim [0], lim [1], 1.5*lim [2], 1.5*lim [3])
 
591
      palette ("gray.gp")
 
592
      paws ()
 
593
 
 
594
      print "Try new slicing function (slice3mesh) to get color graph"
 
595
      [nv, xyzv, col] = slice3mesh (nxnynz [0:2], dxdydz [0:2], x0y0z0 [0:2],
 
596
         zmult * z, color = zmult * z)
 
597
      pl3surf (nv, xyzv, values = col)
 
598
      lim = draw3 (1)
 
599
      dif = 0.5 * (lim [3] - lim [2])
 
600
      limits (lim [0], lim [1], lim [2] - dif, lim [3] + dif)
 
601
 
 
602
      palette ("rainbow.gp")
 
603
      paws ()
 
604
 
 
605
      print "Try plzcont -- see if smooth mode is possible"
 
606
      print "plzcont (nv, xyzv)"
 
607
      palette ("heat.gp")
 
608
      plzcont (nv, xyzv)
 
609
      draw3 (1)
 
610
      paws ()
 
611
 
 
612
      print "plzcont (nv, xyzv, contours = 20)"
 
613
      plzcont (nv, xyzv, contours = 20)
 
614
      draw3 (1)
 
615
      paws ()
 
616
 
 
617
      print 'plzcont (nv, xyzv, contours = 20, scale = "log")'
 
618
      plzcont (nv, xyzv, contours = 20, scale = "log")
 
619
      draw3(1)
 
620
      paws ()
 
621
 
 
622
      print 'plzcont (nv, xyzv, contours = 20, scale = "normal")'
 
623
      plzcont (nv, xyzv, contours = 20, scale = "normal")
 
624
      draw3(1)
 
625
      paws ()
 
626
 
 
627
   if len (itest) == 0 or itest [0] == 9 :
 
628
 
 
629
      vsf = 0.
 
630
      c = 1
 
631
      s = 1000.
 
632
      kmax = 25
 
633
      lmax = 35
 
634
 
 
635
      # The following computations define an interesting 3d surface.
 
636
 
 
637
      xr = multiply.outer (
 
638
         arange (1, kmax + 1, dtype = Float), ones (lmax, Float))
 
639
      yr = multiply.outer (
 
640
         ones (kmax, Float), arange (1, lmax + 1, dtype = Float))
 
641
      zt = 5. + xr + .2 * random_sample (kmax, lmax)   # ranf (xr)
 
642
      rt = 100. + yr + .2 * random_sample (kmax, lmax)   # ranf (yr)
 
643
      z = s * (rt + zt)
 
644
      z = z + .02 * z * random_sample (kmax, lmax)   # ranf (z)
 
645
      ut = rt/sqrt (rt ** 2 + zt ** 2)
 
646
      vt = zt/sqrt (rt ** 2 + zt ** 2)
 
647
      ireg =  multiply.outer ( ones (kmax, Float), ones (lmax, Float))
 
648
      ireg [0:1, 0:lmax]=0
 
649
      ireg [0:kmax, 0:1]=0
 
650
      ireg [1:15, 7:12]=2
 
651
      ireg [1:15, 12:lmax]=3
 
652
      ireg [3:7, 3:7]=0
 
653
      freg=ireg + .2 * (1. - random_sample (kmax, lmax))  # ranf (ireg))
 
654
      freg=array (freg, Float)
 
655
      #rt [4:6, 4:6] = -1.e8
 
656
      z [3:10, 3:12] = z [3:10, 3:12] * .9
 
657
      z [5, 5] = z [5, 5] * .9
 
658
      z [17:22, 15:18] = z [17:22, 15:18] * 1.2
 
659
      z [16, 16] = z [16, 16] * 1.1
 
660
      orient3 ()
 
661
      print "plwf (freg, shade = 1, edges = 0)"
 
662
      plwf (freg, shade = 1, edges = 0)
 
663
      [xmin, xmax, ymin, ymax] = draw3 (1)
 
664
      limits (xmin, xmax, ymin, ymax)
 
665
      paws ()
 
666
 
 
667
      print "two slice3mesh and two pl3tree"
 
668
      nxny = array ( [kmax - 1, lmax - 1])
 
669
      x0y0 = array ( [0., 0.])
 
670
      dxdy = array ( [1., 1.])
 
671
      [nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, freg)
 
672
      [nw, xyzw, col] = slice3mesh (nxny, dxdy, x0y0, freg + ut)
 
673
      pl3tree (nv, xyzv)
 
674
      pl3tree (nw, xyzw)
 
675
      draw3 (1)
 
676
      limits ( )
 
677
      paws ()
 
678
 
 
679
      print "light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])"
 
680
      light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
 
681
      demo5_light (1)
 
682
      paws ()
 
683
 
 
684
      print "[nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, freg, color = freg)"
 
685
      [nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, freg, color = freg)
 
686
      pl3surf (nv, xyzv, values = col)
 
687
      draw3 (1)
 
688
      palette ("rainbow.gp")
 
689
      paws ()
 
690
 
 
691
      print 'palette ("rainbow.gp")'
 
692
      [nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, freg, color = z)
 
693
      pl3surf (nv, xyzv, values = col)
 
694
      draw3 (1)
 
695
      paws ()
 
696
 
 
697
      print 'palette ("stern.gp")'
 
698
      palette ("stern.gp")
 
699
      paws ()
 
700
 
 
701
      print "[nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, z, color = z)"
 
702
      [nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, z, color = z)
 
703
      pl3surf (nv, xyzv, values = col)
 
704
      orient3(phi=0,theta=0)
 
705
      draw3 (1)
 
706
      paws ()
 
707
 
 
708
      print 'palette ("gray.gp")'
 
709
      print "light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))"
 
710
      set_draw3_ (0)
 
711
      palette ("gray.gp")
 
712
      light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))
 
713
      pl3surf (nv, xyzv)
 
714
      draw3 (1)
 
715
      paws ()
 
716
 
 
717
#     spin3 ()
 
718
#     paws ()
 
719
 
 
720
   hcp_finish ()