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