Root/nanonote-files/example-files/data/Examples/lua-plplot-examples/x09.lua

1--[[ $Id: x09.lua 9533 2009-02-16 22:18:37Z smekal $
2
3    Contour plot demo.
4
5  Copyright (C) 2008 Werner Smekal
6
7  This file is part of PLplot.
8
9  PLplot is free software you can redistribute it and/or modify
10  it under the terms of the GNU General Library Public License as published
11  by the Free Software Foundation either version 2 of the License, or
12  (at your option) any later version.
13
14  PLplot is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU Library General Public License for more details.
18
19  You should have received a copy of the GNU Library General Public License
20  along with PLplot if not, write to the Free Software
21  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22--]]
23
24-- initialise Lua bindings for PLplot examples.
25dofile("plplot_examples.lua")
26
27XPTS = 35 -- Data points in x
28YPTS = 46 -- Data points in y
29
30XSPA = 2/(XPTS-1)
31YSPA = 2/(YPTS-1)
32
33-- polar plot data
34PERIMETERPTS = 100
35RPTS = 40
36THETAPTS = 40
37
38-- potential plot data
39PPERIMETERPTS = 100
40PRPTS = 40
41PTHETAPTS = 64
42PNLEVEL = 20
43
44clevel = { -1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1}
45
46-- Transformation function
47tr = { XSPA, 0, -1, 0, YSPA, -1 }
48
49function mypltr(x, y)
50    tx = tr[1] * x + tr[2] * y + tr[3]
51    ty = tr[4] * x + tr[5] * y + tr[6]
52    
53    return tx, ty
54end
55
56--polar contour plot example.
57function polar()
58    px = {}
59    py = {}
60    lev = {}
61
62    pl.env(-1, 1, -1, 1, 0, -2)
63    pl.col0(1)
64         
65    --Perimeter
66    for i=1, PERIMETERPTS do
67        t = (2*math.pi/(PERIMETERPTS-1))*(i-1)
68        px[i] = math.cos(t)
69        py[i] = math.sin(t)
70  end
71  pl.line(px, py)
72           
73    --create data to be contoured.
74  cgrid2["xg"] = {}
75  cgrid2["yg"] = {}
76  cgrid2["nx"] = RPTS
77  cgrid2["ny"] = THETAPTS
78  z = {}
79   
80  for i = 1, RPTS do
81    r = (i-1)/(RPTS-1)
82    cgrid2["xg"][i] = {}
83    cgrid2["yg"][i] = {}
84    z[i] = {}
85    for j = 1, THETAPTS do
86      theta = (2*math.pi/(THETAPTS-1))*(j-1)
87      cgrid2["xg"][i][j] = r*math.cos(theta)
88      cgrid2["yg"][i][j] = r*math.sin(theta)
89      z[i][j] = r
90    end
91  end
92
93  for i = 1, 10 do
94    lev[i] = 0.05 + 0.10*(i-1)
95  end
96
97  pl.col0(2)
98  pl.cont(z, 1, RPTS, 1, THETAPTS, lev, "pltr2", cgrid2)
99  pl.col0(1)
100  pl.lab("", "", "Polar Contour Plot")
101end
102
103
104----------------------------------------------------------------------------
105-- f2mnmx
106--
107-- Returns min & max of input 2d array.
108----------------------------------------------------------------------------
109function f2mnmx(f, nx, ny)
110    fmax = f[1][1]
111    fmin = fmax
112
113    for i=1, nx do
114        for j=1, ny do
115        fmax = math.max(fmax, f[i][j])
116        fmin = math.min(fmin, f[i][j])
117        end
118    end
119        
120    return fmin, fmax
121end
122
123
124--shielded potential contour plot example.
125function potential()
126  clevelneg = {}
127  clevelpos = {}
128  px = {}
129  py = {}
130
131  --create data to be contoured.
132  cgrid2["xg"] = {}
133  cgrid2["yg"] = {}
134  cgrid2["nx"] = PRPTS
135  cgrid2["ny"] = PTHETAPTS
136  z = {}
137
138  for i = 1, PRPTS do
139    r = 0.5 + (i-1)
140    cgrid2["xg"][i] = {}
141    cgrid2["yg"][i] = {}
142    for j = 1, PTHETAPTS do
143      theta = 2*math.pi/(PTHETAPTS-1)*(j-0.5)
144      cgrid2["xg"][i][j] = r*math.cos(theta)
145      cgrid2["yg"][i][j] = r*math.sin(theta)
146    end
147  end
148
149  rmax = PRPTS-0.5
150  xmin, xmax = f2mnmx(cgrid2["xg"], PRPTS, PTHETAPTS)
151  ymin, ymax = f2mnmx(cgrid2["yg"], PRPTS, PTHETAPTS)
152  x0 = (xmin + xmax)/2
153  y0 = (ymin + ymax)/2
154
155  -- Expanded limits
156  peps = 0.05
157  xpmin = xmin - math.abs(xmin)*peps
158  xpmax = xmax + math.abs(xmax)*peps
159  ypmin = ymin - math.abs(ymin)*peps
160  ypmax = ymax + math.abs(ymax)*peps
161
162  -- Potential inside a conducting cylinder (or sphere) by method of images.
163  -- Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
164  -- Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
165  -- Also put in smoothing term at small distances.
166  eps = 2
167  q1 = 1
168  d1 = rmax/4
169
170  q1i = - q1*rmax/d1
171  d1i = rmax^2/d1
172
173  q2 = -1
174  d2 = rmax/4
175
176  q2i = - q2*rmax/d2
177  d2i = rmax^2/d2
178
179  for i = 1, PRPTS do
180    z[i] = {}
181    for j = 1, PTHETAPTS do
182      div1 = math.sqrt((cgrid2.xg[i][j]-d1)^2 + (cgrid2.yg[i][j]-d1)^2 + eps^2)
183      div1i = math.sqrt((cgrid2.xg[i][j]-d1i)^2 + (cgrid2.yg[i][j]-d1i)^2 + eps^2)
184      div2 = math.sqrt((cgrid2.xg[i][j]-d2)^2 + (cgrid2.yg[i][j]+d2)^2 + eps^2)
185      div2i = math.sqrt((cgrid2.xg[i][j]-d2i)^2 + (cgrid2.yg[i][j]+d2i)^2 + eps^2)
186      z[i][j] = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
187    end
188  end
189  zmin, zmax = f2mnmx(z, PRPTS, PTHETAPTS)
190
191  -- Positive and negative contour levels.
192  dz = (zmax-zmin)/PNLEVEL
193  nlevelneg = 1
194  nlevelpos = 1
195  for i = 1, PNLEVEL do
196    clevel = zmin + (i-0.5)*dz
197    if clevel <= 0 then
198      clevelneg[nlevelneg] = clevel
199      nlevelneg = nlevelneg + 1
200    else
201      clevelpos[nlevelpos] = clevel
202      nlevelpos = nlevelpos + 1
203    end
204  end
205  
206  -- Colours!
207  ncollin = 11
208  ncolbox = 1
209  ncollab = 2
210
211  -- Finally start plotting this page!
212  pl.adv(0)
213  pl.col0(ncolbox)
214
215  pl.vpas(0.1, 0.9, 0.1, 0.9, 1)
216  pl.wind(xpmin, xpmax, ypmin, ypmax)
217  pl.box("", 0, 0, "", 0, 0)
218
219  pl.col0(ncollin)
220  if nlevelneg>1 then
221    -- Negative contours
222    pl.lsty(2)
223    pl.cont(z, 1, PRPTS, 1, PTHETAPTS, clevelneg, "pltr2", cgrid2)
224  end
225
226  if nlevelpos>1 then
227    -- Positive contours
228    pl.lsty(1)
229    pl.cont(z, 1, PRPTS, 1, PTHETAPTS, clevelpos, "pltr2", cgrid2)
230  end
231   
232  -- Draw outer boundary
233  for i = 1, PPERIMETERPTS do
234    t = (2*math.pi/(PPERIMETERPTS-1))*(i-1)
235    px[i] = x0 + rmax*math.cos(t)
236    py[i] = y0 + rmax*math.sin(t)
237  end
238
239  pl.col0(ncolbox)
240  pl.line(px, py)
241       
242  pl.col0(ncollab)
243  pl.lab("", "", "Shielded potential of charges in a conducting sphere")
244end
245  
246
247----------------------------------------------------------------------------
248-- main
249--
250-- Does several contour plots using different coordinate mappings.
251----------------------------------------------------------------------------
252mark = { 1500 }
253space = { 1500 }
254
255-- Parse and process command line arguments
256pl.parseopts(arg, pl.PL_PARSE_FULL)
257
258-- Initialize plplot
259pl.init()
260
261-- Set up function arrays
262z = {}
263w = {}
264
265for i = 1, XPTS do
266    xx = (i-1 - math.floor(XPTS/2)) / math.floor(XPTS/2)
267  z[i] = {}
268  w[i] = {}
269    for j = 1, YPTS do
270        yy = (j-1 - math.floor(YPTS/2)) / math.floor(YPTS/2) - 1
271        z[i][j] = xx^2 - yy^2
272        w[i][j] = 2 * xx * yy
273    end
274end
275
276-- Set up grids
277cgrid1 = {}
278cgrid1["xg"] = {}
279cgrid1["yg"] = {}
280cgrid1["nx"] = XPTS
281cgrid1["ny"] = YPTS
282cgrid2 = {}
283cgrid2["xg"] = {}
284cgrid2["yg"] = {}
285cgrid2["nx"] = XPTS
286cgrid2["ny"] = YPTS
287
288for i = 1, XPTS do
289  cgrid2["xg"][i] = {}
290  cgrid2["yg"][i] = {}
291    for j = 1, YPTS do
292    xx, yy = mypltr(i-1, j-1)
293
294    argx = xx * math.pi/2
295    argy = yy * math.pi/2
296    distort = 0.4
297
298    cgrid1["xg"][i] = xx + distort * math.cos(argx)
299    cgrid1["yg"][j] = yy - distort * math.cos(argy)
300
301    cgrid2["xg"][i][j] = xx + distort * math.cos(argx) * math.cos(argy)
302    cgrid2["yg"][i][j] = yy - distort * math.cos(argx) * math.cos(argy)
303  end
304end
305
306-- Plot using identity transform
307pl.setcontlabelformat(4, 3)
308pl.setcontlabelparam(0.006, 0.3, 0.1, 1)
309pl.env(-1, 1, -1, 1, 0, 0)
310pl.col0(2)
311pl.cont(z, 1, XPTS, 1, YPTS, clevel, "mypltr")
312pl.styl(mark, space)
313pl.col0(3)
314pl.cont(w, 1, XPTS, 1, YPTS, clevel, "mypltr")
315pl.styl({}, {})
316pl.col0(1)
317pl.lab("X Coordinate", "Y Coordinate", "Streamlines of flow")
318pl.setcontlabelparam(0.006, 0.3, 0.1, 0)
319    
320-- Plot using 1d coordinate transform
321pl.env(-1, 1, -1, 1, 0, 0)
322pl.col0(2)
323pl.cont(z, 1, XPTS, 1, YPTS, clevel, "pltr1", cgrid1)
324pl.styl(mark, space)
325pl.col0(3)
326pl.cont(w, 1, XPTS, 1, YPTS, clevel, "pltr1", cgrid1)
327pl.styl({}, {})
328pl.col0(1)
329pl.lab("X Coordinate", "Y Coordinate", "Streamlines of flow")
330    
331-- Plot using 2d coordinate transform
332pl.env(-1, 1, -1, 1, 0, 0)
333pl.col0(2)
334pl.cont(z, 1, XPTS, 1, YPTS, clevel, "pltr2", cgrid2)
335
336pl.styl(mark, space)
337pl.col0(3)
338pl.cont(w, 1, XPTS, 1, YPTS, clevel, "pltr2", cgrid2)
339pl.styl({}, {})
340pl.col0(1)
341pl.lab("X Coordinate", "Y Coordinate", "Streamlines of flow")
342
343pl.setcontlabelparam(0.006, 0.3, 0.1, 0)
344polar()
345
346pl.setcontlabelparam(0.006, 0.3, 0.1, 0)
347potential()
348
349-- Clean up
350pl.plend()
351

Archive Download this file

Branches:
backfire
master
release_2010-11-17
release_2010-12-14
release_2011-02-23
release_2011-05-28
release_2011-08-27
release_2011-11-13
release_2012-03-18
release_2012-04-09
release_2012-07-11
release_2012-10-11



interactive