1
## Automatically adapted for scipy Oct 31, 2005 by
3
# $Id: gistdemo3d.py 1592 2006-02-01 16:49:05Z pearu $
4
# ---------------------------------------------------------------------
8
# PURPOSE: Demonstrate 3D functions in pygist.
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
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
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.
31
from arrayfns import *
32
from numpy.random import *
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."
37
window3 (hcp = "talk.ps", dump = 1)
41
demo5_n = 20 * ones (3)
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
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
61
draw3 ( not making_movie )
65
i = raw_input ("Type in any string to continue; ^C to return to prompt. ")
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.
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
79
run (3) demonstrates the slice3, slice2, and pl3tree functions,
80
as well as changing the orientation of the 3D object
84
if len (itest) == 0 or itest [0] == 1 :
86
x = span (-1, 1, 64, 64)
88
z = (x + y) * exp (-6.*(x*x+y*y))
90
print "Plot wire frame: 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)
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)
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)
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)
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)
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)
131
if len (itest) == 0 or itest [0] == 2 :
134
x = span (-1, 1, 64, 64)
136
z = (x + y) * exp (-6.*(x*x+y*y))
137
print "Default lighting: 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)
145
print "light3(diffuse=.2,specular=1)"
146
light3(diffuse=.2,specular=1)
148
[xmin, xmax, ymin, ymax] = draw3(1) # not necessary interactively
149
limits (xmin, xmax, ymin, ymax)
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))"
156
movie(demo5_light, lims = [xmin, xmax, ymin, ymax])
163
if len (itest) == 0 or itest [0] == 3 :
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
180
print " test uses " + `(nx - 1) * (ny - 1) * (nz - 1)` + " cells"
181
elapsed = [0., 0., 0.]
182
elapsed = timer_ (elapsed)
185
[nv, xyzv, dum] = slice3 (m3, 1, None, None, value = .50)
187
[nw, xyzw, dum] = slice3 (m3, 1, None, None, value = 1.)
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)
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)
204
slice2 (- pyz, nw, xyzw, None)
206
elapsed = timer_ (elapsed)
207
timer_print ("slicing time", elapsed - elapsed0)
210
print 'Generate palette for pl3tree: split_palette ("earth.gp")'
211
split_palette ("earth.gp")
212
print "gnomon -- turn on gnomon"
215
print "pl3tree with 1 slicing plane, 2 isosurfaces"
217
# Make sure we don't draw till ready
219
pl3tree (np, xyzp, vp, pyz)
225
light3 (diffuse = .2, specular = 1)
232
print "spin3 animated rotation, use rot3 or orient3 for one frame"
233
# don't want limits to autoscale during animation
236
limits ( ) # back to autoscaling
245
# .. PR is not available. Comment out this section.
246
# if len (itest) == 0 or itest [0] == 4 :
248
# f = PR ('./bills_plot')
249
# n_nodes = f.NumNodes
250
# n_z = f.NodesOnZones
254
# c = f.ZNodeVelocity
255
# n_zones = f.NumZones
256
# # Put vertices in right order for Gist
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))
264
# # draw a colored plane first
267
# # Make sure we don't draw till ready
269
# [np, xyzp, vp] = slice3 (m3, pyz, None, None, 1)
270
# pl3tree (np, xyzp, vp, pyz, split = 0)
271
# palette ("rainbow.gp")
278
# slice2 (- pyz, nv, xyzv, None)
279
# [nw, xyzw, cw] = slice3 (m3, 1, None, None, 1, value = .9 * min (c) )
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)
290
# print "gnomon -- turn on gnomon"
293
# # Make sure we don't draw till ready
295
# pl3tree (nv, xyzv) # , cv)
296
# pl3tree (nw, xyzw) # , cw)
297
# pl3tree (nvi, xyzvi) # , cvi)
298
# pl3tree (nvj, xyzvj) # , cvi)
300
# light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
302
# palette ("gray.gp")
306
# print "spin3 animated rotation, use rot3 or orient3 for one frame"
307
# # don't want limits to autoscale during animation
309
# limits ( ) # back to autoscaling
315
# palette ("gray.gp")
320
# del nv, xyzv, cv, nw, xyzw, cw, nvi, xyzvi, cvi, nvj, xyzvj, cvj
321
# # Make sure we don't draw till ready
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.)
327
# slice2 (pxz, nv, xyzv, None)
330
# light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
332
# palette ("heat.gp")
336
# limits ( ) # back to autoscaling
342
# if len (itest) == 0 or itest [0] == 5 :
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 )
351
# m3 = mesh3 (nums, dxs, x0s, funcs = [transpose (c)])
352
# [nv, xyzv, dum] = slice3 (m3, 1, None, None, value = 6.5)
355
# print "gnomon -- turn on gnomon"
357
# # Make sure we don't draw till ready
359
# palette ("rainbow.gp")
362
# light3 (diffuse = .2, specular = 1)
369
# if len (itest) == 0 or itest [0] == 6 :
370
# # Try Bill's irregular mesh
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
380
# # Now we need to convert this information to avs-style data
381
# istart = 0 # beginning index into ZLstruct_nodelist
382
# NodeError = "NodeError"
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]))
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]))
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]))
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])))
422
# istart = istart + ZLss [i] * ZLsc [i]
424
# raise NodeError, `ZLss[i]` + "is an incorrect number of nodes."
426
# m3 = mesh3 (x, y, z, funcs = [c], verts = [nz_tet, nz_pyr, nz_prism,
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))
432
# # draw a colored plane first
435
# # Make sure we don't draw till ready
437
# [np, xyzp, vp] = slice3 (m3, pyz, None, None, 1)
438
# pl3tree (np, xyzp, vp, pyz, split = 0)
439
# palette ("rainbow.gp")
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)
462
# print "gnomon -- turn on gnomon"
465
# # Make sure we don't draw till ready
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)
474
# light3 (ambient = 0, diffuse = .5, specular = 1, sdir = [0, 0, -1])
476
# palette ("gray.gp")
479
# palette ("heat.gp")
482
if len (itest) == 0 or itest [0] == 7 :
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
494
# Make sure we don't draw till ready
496
palette ("rainbow.gp")
500
plwf (z, fill = z, ecolor = "black")
501
[xmin, xmax, ymin, ymax] = draw3 (1)
502
limits (xmin, xmax, ymin, ymax)
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)
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)
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)
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)
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)
537
movie(demo5_light, lims = [xmin, xmax, ymin, ymax])
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)
547
limits (xmin, xmax, ymin, ymax)
551
if len (itest) == 0 or itest [0] == 8 :
553
print "Test pl3surf on the sombrero function"
554
# compute sombrero function
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
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.
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)))
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])
581
m3 = mesh3 (nxnynz, dxdydz, x0y0z0, funcs = [meshf])
583
# Make sure we don't draw till ready
586
[nv, xyzv, col] = slice3 (m3, 1, None, None, value = 0.)
590
limits (lim [0], lim [1], 1.5*lim [2], 1.5*lim [3])
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)
599
dif = 0.5 * (lim [3] - lim [2])
600
limits (lim [0], lim [1], lim [2] - dif, lim [3] + dif)
602
palette ("rainbow.gp")
605
print "Try plzcont -- see if smooth mode is possible"
606
print "plzcont (nv, xyzv)"
612
print "plzcont (nv, xyzv, contours = 20)"
613
plzcont (nv, xyzv, contours = 20)
617
print 'plzcont (nv, xyzv, contours = 20, scale = "log")'
618
plzcont (nv, xyzv, contours = 20, scale = "log")
622
print 'plzcont (nv, xyzv, contours = 20, scale = "normal")'
623
plzcont (nv, xyzv, contours = 20, scale = "normal")
627
if len (itest) == 0 or itest [0] == 9 :
635
# The following computations define an interesting 3d surface.
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)
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))
651
ireg [1:15, 12:lmax]=3
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
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)
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)
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])
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)
688
palette ("rainbow.gp")
691
print 'palette ("rainbow.gp")'
692
[nv, xyzv, col] = slice3mesh (nxny, dxdy, x0y0, freg, color = z)
693
pl3surf (nv, xyzv, values = col)
697
print 'palette ("stern.gp")'
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)
708
print 'palette ("gray.gp")'
709
print "light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))"
712
light3 ( diffuse=.1, specular=1., sdir=array([0,0,-1]))