1# x08.py PLplot demo for Python
2#
3# Copyright (C) 2004-2016  Alan W. Irwin
4#
5# This file is part of PLplot.
6#
7# PLplot is free software; you can redistribute it and/or modify
8# it under the terms of the GNU Library General Public License as published
9# by the Free Software Foundation; either version 2 of the License, or
10# (at your option) any later version.
11#
12# PLplot is distributed in the hope that it will be useful,
13# but WITHOUT ANY WARRANTY; without even the implied warranty of
14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15# GNU Library General Public License for more details.
16#
17# You should have received a copy of the GNU Library General Public License
18# along with PLplot; if not, write to the Free Software
19# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20
21from numpy import *
22
23# These values must be odd, for the middle
24# of the index range to be an integer, and thus
25# to correspond to the exact floating point centre
26# of the sombrero.
27
28XPTS = 35               # Data points in x
29YPTS = 45               # Data points in y
30
31alt = [60.0, 40.0]
32
33az = [30.0, -30.0]
34
35title = ["#frPLplot Example 8 - Alt=60, Az=30",
36         "#frPLplot Example 8 - Alt=40, Az=-30"]
37
38# Routine for defining a specific color map 1 in HLS space.
39# if gray is true, use basic grayscale variation from half-dark to light.
40# otherwise use false color variation from blue (240 deg) to red (360 deg).
41def cmap1_init(w, gray):
42    # Independent variable of control points.
43    i = array((0., 1.))
44    if gray:
45        # Hue for control points.  Doesn't matter since saturation is zero.
46        h = array((0., 0.))
47        # Lightness ranging from half-dark (for interest) to light.
48        l = array((0.5, 1.))
49        # Gray scale has zero saturation
50        s = array((0., 0.))
51    else:
52        # Hue ranges from blue (240 deg) to red (0 or 360 deg)
53        h = array((240., 0.))
54        # Lightness and saturation are constant (values taken from C example).
55        l = array((0.6, 0.6))
56        s = array((0.8, 0.8))
57
58    # number of cmap1 colours is 256 in this case.
59    w.plscmap1n(256)
60    # Interpolate between control points to set up cmap1.
61    w.plscmap1l(0, i, h, l, s)
62# main
63#
64# Does a series of 3-d plots for a given data set, with different
65# viewing options in each plot.
66
67def main(w):
68
69    rosen = 0
70    dx = 2. / float( XPTS - 1 )
71    dy = 2. / float( YPTS - 1 )
72
73    x = -1. + dx*arange(XPTS)
74    y = -1. + dy*arange(YPTS)
75    if rosen == 1:
76        x = 1.5*x
77        y = 0.5 + y
78    x.shape = (-1,1)
79    r2 = (x*x) + (y*y)
80    if rosen == 1:
81        z = (1. - x)*(1. - x) + 100 * (x*x - y)*(x*x - y)
82        # The log argument might be zero for just the right grid.
83        z = log(choose(greater(z,0.), (exp(-5.), z)))
84    else:
85        z = exp(-r2)*cos((2.0*pi)*sqrt(r2))
86
87    x.shape = (-1,)
88    zmin = min(z.flat)
89    zmax = max(z.flat)
90    nlevel = 10
91    step = (zmax-zmin)/(nlevel+1)
92    clevel = zmin + step + arange(nlevel)*step
93
94    # Set up data and arrays for w.plsurf3dl call below.
95    indexxmin = 0
96    indexxmax = XPTS
97    # Must be same shape as z, and a row of z.
98    zlimited = empty(z.shape)
99    indexymin = empty(z.shape[0], dtype=int)
100    indexymax = empty(z.shape[0], dtype=int)
101    # Parameters of ellipse that limits the data.
102    x0 = 0.5*(XPTS - 1)
103    a = 0.9*x0
104    y0 = 0.5*(YPTS - 1)
105    b = 0.7*y0
106    for i in range(indexxmin, indexxmax):
107         square_root = sqrt(1. - min(1., ((double(i) - x0)/a)**2))
108         # Add 0.5 to find nearest integer and therefore preserve symmetry
109         # with regard to lower and upper bound of y range.
110         indexymin[i] = max(0, int(0.5 + y0 - b*square_root))
111         # indexymax calculated with the convention that it is 1
112         # greater than highest valid index.
113         indexymax[i] = min(YPTS, 1 + int(0.5 + y0 + b*square_root))
114         zlimited[i][indexymin[i]:indexymax[i]] = z[i][indexymin[i]:indexymax[i]]
115
116    w.pllightsource(1., 1., 1.)
117
118    for k in range(2):
119        for ifshade in range(5):
120            w.pladv(0)
121            w.plvpor(0.0, 1.0, 0.0, 0.9)
122            w.plwind(-1.0, 1.0, -0.9, 1.1)
123            w.plcol0(3)
124            w.plmtex("t", 1.0, 0.5, 0.5, title[k])
125            w.plcol0(1)
126            if rosen == 1:
127                w.plw3d(1.0, 1.0, 1.0, -1.5, 1.5, -0.5, 1.5, zmin, zmax,
128                alt[k], az[k])
129            else:
130                w.plw3d(1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, zmin, zmax,
131                alt[k], az[k])
132            w.plbox3("bnstu", "x axis", 0.0, 0,
133            "bnstu", "y axis", 0.0, 0,
134            "bcdmnstuv", "z axis", 0.0, 0)
135
136            w.plcol0(2)
137            if ifshade == 0:
138                # diffuse light surface plot.
139                # set up modified gray scale cmap1.
140                cmap1_init(w, 1)
141                w.plsurf3d(x, y, z, 0, ())
142            elif ifshade == 1:
143                # magnitude colored plot.
144                cmap1_init(w, 0)
145                w.plsurf3d(x, y, z, w.MAG_COLOR, ())
146            elif ifshade == 2:
147                # magnitude colored plot with faceted squares
148                cmap1_init(w, 0)
149                w.plsurf3d(x, y, z, w.MAG_COLOR | w.FACETED, ())
150            elif ifshade == 3:
151                # magnitude colored plot with contours
152                cmap1_init(w, 0)
153                w.plsurf3d(x, y, z, w.MAG_COLOR | w.SURF_CONT | w.BASE_CONT, clevel)
154            elif ifshade == 4:
155                # magnitude colored plot with contoursmagnitude colored plot and index limits
156                cmap1_init(w, 0)
157                w.plsurf3dl(x, y, zlimited, w.MAG_COLOR | w.SURF_CONT | w.BASE_CONT, clevel, indexxmin, indexymin, indexymax)
158
159    # Restore defaults
160    # cmap1 default color palette.
161    w.plspal1("cmap1_default.pal",1)
162
163    # Must be done independently because otherwise this changes output files
164    # and destroys agreement with C examples.
165    #w.plcol0(1)
166