1#  Copyright (C) 2001-2017 Alan W. Irwin
2
3#  Contour plot demo.
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#
21
22from numpy import *
23
24XPTS = 35
25YPTS = 46
26XSPA = 2./(XPTS-1)
27YSPA = 2./(YPTS-1)
28#polar plot data
29PERIMETERPTS = 100
30RPTS = 40
31THETAPTS = 40
32
33#potential plot data
34PPERIMETERPTS = 100
35PRPTS = 40
36PTHETAPTS = 64
37
38tr = array((XSPA, 0.0, -1.0, 0.0, YSPA, -1.0))
39
40def mypltr(x, y, data):
41    result0 = data[0] * x + data[1] * y + data[2]
42    result1 = data[3] * x + data[4] * y + data[5]
43    return array((result0, result1))
44
45def polar(w):
46    #polar contour plot example.
47    w.plenv(-1., 1., -1., 1., 0, -2,)
48    w.plcol0(1)
49
50    # Perimeter
51    t = (2.*pi/(PERIMETERPTS-1))*arange(PERIMETERPTS)
52    px = cos(t)
53    py = sin(t)
54    w.plline(px, py)
55
56    # create data to be contoured.
57    r = arange(RPTS)/float(RPTS-1)
58    r.shape = (-1,1)
59    theta = (2.*pi/float(THETAPTS-1))*arange(THETAPTS-1)
60    xg = r*cos(theta)
61    yg = r*sin(theta)
62    zg = r*ones(THETAPTS-1)
63
64    lev = 0.05 + 0.10*arange(10)
65
66    w.plcol0(2)
67    w.plcont(zg, lev, w.pltr2, xg, yg, 2)
68    #                              ^-- :-).  Means: "2nd coord is wrapped."
69    w.plcol0(1)
70    w.pllab("", "", "Polar Contour Plot")
71
72def potential(w):
73    #shielded potential contour plot example.
74
75    # create data to be contoured.
76    r = 0.5 + arange(PRPTS)
77    r.shape = (-1,1)
78    theta = (2.*pi/float(PTHETAPTS-1))*(0.5+arange(PTHETAPTS-1))
79    xg = r*cos(theta)
80    yg = r*sin(theta)
81
82    rmax = max(r.flat)
83    xmin = min(xg.flat)
84    xmax = max(xg.flat)
85    ymin = min(yg.flat)
86    ymax = max(yg.flat)
87    x0 = (xmin + xmax)/2.
88    y0 = (ymin + ymax)/2.
89    #perturbed (expanded) limits
90
91    peps = 0.05
92    xpmin = xmin - abs(xmin)*peps
93    xpmax = xmax + abs(xmax)*peps
94    ypmin = ymin - abs(ymin)*peps
95    ypmax = ymax + abs(ymax)*peps
96
97
98    # Potential inside a conducting cylinder (or sphere) by method of images.
99    # Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
100    # Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
101    # Also put in smoothing term at small distances.
102
103    eps = 2.
104
105    q1 = 1.
106    d1 = rmax/4.
107
108    q1i = - q1*rmax/d1
109    d1i = rmax**2/d1
110
111    q2 = -1.
112    d2 = rmax/4.
113
114    q2i = - q2*rmax/d2
115    d2i = rmax**2/d2
116
117    div1 = sqrt((xg-d1)**2 + (yg-d1)**2 + eps**2)
118    div1i = sqrt((xg-d1i)**2 + (yg-d1i)**2 + eps**2)
119
120    div2 = sqrt((xg-d2)**2 + (yg+d2)**2 + eps**2)
121    div2i = sqrt((xg-d2i)**2 + (yg+d2i)**2 + eps**2)
122
123    zg = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
124
125    zmin = min(zg.flat)
126    zmax = max(zg.flat)
127#    print "%.15g %.15g %.15g %.15g %.15g %.15g %.15g %.15g \n" % \
128#    (q1, d1, q1i, d1i, q2, d2, q2i, d2i)
129#    print "%.15g %.15g %.15g %.15g %.15g %.15g \n" % \
130#    (xmin, xmax, ymin, ymax, zmin, zmax)
131
132    # Positive and negative contour levels.
133    nlevel = 20
134    dz = (zmax-zmin)/float(nlevel)
135    clevel = zmin + (arange(20)+0.5)*dz
136    clevelpos = compress(clevel > 0., clevel)
137    clevelneg = compress(clevel <= 0., clevel)
138
139    #Colours!
140    ncollin = 11
141    ncolbox = 1
142    ncollab = 2
143
144    #Finally start plotting this page!
145    w.pladv(0)
146    w.plcol0(ncolbox)
147
148    w.plvpas(0.1, 0.9, 0.1, 0.9, 1.0)
149    w.plwind(xpmin, xpmax, ypmin, ypmax)
150    w.plbox("", 0., 0, "", 0., 0)
151
152    w.plcol0(ncollin)
153    # Negative contours
154    w.pllsty(2)
155    w.plcont(zg, clevelneg, w.pltr2, xg, yg, 2)
156
157    # Positive contours
158    w.pllsty(1)
159    w.plcont(zg, clevelpos, w.pltr2, xg, yg, 2)
160
161
162    # Draw outer boundary
163    t = (2.*pi/(PPERIMETERPTS-1))*arange(PPERIMETERPTS)
164    px = x0 + rmax*cos(t)
165    py = y0 + rmax*sin(t)
166
167    w.plcol0(ncolbox)
168    w.plline(px, py)
169
170    w.plcol0(ncollab)
171    w.pllab("", "", "Shielded potential of charges in a conducting sphere")
172
173def main(w):
174
175    mark = 1500
176    space = 1500
177
178    clevel = -1. + 0.2*arange(11)
179
180    xx = (arange(XPTS) - XPTS//2) / float((XPTS//2))
181    yy = (arange(YPTS) - YPTS//2) / float((YPTS//2)) - 1.
182    xx.shape = (-1,1)
183    z = (xx*xx)-(yy*yy)
184    # 2.*outerproduct(xx,yy) for new versions of Numeric which have outerproduct.
185    w_array = 2.*xx*yy
186
187    # Set up grids.
188
189    # Note *for the given* tr, mypltr(i,j,tr)[0] is only a function of i
190    # and mypltr(i,j,tr)[1] is only function of j.
191    xg0 = mypltr(arange(XPTS),0,tr)[0]
192    yg0 = mypltr(0,arange(YPTS),tr)[1]
193
194    distort = 0.4
195    cos_x = cos((pi/2.)*xg0)
196    cos_y = cos((pi/2.)*yg0)
197    xg1 = xg0 + distort*cos_x
198    yg1 = yg0 - distort*cos_y
199    # Need independent copy here so the shape changes for xg0t do not affect
200    # xg0.
201    xg0t = xg0.copy()
202    cos_x.shape = (-1,1)
203    xg0t.shape = (-1,1)
204    xg2 = xg0t + distort*cos_x*cos_y
205    yg2 = yg0 - distort*cos_x*cos_y
206
207    # Plot using mypltr (scaled identity) transformation used to create
208    # xg0 and yg0
209#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
210#    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
211#    w.plcol0(2)
212#    w.plcont(z, clevel, mypltr, tr)
213#    w.plstyl([mark], [space])
214#    w.plcol0(3)
215#    w.plcont(w, clevel, mypltr, tr)
216#    w.plstyl([], [])
217#    w.plcol0(1)
218#    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
219
220    w.pl_setcontlabelformat(4,3)
221    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 1)
222    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
223    w.plcol0(2)
224    w.plcont(z, clevel, mypltr, tr)
225    w.plstyl([mark], [space])
226    w.plcol0(3)
227    w.plcont(w_array, clevel, mypltr, tr)
228    w.plstyl([], [])
229    w.plcol0(1)
230    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
231
232    # Plot using 1D coordinate transformation.
233    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
234    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
235    w.plcol0(2)
236    w.plcont(z, clevel, w.pltr1, xg1, yg1)
237    w.plstyl([mark], [space])
238    w.plcol0(3)
239    w.plcont(w_array, clevel, w.pltr1, xg1, yg1)
240    w.plstyl([], [])
241    w.plcol0(1)
242    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
243
244#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 1)
245#    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
246#    w.plcol0(2)
247#    w.plcont(z, clevel, w.pltr1, xg1, yg1)
248#    w.plstyl([mark], [space])
249#    w.plcol0(3)
250#    w.plcont(w, clevel, w.pltr1, xg1, yg1)
251#    w.plstyl([], [])
252#    w.plcol0(1)
253#    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
254#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
255#
256    # Plot using 2D coordinate transformation.
257    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
258    w.plcol0(2)
259    w.plcont(z, clevel, w.pltr2, xg2, yg2)
260    w.plstyl([mark], [space])
261    w.plcol0(3)
262    w.plcont(w_array, clevel, w.pltr2, xg2, yg2)
263    w.plstyl([], [])
264    w.plcol0(1)
265    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
266
267#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 1)
268#    w.plenv(-1.0, 1.0, -1.0, 1.0, 0, 0)
269#    w.plcol0(2)
270#    w.plcont(z, clevel, w.pltr2, xg2, yg2)
271#    w.plstyl([mark], [space])
272#    w.plcol0(3)
273#    w.plcont(w, clevel, w.pltr2, xg2, yg2)
274#    w.plstyl([], [])
275#    w.plcol0(1)
276#    w.pllab("X Coordinate", "Y Coordinate", "Streamlines of flow")
277#
278#   polar contour examples.
279    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
280    polar(w)
281#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 1)
282#    polar(w)
283
284#   potential contour examples.
285    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
286    potential(w)
287#    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 1)
288#    potential(w)
289
290    # Restore defaults
291    w.pl_setcontlabelparam(0.006, 0.3, 0.1, 0)
292    # Must be done independently because otherwise this changes output files
293    # and destroys agreement with C examples.
294    #w.plcol0(1)
295