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