1 //--------------------------------------------------------------------------
2 // Simple vector plot example
3 //--------------------------------------------------------------------------
4 //
5 //--------------------------------------------------------------------------
6 // Copyright (C) 2004 Andrew Ross
7 //
8 // This file is part of PLplot.
9 //
10 // PLplot is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Library General Public License as published by
12 // the Free Software Foundation; version 2 of the License.
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 //--------------------------------------------------------------------------
25 // Implementation of PLplot example 22 in C++.
26 //--------------------------------------------------------------------------
27
28 #include "plc++demos.h"
29
30 #ifdef PL_USE_NAMESPACE
31 using namespace std;
32 #endif
33
34 //
35 // Global transform function for a constriction using data passed in
36 // This is the same transformation used in constriction.
37 //
38 void
transform(PLFLT x,PLFLT y,PLFLT * xt,PLFLT * yt,PLPointer data)39 transform( PLFLT x, PLFLT y, PLFLT *xt, PLFLT *yt, PLPointer data )
40 {
41 PLFLT *trdata;
42 PLFLT xmax;
43
44 trdata = (PLFLT *) data;
45 xmax = *trdata;
46
47 *xt = x;
48 *yt = y / 4.0 * ( 3 - cos( M_PI * x / xmax ) );
49 }
50
51
52 class x22 {
53 public:
54 x22( int, char ** );
55
56 private:
57 void circulation();
58 void constriction( int astyle );
59 void constriction2();
60 void potential();
61 void f2mnmx( PLFLT **f, PLINT nx, PLINT ny, PLFLT *fmin, PLFLT *fmax );
62
MIN(PLFLT x,PLFLT y)63 PLFLT MIN( PLFLT x, PLFLT y ) { return ( x < y ? x : y ); };
MAX(PLFLT x,PLFLT y)64 PLFLT MAX( PLFLT x, PLFLT y ) { return ( x > y ? x : y ); };
65
66 plstream *pls;
67
68 PLFLT **u, **v;
69 PLcGrid2 cgrid2;
70 int nx, ny, nc, nseg;
71 };
72
73 // Vector plot of the circulation about the origin
74 void
circulation()75 x22::circulation()
76 {
77 int i, j;
78 PLFLT dx, dy, x, y;
79 PLFLT xmin, xmax, ymin, ymax;
80
81 dx = 1.0;
82 dy = 1.0;
83
84 xmin = -nx / 2 * dx;
85 xmax = nx / 2 * dx;
86 ymin = -ny / 2 * dy;
87 ymax = ny / 2 * dy;
88
89
90 // Create data - cirulation around the origin.
91 for ( i = 0; i < nx; i++ )
92 {
93 for ( j = 0; j < ny; j++ )
94 {
95 x = ( i - nx / 2 + 0.5 ) * dx;
96 y = ( j - ny / 2 + 0.5 ) * dy;
97 cgrid2.xg[i][j] = x;
98 cgrid2.yg[i][j] = y;
99 u[i][j] = y;
100 v[i][j] = -x;
101 }
102 }
103
104 // Plot vectors with default arrows
105 pls->env( xmin, xmax, ymin, ymax, 0, 0 );
106 pls->lab( "(x)", "(y)", "#frPLplot Example 22 - circulation" );
107 pls->col0( 2 );
108 pls->vect( u, v, nx, ny, 0.0, plcallback::tr2, (void *) &cgrid2 );
109 pls->col0( 1 );
110 }
111
112 // Vector plot of flow through a constricted pipe
113 void
constriction(int astyle)114 x22::constriction( int astyle )
115 {
116 int i, j;
117 PLFLT dx, dy, x, y;
118 PLFLT xmin, xmax, ymin, ymax;
119 PLFLT Q, b, dbdx;
120 char title[80];
121
122 dx = 1.0;
123 dy = 1.0;
124
125 xmin = -nx / 2 * dx;
126 xmax = nx / 2 * dx;
127 ymin = -ny / 2 * dy;
128 ymax = ny / 2 * dy;
129
130 Q = 2.0;
131 for ( i = 0; i < nx; i++ )
132 {
133 x = ( i - nx / 2 + 0.5 ) * dx;
134 for ( j = 0; j < ny; j++ )
135 {
136 y = ( j - ny / 2 + 0.5 ) * dy;
137 cgrid2.xg[i][j] = x;
138 cgrid2.yg[i][j] = y;
139 b = ymax / 4.0 * ( 3.0 - cos( M_PI * x / xmax ) );
140 if ( fabs( y ) < b )
141 {
142 dbdx = ymax / 4.0 * sin( M_PI * x / xmax ) *
143 M_PI / xmax * y / b;
144 u[i][j] = Q * ymax / b;
145 v[i][j] = dbdx * u[i][j];
146 }
147 else
148 {
149 u[i][j] = 0.0;
150 v[i][j] = 0.0;
151 }
152 }
153 }
154
155 pls->env( xmin, xmax, ymin, ymax, 0, 0 );
156 sprintf( title, "#frPLplot Example 22 - constriction (arrow style %d)", astyle );
157 pls->lab( "(x)", "(y)", title );
158 pls->col0( 2 );
159 pls->vect( u, v, nx, ny, -1.0, plcallback::tr2, (void *) &cgrid2 );
160 pls->col0( 1 );
161 }
162
163 //
164 // Vector plot of flow through a constricted pipe
165 // with a coordinate transform
166 //
167 void
constriction2(void)168 x22::constriction2( void )
169 {
170 int i, j;
171 PLFLT dx, dy, x, y;
172 PLFLT xmin, xmax, ymin, ymax;
173 PLFLT Q, b;
174 #define NC 11
175 int nc = NC;
176 PLFLT clev[NC];
177
178 dx = 1.0;
179 dy = 1.0;
180
181 xmin = -nx / 2 * dx;
182 xmax = nx / 2 * dx;
183 ymin = -ny / 2 * dy;
184 ymax = ny / 2 * dy;
185
186 pls->stransform( transform, ( PLPointer ) & xmax );
187
188 Q = 2.0;
189 for ( i = 0; i < nx; i++ )
190 {
191 x = ( i - nx / 2 + 0.5 ) * dx;
192 for ( j = 0; j < ny; j++ )
193 {
194 y = ( j - ny / 2 + 0.5 ) * dy;
195 cgrid2.xg[i][j] = x;
196 cgrid2.yg[i][j] = y;
197 b = ymax / 4.0 * ( 3 - cos( M_PI * x / xmax ) );
198 u[i][j] = Q * ymax / b;
199 v[i][j] = 0.0;
200 }
201 }
202
203 for ( i = 0; i < nc; i++ )
204 {
205 clev[i] = Q + i * Q / ( nc - 1 );
206 }
207
208 pls->env( xmin, xmax, ymin, ymax, 0, 0 );
209 pls->lab( "(x)", "(y)", "#frPLplot Example 22 - constriction with plstransform" );
210 pls->col0( 2 );
211 pls->shades( (const PLFLT * const *) u, nx, ny, NULL,
212 xmin + dx / 2, xmax - dx / 2, ymin + dy / 2, ymax - dy / 2,
213 clev, nc, 0, 1, 1.0, plcallback::fill, 0, NULL, NULL );
214 pls->vect( (const PLFLT * const *) u, (const PLFLT * const *) v, nx, ny,
215 -1.0, plcallback::tr2, (void *) &cgrid2 );
216 // Plot edges using plpath (which accounts for coordinate transformation) rather than plline
217 pls->path( nseg, xmin, ymax, xmax, ymax );
218 pls->path( nseg, xmin, ymin, xmax, ymin );
219 pls->col0( 1 );
220
221 pls->stransform( NULL, NULL );
222 }
223
224 // Vector plot of the gradient of a shielded potential (see example 9)
225 void
potential()226 x22::potential()
227 {
228 const int nper = 100;
229 const int nlevel = 10;
230
231 int i, j, nr, ntheta;
232 PLFLT eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i;
233 PLFLT div1, div1i, div2, div2i;
234 PLFLT **z, r, theta, x, y, dz;
235 PLFLT xmin, xmax, ymin, ymax, rmax, zmax, zmin;
236 PLFLT px[nper], py[nper], clevel[nlevel];
237
238 nr = nx;
239 ntheta = ny;
240
241 // Create data to be plotted
242 pls->Alloc2dGrid( &z, nr, ntheta );
243
244 // Potential inside a conducting cylinder (or sphere) by method of images.
245 // Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
246 // Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
247 // Also put in smoothing term at small distances.
248
249 rmax = (PLFLT) nr;
250
251 eps = 2.;
252
253 q1 = 1.;
254 d1 = rmax / 4.;
255
256 q1i = -q1 * rmax / d1;
257 d1i = pow( rmax, 2. ) / d1;
258
259 q2 = -1.;
260 d2 = rmax / 4.;
261
262 q2i = -q2 * rmax / d2;
263 d2i = pow( rmax, 2. ) / d2;
264
265 for ( i = 0; i < nr; i++ )
266 {
267 r = 0.5 + (PLFLT) i;
268 for ( j = 0; j < ntheta; j++ )
269 {
270 theta = 2. * M_PI / ( ntheta - 1 ) * ( 0.5 + (PLFLT) j );
271 x = r * cos( theta );
272 y = r * sin( theta );
273 cgrid2.xg[i][j] = x;
274 cgrid2.yg[i][j] = y;
275 div1 = sqrt( pow( ( x - d1 ), 2. ) + pow( ( y - d1 ), 2. ) + pow( eps, 2. ) );
276 div1i = sqrt( pow( ( x - d1i ), 2. ) + pow( ( y - d1i ), 2. ) + pow( eps, 2. ) );
277 div2 = sqrt( pow( ( x - d2 ), 2. ) + pow( ( y + d2 ), 2. ) + pow( eps, 2. ) );
278 div2i = sqrt( pow( ( x - d2i ), 2. ) + pow( ( y + d2i ), 2. ) + pow( eps, 2. ) );
279 z[i][j] = q1 / div1 + q1i / div1i + q2 / div2 + q2i / div2i;
280 u[i][j] = -q1 * ( x - d1 ) / pow( div1, 3. ) - q1i * ( x - d1i ) / pow( div1i, 3.0 )
281 - q2 * ( x - d2 ) / pow( div2, 3. ) - q2i * ( x - d2i ) / pow( div2i, 3. );
282 v[i][j] = -q1 * ( y - d1 ) / pow( div1, 3. ) - q1i * ( y - d1i ) / pow( div1i, 3.0 )
283 - q2 * ( y + d2 ) / pow( div2, 3. ) - q2i * ( y + d2i ) / pow( div2i, 3. );
284 }
285 }
286
287 f2mnmx( cgrid2.xg, nr, ntheta, &xmin, &xmax );
288 f2mnmx( cgrid2.yg, nr, ntheta, &ymin, &ymax );
289 f2mnmx( z, nr, ntheta, &zmin, &zmax );
290
291 pls->env( xmin, xmax, ymin, ymax, 0, 0 );
292 pls->lab( "(x)", "(y)", "#frPLplot Example 22 - potential gradient vector plot" );
293 // Plot contours of the potential
294 dz = ( zmax - zmin ) / (PLFLT) nlevel;
295 for ( i = 0; i < nlevel; i++ )
296 {
297 clevel[i] = zmin + ( (PLFLT) i + 0.5 ) * dz;
298 }
299 pls->col0( 3 );
300 pls->lsty( 2 );
301 pls->cont( z, nr, ntheta, 1, nr, 1, ntheta, clevel, nlevel, plcallback::tr2, (void *) &cgrid2 );
302 pls->lsty( 1 );
303 pls->col0( 1 );
304
305 // Plot the vectors of the gradient of the potential
306 pls->col0( 2 );
307 pls->vect( u, v, nr, ntheta, 25.0, plcallback::tr2, (void *) &cgrid2 );
308 pls->col0( 1 );
309
310 // Plot the perimeter of the cylinder
311 for ( i = 0; i < nper; i++ )
312 {
313 theta = ( 2. * M_PI / ( nper - 1 ) ) * (PLFLT) i;
314 px[i] = rmax * cos( theta );
315 py[i] = rmax * sin( theta );
316 }
317 pls->line( nper, px, py );
318
319 pls->Free2dGrid( z, nr, ntheta );
320 }
321
322 void
f2mnmx(PLFLT ** f,PLINT nx,PLINT ny,PLFLT * fmin,PLFLT * fmax)323 x22::f2mnmx( PLFLT **f, PLINT nx, PLINT ny, PLFLT *fmin, PLFLT *fmax )
324 {
325 int i, j;
326
327 *fmax = f[0][0];
328 *fmin = *fmax;
329
330 for ( i = 0; i < nx; i++ )
331 {
332 for ( j = 0; j < ny; j++ )
333 {
334 *fmax = MAX( *fmax, f[i][j] );
335 *fmin = MIN( *fmin, f[i][j] );
336 }
337 }
338 }
339
340
x22(int argc,char ** argv)341 x22::x22( int argc, char ** argv )
342 {
343 PLINT narr;
344 bool fill;
345
346 // Set of points making a polygon to use as the arrow
347 PLFLT arrow_x[6] = { -0.5, 0.5, 0.3, 0.5, 0.3, 0.5 };
348 PLFLT arrow_y[6] = { 0.0, 0.0, 0.2, 0.0, -0.2, 0.0 };
349 PLFLT arrow2_x[6] = { -0.5, 0.3, 0.3, 0.5, 0.3, 0.3 };
350 PLFLT arrow2_y[6] = { 0.0, 0.0, 0.2, 0.0, -0.2, 0.0 };
351
352 // Create new plstream
353 pls = new plstream();
354
355 // Parse and process command line arguments
356
357 pls->parseopts( &argc, argv, PL_PARSE_FULL );
358
359 // Initialize plplot
360
361 pls->init();
362
363 nx = 20;
364 ny = 20;
365 nc = 11;
366 nseg = 20;
367
368 // Allocate arrays
369 pls->Alloc2dGrid( &cgrid2.xg, nx, ny );
370 pls->Alloc2dGrid( &cgrid2.yg, nx, ny );
371 pls->Alloc2dGrid( &u, nx, ny );
372 pls->Alloc2dGrid( &v, nx, ny );
373
374 cgrid2.nx = nx;
375 cgrid2.ny = ny;
376
377 circulation();
378
379 narr = 6;
380 fill = false;
381
382 // Set arrow style using arrow_x and arrow_y then
383 // plot using these arrows.
384 pls->svect( arrow_x, arrow_y, narr, fill );
385 constriction( 1 );
386
387 // Set arrow style using arrow2_x and arrow2_y then
388 // plot using these filled arrows.
389 fill = true;
390 pls->svect( arrow2_x, arrow2_y, narr, fill );
391 constriction( 2 );
392
393 constriction2();
394
395 // Reset arrow style to the default by passing two
396 // NULL arrays (this are the default arguments)
397 pls->svect( );
398
399 potential();
400
401 pls->Free2dGrid( cgrid2.xg, nx, ny );
402 pls->Free2dGrid( cgrid2.yg, nx, ny );
403 pls->Free2dGrid( u, nx, ny );
404 pls->Free2dGrid( v, nx, ny );
405
406 delete pls;
407 }
408
main(int argc,char ** argv)409 int main( int argc, char ** argv )
410 {
411 x22 *x = new x22( argc, argv );
412 delete x;
413 }
414
415
416 //--------------------------------------------------------------------------
417 // End of x22.cc
418 //--------------------------------------------------------------------------
419
420