1#  Copyright (C) 2004, 2005, 2006, 2007, 2008 Andrew Ross
2#  Copyright (C) 2004-2016 Alan W. Irwin
3
4#  Simple vector plot example.
5#
6#  This file is part of PLplot.
7#
8#  PLplot is free software; you can redistribute it and/or modify
9#  it under the terms of the GNU Library General Public License as published
10#  by the Free Software Foundation; either version 2 of the License, or
11#  (at your option) any later version.
12#
13#  PLplot is distributed in the hope that it will be useful,
14#  but WITHOUT ANY WARRANTY; without even the implied warranty of
15#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16#  GNU Library General Public License for more details.
17#
18#  You should have received a copy of the GNU Library General Public License
19#  along with PLplot; if not, write to the Free Software
20#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21#
22
23from numpy import *
24
25# Pairs of points making the line segments used to plot the user defined arrow
26arrow_x = [-0.5, 0.5, 0.3, 0.5, 0.3, 0.5]
27arrow_y = [0.0, 0.0, 0.2, 0.0, -0.2, 0.0]
28arrow2_x = [-0.5, 0.3, 0.3, 0.5, 0.3, 0.3]
29arrow2_y = [0.0, 0.0,   0.2, 0.0, -0.2, 0.0]
30xmax = 0.0
31
32def circulation(w):
33
34    nx = 20
35    ny = 20
36
37    dx = 1.0
38    dy = 1.0
39
40    xmin = -nx/2*dx
41    xmax = nx/2*dx
42    ymin = -ny/2*dy
43    ymax = ny/2*dy
44
45    # Create data - circulation around the origin.
46    ix = ones(nx)
47    iy = ones(ny)
48    x = arange(nx)+0.5-nx/2.0
49    y = arange(ny)+0.5-ny/2.0
50    xg = multiply.outer(x,iy)
51    yg = multiply.outer(ix,y)
52    u = yg
53    v = -xg
54
55    # Plot vectors with default arrows
56    w.plenv(xmin, xmax, ymin, ymax, 0, 0)
57    w.pllab("(x)", "(y)", "#frPLplot Example 22 - circulation")
58    w.plcol0(2)
59    scaling = 0.0
60    w.plvect(u,v,scaling,w.pltr2,xg,yg)
61    w.plcol0(1)
62
63
64# Vector plot of flow through a constricted pipe
65def constriction( w, astyle ):
66
67    nx = 20
68    ny = 20
69
70    dx = 1.0
71    dy = 1.0
72
73    xmin = -nx/2*dx
74    xmax = nx/2*dx
75    ymin = -ny/2*dy
76    ymax = ny/2*dy
77
78    Q = 2.0
79    ix = ones(nx)
80    iy = ones(ny)
81    x = (arange(nx)-nx/2+0.5)*dx
82    y = (arange(ny)-ny/2+0.5)*dy
83    xg = multiply.outer(x,iy)
84    yg = multiply.outer(ix,y)
85    b = ymax/4.0*(3-cos(pi*x/xmax))
86    b2 = multiply.outer(b,iy)
87    mask = greater.outer(b,abs(y))
88    dbdx = ymax/4.0*(sin(pi*xg/xmax)*pi/xmax*yg/b2)
89    u = Q*ymax/b2*mask
90    v = dbdx*u
91
92    w.plenv(xmin, xmax, ymin, ymax, 0, 0)
93    w.pllab("(x)", "(y)", "#frPLplot Example 22 - constriction (arrow style "+str(astyle)+")")
94    w.plcol0(2)
95    scaling=-1.0
96    w.plvect(u,v,scaling,w.pltr2,xg,yg)
97    w.plcol0(1)
98
99def transform( x, y, xt, yt, data ):
100
101    xt[0] = x
102    yt[0] = y / 4.0 * ( 3 - cos( pi * x / xmax ) )
103
104
105# Vector plot of flow through a constricted pipe
106def constriction2(w):
107
108    global xmax
109
110    nx = 20
111    ny = 20
112    nc = 11
113    nseg = 20
114
115    dx = 1.0
116    dy = 1.0
117
118    xmin = -nx/2*dx
119    xmax = nx/2*dx
120    ymin = -ny/2*dy
121    ymax = ny/2*dy
122
123    w.plstransform( transform, None )
124
125    Q = 2.0
126    ix = ones(nx)
127    iy = ones(ny)
128    x = (arange(nx)-nx/2+0.5)*dx
129    y = (arange(ny)-ny/2+0.5)*dy
130    xg = multiply.outer(x,iy)
131    yg = multiply.outer(ix,y)
132    b = ymax/4.0*(3-cos(pi*x/xmax))
133    b2 = multiply.outer(b,iy)
134    u = Q*ymax/b2
135    v = multiply.outer(zeros(nx),iy)
136
137    clev = Q + arange(nc)*Q/(nc-1)
138
139    w.plenv(xmin, xmax, ymin, ymax, 0, 0)
140    w.pllab("(x)", "(y)", "#frPLplot Example 22 - constriction with plstransform")
141
142    w.plcol0(2)
143    w.plshades(u,xmin+dx/2,xmax-dx/2,ymin+dy/2,ymax-dy/2,clev,0.0,1,1.0,0,None,None)
144    scaling=-1.0
145    w.plvect(u,v,scaling,w.pltr2,xg,yg)
146    w.plpath(nseg,xmin,ymax,xmax,ymax)
147    w.plpath(nseg,xmin,ymin,xmax,ymin)
148    w.plcol0(1)
149
150    w.plstransform(None,None)
151
152# Vector plot of the gradient of a shielded potential (see example 9)
153def potential(w):
154    nper = 100
155    nlevel = 10
156    nr = 20
157    ntheta = 20
158
159    # Create data to be contoured
160    r = 0.5+arange(nr)
161    r.shape = (-1,1)
162    theta = (2.*pi/float(ntheta-1))*(0.5+arange(ntheta))
163    xg = r*cos(theta)
164    yg = r*sin(theta)
165
166    rmax = nr
167    xmin = min(xg.flat)
168    xmax = max(xg.flat)
169    ymin = min(yg.flat)
170    ymax = max(yg.flat)
171
172    x = xg
173    y = yg
174
175    # Potential inside a conducting cylinder (or sphere) by method of images.
176    # Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
177    # Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
178    # Also put in smoothing term at small distances.
179
180    eps = 2.
181
182    q1 = 1.
183    d1 = rmax/4.
184
185    q1i = - q1*rmax/d1
186    d1i = rmax**2/d1
187
188    q2 = -1.
189    d2 = rmax/4.
190
191    q2i = - q2*rmax/d2
192    d2i = rmax**2/d2
193
194    div1 = sqrt((x-d1)**2 + (y-d1)**2 + eps**2)
195    div1i = sqrt((x-d1i)**2 + (y-d1i)**2 + eps**2)
196    div2 = sqrt((x-d2)**2 + (y+d2)**2 + eps**2)
197    div2i = sqrt((x-d2i)**2 + (y+d2i)**2 + eps**2)
198    zg = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
199    zmin = min(zg.flat)
200    zmax = max(zg.flat)
201
202    ug = (-q1*(x-d1)/div1**3 - q1i*(x-d1i)/div1i**3
203          - q2*(x-d2)/div2**3 - q2i*(x-d2i)/div2i**3 )
204    vg = (-q1*(y-d1)/div1**3 - q1i*(y-d1i)/div1i**3
205          - q2*(y+d2)/div2**3 - q2i*(y+d2i)/div2i**3 )
206
207    umin = min(ug.flat)
208    umax = max(ug.flat)
209    vmin = min(vg.flat)
210    vmax = max(vg.flat)
211
212    w.plenv(xmin, xmax, ymin, ymax, 0, 0)
213    w.pllab("(x)", "(y)", "#frPLplot Example 22 - potential gradient vector plot")
214    # Plot contours of the potential
215    dz = (zmax-zmin)/float(nlevel)
216    clevel = zmin + (arange(nlevel)+0.5)*dz
217    du = (umax-umin)/float(nlevel)
218    clevelu = umin + (arange(nlevel)+0.5)*du
219    dv = (vmax-vmin)/float(nlevel)
220    clevelv = vmin + (arange(nlevel)+0.5)*dv
221
222    w.plcol0(3)
223    w.pllsty(2)
224    w.plcont(zg,clevel,w.pltr2,xg,yg)
225    w.pllsty(1)
226    w.plcol0(1)
227
228    # Plot the vectors of the gradient of the potential
229    w.plcol0(2)
230    scaling = 25.0
231    w.plvect(ug,vg,scaling,w.pltr2,xg,yg)
232    w.plcol0(1)
233
234    # Perimeter
235    t = (2.*pi/(nper-1))*arange(nper)
236    px = rmax*cos(t)
237    py = rmax*sin(t)
238    w.plline(px,py)
239
240# main
241#
242# Does a series of vector plots
243#
244def main(w):
245
246    circulation(w)
247
248    narr = 6
249    fill = 0
250
251# Set arrow style using arrow_x and arrow_y then
252# plot using these arrows.
253    w.plsvect(arrow_x, arrow_y, fill)
254    constriction(w, 1)
255
256# Set arrow style using arrow2_x and arrow2_y then
257# plot using these filled arrows.
258    fill = 1
259    w.plsvect(arrow2_x, arrow2_y, fill)
260    constriction(w, 2)
261
262    constriction2(w)
263
264    w.plsvect( None, None, 0)
265
266    potential(w)
267
268    # Restore defaults
269    # Must be done independently because otherwise this changes output files
270    # and destroys agreement with C examples.
271    #w.plcol0(1)
272
273