1-- Simple vector plot example 2 3-- Copyright (C) 2008, 2013 Jerry Bauck 4 5-- This file is part of PLplot. 6 7-- PLplot is free software; you can redistribute it and/or modify 8-- it under the terms of the GNU Library General Public License as published 9-- by the Free Software Foundation; either version 2 of the License, or 10-- (at your option) any later version. 11 12-- PLplot is distributed in the hope that it will be useful, 13-- but WITHOUT ANY WARRANTY; without even the implied warranty of 14-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15-- GNU Library General Public License for more details. 16 17-- You should have received a copy of the GNU Library General Public License 18-- along with PLplot; if not, write to the Free Software 19-- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 20 21with 22 System, 23 System.Address_To_Access_Conversions, 24 Ada.Numerics, 25 Ada.Numerics.Long_Elementary_Functions, 26 PLplot_Auxiliary, 27 PLplot_Standard; 28use 29 Ada.Numerics, 30 Ada.Numerics.Long_Elementary_Functions, 31 System, 32 PLplot_Auxiliary, 33 PLplot_Standard; 34 35procedure xstandard22a is 36 -- Pairs of points making the line segments used to plot the user defined arrow 37 arrow_x : Real_Vector(0 .. 5) := (-0.5, 0.5, 0.3, 0.5, 0.3, 0.5); 38 arrow_y : Real_Vector(0 .. 5) := ( 0.0, 0.0, 0.2, 0.0, -0.2, 0.0); 39 arrow2_x : Real_Vector(0 .. 5) := (-0.5, 0.3, 0.3, 0.5, 0.3, 0.3); 40 arrow2_y : Real_Vector(0 .. 5) := ( 0.0, 0.0, 0.2, 0.0, -0.2, 0.0); 41 42 43 -- Vector plot of the circulation about the origin 44 procedure circulation is 45 dx, dy, x, y : Long_Float; 46 nx : constant Integer := 20; 47 ny : constant Integer := 20; 48 xmin, xmax, ymin, ymax : Long_Float; 49 u, v : Real_Matrix(0 .. nx - 1, 0 .. ny -1); 50 cgrid2 : aliased Transformation_Data_Type_2 51 (x_Last => nx - 1, 52 y_Last => ny - 1); 53 begin 54 dx := 1.0; 55 dy := 1.0; 56 57 xmin := Long_Float(-nx / 2) * dx; 58 xmax := Long_Float( nx / 2) * dx; 59 ymin := Long_Float(-ny / 2) * dy; 60 ymax := Long_Float( ny / 2) * dy; 61 62 -- Create data - circulation around the origin. 63 for i in 0 .. nx - 1 loop 64 x := (Long_Float(i - nx / 2) + 0.5) * dx; 65 for j in 0 .. ny - 1 loop 66 y := (Long_Float(j - ny / 2) + 0.5) * dy; 67 cgrid2.xg(i, j) := x; 68 cgrid2.yg(i, j) := y; 69 u(i, j) := y; 70 v(i, j) := -x; 71 end loop; 72 end loop; 73 74 -- Plot vectors with default arrows 75 Set_Environment(xmin, xmax, ymin, ymax, Not_Justified, Linear_Box_Plus); 76 Write_Labels("(x)", "(y)", "#frPLplot Example 22 - circulation"); 77 Set_Pen_Color(Yellow); 78 Vector_Plot(u, v, 0.0, Plot_Transformation_2'access, cgrid2'Address); 79 Set_Pen_Color(Red); 80 end circulation; 81 82 83 --Vector plot of flow through a constricted pipe 84 procedure constriction(astyle : Integer) is 85 dx, dy, x, y : Long_Float; 86 xmin, xmax, ymin, ymax : Long_Float; 87 Q, b, dbdx : Long_Float; 88 nx : constant Integer := 20; 89 ny : constant Integer := 20; 90 u, v : Real_Matrix(0 .. nx - 1, 0 .. ny -1); 91 cgrid2 : aliased Transformation_Data_Type_2 92 (x_Last => nx - 1, 93 y_Last => ny - 1); 94 begin 95 dx := 1.0; 96 dy := 1.0; 97 98 xmin := Long_Float(-nx / 2) * dx; 99 xmax := Long_Float( nx / 2) * dx; 100 ymin := Long_Float(-ny / 2) * dy; 101 ymax := Long_Float( ny / 2) * dy; 102 103 Q := 2.0; 104 for i in 0 .. nx - 1 loop 105 x := (Long_Float(i - nx / 2) + 0.5) * dx; 106 for j in 0 .. ny - 1 loop 107 y := (Long_Float(j - ny / 2) + 0.5) * dy; 108 cgrid2.xg(i, j) := x; 109 cgrid2.yg(i, j) := y; 110 b := ymax / 4.0 * (3.0 - cos(pi * x / xmax)); 111 if abs(y) < b then 112 dbdx := ymax / 4.0 * sin(pi * x / xmax) * pi / xmax * y / b; 113 u(i, j) := Q * ymax / b; 114 v(i, j) := dbdx * u(i, j); 115 else 116 u(i, j) := 0.0; 117 v(i, j) := 0.0; 118 end if; 119 end loop; 120 end loop; 121 122 Set_Environment(xmin, xmax, ymin, ymax, Not_Justified, Linear_Box_Plus); 123 Write_Labels("(x)", "(y)", "#frPLplot Example 22 - constriction (arrow style"&Integer'image(astyle)&")"); 124 Set_Pen_Color(Yellow); 125 Vector_Plot(u, v, -1.0, Plot_Transformation_2'access, cgrid2'Address); 126 Set_Pen_Color(Red); 127 end constriction; 128 129 130 -- This spec is necessary in order to enforce C calling conventions, used 131 -- in the callback by intervening C code. 132 procedure transform 133 (x, y : Long_Float; 134 xt, yt : out Long_Float; 135 data : PL_Pointer); 136 pragma Convention(C, transform); 137 138 -- Global transform function for a constriction using data passed in 139 -- This is the same transformation used in constriction. 140 procedure transform(x, y : Long_Float; xt, yt : out Long_Float; Data : PL_Pointer) is 141 142 -- Convert the generic pointer represented as System.Address to a proper Ada pointer aka 143 -- access variable. Recall that PL_Pointer is a subtype of System.Address. 144 package Data_Address_Conversions is new System.Address_To_Access_Conversions(Long_Float); 145 Data_Pointer : Data_Address_Conversions.Object_Pointer; -- An Ada access variable 146 xmax : Long_Float; 147 begin 148 Data_Pointer := Data_Address_Conversions.To_Pointer(Data); 149 xmax := Data_Pointer.all; 150 xt := x; 151 yt := y / 4.0 * (3.0 - cos(Pi * x / xmax)); 152 end transform; 153 154 155 -- Vector plot of flow through a constricted pipe with a coordinate transform 156 procedure constriction2 is 157 dx, dy, x, y : Long_Float; 158 xmin, xmax, ymin, ymax : Long_Float; 159 Q, b : Long_Float; 160 nx : constant Integer := 20; 161 ny : constant Integer := 20; 162 cgrid2 : aliased Transformation_Data_Type_2 163 (x_Last => nx - 1, 164 y_Last => ny - 1); 165 u, v : Real_Matrix(0 .. nx - 1, 0 .. ny - 1); 166 nc : constant Integer := 11; 167 nseg : constant Integer := 20; 168 clev : Real_Vector(0 .. nc - 1); 169 begin 170 dx := 1.0; 171 dy := 1.0; 172 173 xmin := Long_Float(-nx / 2) * dx; -- Careful; Ada / rounds, C / truncates. 174 xmax := Long_Float( nx / 2) * dx; 175 ymin := Long_Float(-ny / 2) * dy; 176 ymax := Long_Float( ny / 2) * dy; 177 178 Set_Custom_Coordinate_Transform(transform'Unrestricted_Access, xmax'Address); 179 180 cgrid2.nx := nx; 181 cgrid2.ny := ny; 182 Q := 2.0; 183 184 for i in 0 .. nx - 1 loop 185 x := (Long_Float(i - nx / 2) + 0.5) * dx; 186 for j in 0 .. ny - 1 loop 187 y := (Long_Float(j - ny / 2) + 0.5) * dy; 188 cgrid2.xg(i, j) := x; 189 cgrid2.yg(i, j) := y; 190 b := ymax / 4.0 * (3.0 - cos(Pi * x / xmax)); 191 u(i, j) := Q * ymax / b; 192 v(i, j) := 0.0; 193 end loop; 194 end loop; 195 196 for i in 0 .. nc - 1 loop 197 clev(i) := Q + Long_Float(i) * Q / Long_Float(nc - 1); 198 end loop; 199 200 Set_Environment(xmin, xmax, ymin, ymax, 0, 0); 201 Write_Labels("(x)", "(y)", "#frPLplot Example 22 - constriction with plstransform"); 202 Set_Pen_Color(Yellow); 203 Shade_Regions(u, Null, 204 xmin + dx / 2.0, xmax - dx / 2.0, ymin + dy / 2.0, ymax - dy / 2.0, 205 clev, 0.0, 1, 1.0, Fill_Polygon'access, False, Null, System.Null_Address); 206 Vector_Plot(u, v, 207 -1.0, Plot_Transformation_2'access, cgrid2'Address); 208 209 -- Plot edges using plpath (which accounts for coordinate transformation) rather than plline 210 Draw_Line_With_Transform(nseg, xmin, ymax, xmax, ymax); 211 Draw_Line_With_Transform(nseg, xmin, ymin, xmax, ymin); 212 Set_Pen_Color(Red); 213 214 Clear_Custom_Coordinate_Transform; 215 -- or... 216 -- plstransform(null, System.Null_Address); 217 end constriction2; 218 219 220 -- Vector plot of the gradient of a shielded potential (see example 9) 221 procedure potential is 222 nper : constant Integer := 100; 223 nlevel : constant Integer := 10; 224 nr : constant Integer := 20; 225 ntheta : constant Integer := 20; 226 227 eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i : Long_Float; 228 div1, div1i, div2, div2i : Long_Float; 229 r, theta, x, y, dz : Long_Float; 230 xmin, xmax, ymin, ymax, rmax, zmax, zmin : Long_Float; 231 u, v, z : Real_Matrix(0 .. nr - 1, 0 .. ntheta - 1); 232 px, py : Real_Vector(0 .. nper - 1); 233 clevel : Real_Vector(0 .. nlevel - 1); 234 cgrid2 : aliased Transformation_Data_Type_2 235 (x_Last => nr - 1, 236 y_Last => ntheta - 1); 237 238 function pow(x, y : Long_Float) return Long_Float is 239 Result : Long_Float := 1.0; 240 begin 241 for i in 1 .. Integer(y) loop 242 Result := Result * x; 243 end loop; 244 return Result; 245 end pow; 246 247 begin 248 -- Potential inside a conducting cylinder (or sphere) by method of images. 249 -- Charge 1 is placed at (d1, d1), with image charge at (d2, d2). 250 -- Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2). 251 -- Also put in smoothing term at small distances. 252 rmax := Long_Float(nr); 253 eps := 2.0; 254 255 q1 := 1.0; 256 d1 := rmax / 4.0; 257 258 q1i := - q1 * rmax / d1; 259 d1i := (rmax * rmax) / d1; 260 261 q2 := -1.0; 262 d2 := rmax / 4.0; 263 264 q2i := - q2 * rmax / d2; 265 d2i := (rmax * rmax) / d2; 266 267 for i in 0 .. nr - 1 loop 268 r := 0.5 + Long_Float(i); 269 for j in 0 .. ntheta - 1 loop 270 theta := 2.0 * pi / Long_Float(ntheta - 1) * (0.5 + Long_Float(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.0) + pow(y-d1, 2.0) + pow(eps, 2.0)); 276 div1i := sqrt(pow(x-d1i, 2.0) + pow(y-d1i, 2.0) + pow(eps, 2.0)); 277 div2 := sqrt(pow(x-d2, 2.0) + pow(y+d2, 2.0) + pow(eps, 2.0)); 278 div2i := sqrt(pow(x-d2i, 2.0) + pow(y+d2i, 2.0) + pow(eps, 2.0)); 279 280 z(i, j) := q1/div1 + q1i/div1i + q2/div2 + q2i/div2i; 281 u(i, j) := -q1*(x-d1)/pow(div1, 3.0) - q1i*(x-d1i)/pow(div1i, 3.00) 282 - q2*(x-d2)/pow(div2,3.0) - q2i*(x-d2i)/pow(div2i, 3.0); 283 v(i, j) := -q1*(y-d1)/pow(div1, 3.0) - q1i*(y-d1i)/pow(div1i, 3.00) 284 - q2*(y+d2)/pow(div2, 3.0) - q2i*(y+d2i)/pow(div2i, 3.0); 285 end loop; 286 end loop; 287 288 xmin := Matrix_Min(cgrid2.xg); 289 xmax := Matrix_Max(cgrid2.xg); 290 ymin := Matrix_Min(cgrid2.yg); 291 ymax := Matrix_Max(cgrid2.yg); 292 zmin := Matrix_Min(z); 293 zmax := Matrix_Max(z); 294 295 Set_Environment(xmin, xmax, ymin, ymax, Not_Justified, Linear_Box_Plus); 296 Write_Labels("(x)", "(y)", "#frPLplot Example 22 - potential gradient vector plot"); 297 298 -- Plot contours of the potential 299 dz := (zmax - zmin) / Long_Float(nlevel); 300 for i in clevel'range loop 301 clevel(i) := zmin + (Long_Float(i) + 0.5) * dz; 302 end loop; 303 Set_Pen_Color(Green); 304 Select_Line_Style(2); 305 Contour_Plot(z, 1, nr, 1, ntheta, clevel, Plot_Transformation_2'access, cgrid2'Address); 306 Select_Line_Style(1); 307 Set_Pen_Color(Red); 308 309 -- Plot the vectors of the gradient of the potential 310 Set_Pen_Color(Yellow); 311 Vector_Plot(u, v, 25.0, Plot_Transformation_2'access, cgrid2'Address); 312 Set_Pen_Color(Red); 313 314 -- Plot the perimeter of the cylinder 315 for i in px'range loop 316 theta := (2.0 * pi / Long_Float(nper - 1)) * Long_Float(i); 317 px(i) := rmax * cos(theta); 318 py(i) := rmax * sin(theta); 319 end loop; 320 Draw_Curve(px,py); 321 end potential; 322 323 324---------------------------------------------------------------------------- 325-- Generates several simple vector plots. 326---------------------------------------------------------------------------- 327begin 328 329 -- Parse and process command line arguments 330 Parse_Command_Line_Arguments(Parse_Full); 331 332 -- Initialize plplot 333 Initialize_PLplot; 334 335 circulation; 336 337 -- Set arrow style using arrow_x and arrow_y then plot using these arrows. 338 Set_Arrow_Style_For_Vector_Plots(arrow_x, arrow_y, False); 339 constriction(1); 340 341 -- Set arrow style using arrow2_x and arrow2_y then plot using these filled arrows. 342 Set_Arrow_Style_For_Vector_Plots(arrow2_x, arrow2_y, True); 343 constriction(2); 344 345 constriction2; 346 347 -- Reset arrow style to the default by passing two NULL arrays. 348 -- This line uses the awkward method of the C API to reset the default arrow style. 349 -- Set_Arrow_Style_For_Vector_Plots(System.Null_Address, System.Null_Address, False); 350 351 -- This method of resetting the default arrow style is a little more Ada-friendly... 352 -- plsvect; 353 354 -- ... as is this one which is identical but for name. 355 Reset_Vector_Arrow_Style; 356 357 potential; 358 359 End_PLplot; 360end xstandard22a; 361