1--[[ 2 Simple vector plot example 3 4 Copyright (C) 2008 Werner Smekal 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 23-- initialise Lua bindings for PLplot examples. 24dofile("plplot_examples.lua") 25 26-- Pairs of points making the line segments used to plot the user defined arrow 27arrow_x = { -0.5, 0.5, 0.3, 0.5, 0.3, 0.5 } 28arrow_y = { 0, 0, 0.2, 0, -0.2, 0 } 29arrow2_x = { -0.5, 0.3, 0.3, 0.5, 0.3, 0.3 } 30arrow2_y = { 0, 0, 0.2, 0, -0.2, 0 } 31 32 33-- Vector plot of the circulation about the origin 34function circulation() 35 nx = 20 36 ny = 20 37 dx = 1 38 dy = 1 39 40 xmin = -nx/2*dx 41 xmax = nx/2*dx 42 ymin = -ny/2*dy 43 ymax = ny/2*dy 44 45 cgrid2 = {} 46 cgrid2["xg"] = {} 47 cgrid2["yg"] = {} 48 cgrid2["nx"] = nx 49 cgrid2["ny"] = ny 50 u = {} 51 v = {} 52 53 -- Create data - circulation around the origin. 54 for i = 1, nx do 55 x = (i-1-nx/2+0.5)*dx 56 cgrid2["xg"][i] = {} 57 cgrid2["yg"][i] = {} 58 u[i] = {} 59 v[i] = {} 60 for j=1, ny do 61 y = (j-1-ny/2+0.5)*dy 62 cgrid2["xg"][i][j] = x 63 cgrid2["yg"][i][j] = y 64 u[i][j] = y 65 v[i][j] = -x 66 end 67 end 68 69 -- Plot vectors with default arrows 70 pl.env(xmin, xmax, ymin, ymax, 0, 0) 71 pl.lab("(x)", "(y)", "#frPLplot Example 22 - circulation") 72 pl.col0(2) 73 pl.vect(u, v, 0, "pltr2", cgrid2 ) 74 pl.col0(1) 75end 76 77 78-- Vector plot of flow through a constricted pipe 79function constriction( astyle ) 80 nx = 20 81 ny = 20 82 dx = 1 83 dy = 1 84 85 xmin = -nx/2*dx 86 xmax = nx/2*dx 87 ymin = -ny/2*dy 88 ymax = ny/2*dy 89 90 cgrid2 = {} 91 cgrid2["xg"] = {} 92 cgrid2["yg"] = {} 93 cgrid2["nx"] = nx 94 cgrid2["ny"] = ny 95 u = {} 96 v = {} 97 98 Q = 2 99 for i = 1, nx do 100 x = (i-1-nx/2+0.5)*dx 101 cgrid2["xg"][i] = {} 102 cgrid2["yg"][i] = {} 103 u[i] = {} 104 v[i] = {} 105 for j = 1, ny do 106 y = (j-1-ny/2+0.5)*dy 107 cgrid2["xg"][i][j] = x 108 cgrid2["yg"][i][j] = y 109 b = ymax/4*(3-math.cos(math.pi*x/xmax)) 110 if math.abs(y)<b then 111 dbdx = ymax/4*math.sin(math.pi*x/xmax)*math.pi/xmax*y/b 112 u[i][j] = Q*ymax/b 113 v[i][j] = dbdx*u[i][j] 114 else 115 u[i][j] = 0 116 v[i][j] = 0 117 end 118 end 119 end 120 121 pl.env(xmin, xmax, ymin, ymax, 0, 0) 122 title = string.format( "#frPLplot Example 22 - constriction (arrow style %d)", astyle ) 123 pl.lab("(x)", "(y)", title) 124 pl.col0(2) 125 pl.vect(u, v, -1.0, "pltr2", cgrid2) 126 pl.col0(1) 127end 128 129-- Note this function uses the global variable xmax rather than passing 130-- data as in C. 131function transform(x,y) 132 133 xt = x 134 yt = y / 4.0 * ( 3 - math.cos( math.pi * x / xmax )) 135 136 return xt, yt 137end 138 139-- Vector plot of flow through a constricted pipe 140function constriction2() 141 nx = 20 142 ny = 20 143 nc = 11 144 nseg = 20 145 146 dx = 1 147 dy = 1 148 149 xmin = -nx/2*dx 150 xmax = nx/2*dx 151 ymin = -ny/2*dy 152 ymax = ny/2*dy 153 154 pl.stransform( "transform" ) 155 156 cgrid2 = {} 157 cgrid2["xg"] = {} 158 cgrid2["yg"] = {} 159 cgrid2["nx"] = nx 160 cgrid2["ny"] = ny 161 u = {} 162 v = {} 163 164 Q = 2 165 for i = 1, nx do 166 x = (i-1-nx/2+0.5)*dx 167 cgrid2["xg"][i] = {} 168 cgrid2["yg"][i] = {} 169 u[i] = {} 170 v[i] = {} 171 for j = 1, ny do 172 y = (j-1-ny/2+0.5)*dy 173 cgrid2["xg"][i][j] = x 174 cgrid2["yg"][i][j] = y 175 b = ymax/4*(3-math.cos(math.pi*x/xmax)) 176 u[i][j] = Q*ymax/b 177 v[i][j] = 0.0 178 end 179 end 180 181 clev = {} 182 for i = 1, nc do 183 clev[i] = Q + (i-1)*Q/(nc-1) 184 end 185 186 pl.env(xmin, xmax, ymin, ymax, 0, 0) 187 pl.lab("(x)", "(y)", "#frPLplot Example 22 - constriction with plstransform") 188 pl.col0(2) 189 pl.shades(u, xmin+dx/2, xmax-dx/2, ymin+dy/2, ymax-dy/2, clev, 0.0, 1, 1.0, 0 ); 190 pl.vect(u, v, -1.0, "pltr2", cgrid2) 191 pl.path( nseg, xmin, ymax, xmax, ymax ); 192 pl.path( nseg, xmin, ymin, xmax, ymin ); 193 pl.col0(1) 194 195 pl.stransform() 196 197end 198 199 200function f2mnmx(f, nx, ny) 201 fmax = f[1][1] 202 fmin = fmax 203 204 for i=1, nx do 205 for j=1, ny do 206 fmax = math.max(fmax, f[i][j]) 207 fmin = math.min(fmin, f[i][j]) 208 end 209 end 210 211 return fmin, fmax 212end 213 214-- Vector plot of the gradient of a shielded potential (see example 9) 215function potential() 216 nper = 100 217 nlevel = 10 218 nr = 20 219 ntheta = 20 220 221 u = {} 222 v = {} 223 z = {} 224 clevel = {} 225 px = {} 226 py = {} 227 228 cgrid2 = {} 229 cgrid2["xg"] = {} 230 cgrid2["yg"] = {} 231 cgrid2["nx"] = nr 232 cgrid2["ny"] = ntheta 233 234 -- Potential inside a conducting cylinder (or sphere) by method of images. 235 -- Charge 1 is placed at (d1, d1), with image charge at (d2, d2). 236 -- Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2). 237 -- Also put in smoothing term at small distances. 238 rmax = nr 239 240 eps = 2 241 242 q1 = 1 243 d1 = rmax/4 244 245 q1i = -q1*rmax/d1 246 d1i = rmax^2/d1 247 248 q2 = -1 249 d2 = rmax/4 250 251 q2i = -q2*rmax/d2 252 d2i = rmax^2/d2 253 254 for i = 1, nr do 255 r = i - 0.5 256 cgrid2["xg"][i] = {} 257 cgrid2["yg"][i] = {} 258 u[i] = {} 259 v[i] = {} 260 z[i] = {} 261 for j = 1, ntheta do 262 theta = 2*math.pi/(ntheta-1)*(j-0.5) 263 x = r*math.cos(theta) 264 y = r*math.sin(theta) 265 cgrid2["xg"][i][j] = x 266 cgrid2["yg"][i][j] = y 267 div1 = math.sqrt((x-d1)^2 + (y-d1)^2 + eps^2) 268 div1i = math.sqrt((x-d1i)^2 + (y-d1i)^2 + eps^2) 269 div2 = math.sqrt((x-d2)^2 + (y+d2)^2 + eps^2) 270 div2i = math.sqrt((x-d2i)^2 + (y+d2i)^2 + eps^2) 271 z[i][j] = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i 272 u[i][j] = -q1*(x-d1)/div1^3 - q1i*(x-d1i)/div1i^3 273 -q2*(x-d2)/div2^3 - q2i*(x-d2i)/div2i^3 274 v[i][j] = -q1*(y-d1)/div1^3 - q1i*(y-d1i)/div1i^3 275 -q2*(y+d2)/div2^3 - q2i*(y+d2i)/div2i^3 276 end 277 end 278 279 xmin, xmax = f2mnmx(cgrid2["xg"], nr, ntheta) 280 ymin, ymax = f2mnmx(cgrid2["yg"], nr, ntheta) 281 zmin, zmax = f2mnmx(z, nr, ntheta) 282 283 pl.env(xmin, xmax, ymin, ymax, 0, 0) 284 pl.lab("(x)", "(y)", "#frPLplot Example 22 - potential gradient vector plot") 285 286 -- Plot contours of the potential 287 dz = (zmax-zmin)/nlevel 288 for i = 1, nlevel do 289 clevel[i] = zmin + (i-0.5)*dz 290 end 291 292 pl.col0(3) 293 pl.lsty(2) 294 pl.cont(z, 1, nr, 1, ntheta, clevel, "pltr2", cgrid2) 295 pl.lsty(1) 296 pl.col0(1) 297 298 -- Plot the vectors of the gradient of the potential 299 pl.col0(2) 300 pl.vect(u, v, 25, "pltr2", cgrid2) 301 pl.col0(1) 302 303 -- Plot the perimeter of the cylinder 304 for i=1, nper do 305 theta = 2*math.pi/(nper-1)*(i-1) 306 px[i] = rmax*math.cos(theta) 307 py[i] = rmax*math.sin(theta) 308 end 309 310 pl.line(px, py) 311end 312 313 314---------------------------------------------------------------------------- 315-- main 316-- 317-- Generates several simple vector plots. 318---------------------------------------------------------------------------- 319 320-- Parse and process command line arguments 321pl.parseopts(arg, pl.PL_PARSE_FULL) 322 323-- Initialize plplot 324pl.init() 325 326circulation() 327 328fill = 0 329 330-- Set arrow style using arrow_x and arrow_y then 331-- plot using these arrows. 332pl.svect(arrow_x, arrow_y, fill) 333constriction(1) 334 335-- Set arrow style using arrow2_x and arrow2_y then 336-- plot using these filled arrows. 337fill = 1 338pl.svect(arrow2_x, arrow2_y, fill) 339constriction(2) 340 341constriction2() 342 343pl.svect(nil, nil, 0) 344 345potential() 346 347pl.plend() 348