1-- Drawing "spirograph" curves - epitrochoids, cycolids, roulettes
2
3-- Copyright (C) 2008 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    Ada.Numerics,
23    Ada.Numerics.Long_Elementary_Functions,
24    PLplot_Auxiliary,
25    PLplot_Traditional;
26use
27    Ada.Numerics,
28    Ada.Numerics.Long_Elementary_Functions,
29    PLplot_Auxiliary,
30    PLplot_Traditional;
31
32------------------------------------------------------------------------------
33-- Generates two kinds of plots:
34--   - construction of a cycloid (animated)
35--   - series of epitrochoids and hypotrochoids
36------------------------------------------------------------------------------
37
38procedure xtraditional27a is
39    -- R, r, p, N
40    -- R and r should be integers to give correct termination of the
41    -- angle loop using gcd.
42    -- N.B. N is just a place holder since it is no longer used
43    -- (because we now have proper termination of the angle loop).
44    params : Real_Matrix(0 .. 8, 0 .. 3) :=
45      ((21.0,  7.0,  7.0,  3.0),
46       (21.0,  7.0, 10.0,  3.0),
47       (21.0, -7.0, 10.0,  3.0),
48       (20.0,  3.0,  7.0, 20.0),
49       (20.0,  3.0, 10.0, 20.0),
50       (20.0, -3.0, 10.0, 20.0),
51       (20.0, 13.0,  7.0, 20.0),
52       (20.0, 13.0, 20.0, 20.0),
53       (20.0,-13.0, 20.0, 20.0));
54
55    fill : Boolean;
56
57    -- To understand why spiro is written this way you need to understand the
58    -- C code from which this was derived. In the main C program, params
59    -- is a two-dimensional array with 9 rows numbered 0 .. 8 and 4 columns
60    -- numbered 0 .. 3. When spiro is called, it is passed the _address_ of the
61    -- element of params's ith row, 0th column--nothing else. Then, inside spiro,
62    -- the corresponding entity (also called params!) appears as a
63    -- _one-dimensional_ array whose 0th element shares the same address as what
64    -- was passed from the main program. So memory locations starting there,
65    -- and numbered from 0, represent the 1D array equivalent to the ith row of
66    -- params in the main program. Wilma, call Barney--we're programming a
67    -- micaprocessor here.
68    procedure spiro(params : Real_Matrix; row : Integer; fill : Boolean) is
69        NPNT : constant Integer := 2000;
70        xcoord, ycoord : Real_Vector(0 .. NPNT);
71        windings : Integer;
72        steps    : Integer;
73        phi      : Long_Float;
74        phiw     : Long_Float;
75        dphi     : Long_Float;
76        xmin     : Long_Float;
77        xmax     : Long_Float;
78        xrange_adjust : Long_Float;
79        ymin     : Long_Float;
80        ymax     : Long_Float;
81        yrange_adjust : Long_Float;
82
83        function Trunc(a : Long_Float) return Integer renames PLplot_Auxiliary.Trunc;
84
85        -- Calculate greatest common divisor following pseudo-code for the
86        -- Euclidian algorithm at http://en.wikipedia.org/wiki/Euclidean_algorithm
87        function gcd(a, b : Integer) return Integer is
88            t : Integer;
89            aa : Integer := a;
90            bb : Integer := b;
91        begin
92            aa := abs(aa);
93            bb := abs(bb);
94            while bb /= 0 loop
95                t := bb;
96                bb := aa mod bb;
97                aa := t;
98            end loop;
99            return aa;
100        end gcd;
101
102    begin -- spiro
103        -- Fill the coordinates.
104        -- Proper termination of the angle loop very near the beginning
105        -- point, see http://mathforum.org/mathimages/index.php/Hypotrochoid
106        windings := Trunc(abs(params(row, 1)) /
107            Long_Float(gcd(Trunc(params(row, 0)), Trunc(params(row, 1)))));
108        steps    := NPNT / windings;
109        dphi     := 2.0 * pi / Long_Float(steps);
110        for i in 0 .. windings * steps loop
111            phi       := Long_Float(i) * dphi;
112            phiw      := (params(row, 0) - params(row, 1)) / params(row, 1) * phi;
113            xcoord(i) := (params(row, 0)-params(row, 1))*cos(phi)+params(row, 2)*cos(phiw);
114            ycoord(i) := (params(row, 0)-params(row, 1))*sin(phi)-params(row, 2)*sin(phiw);
115            if i = 0 then
116                xmin := xcoord(i);
117                xmax := xcoord(i);
118                ymin := ycoord(i);
119                ymax := ycoord(i);
120            end if;
121            if xmin > xcoord(i) then xmin := xcoord(i); end if;
122            if xmax < xcoord(i) then xmax := xcoord(i); end if;
123            if ymin > ycoord(i) then ymin := ycoord(i); end if;
124            if ymax < ycoord(i) then ymax := ycoord(i); end if;
125        end loop;
126
127        xrange_adjust := 0.15 * (xmax - xmin);
128        xmin := xmin - xrange_adjust;
129        xmax := xmax + xrange_adjust;
130        yrange_adjust := 0.15 * (ymax - ymin);
131        ymin := ymin - yrange_adjust;
132        ymax := ymax + yrange_adjust;
133
134        plwind(xmin, xmax, ymin, ymax);
135        plcol0(1);
136
137        declare
138            xcoord_local, ycoord_local : Real_Vector(0 .. steps * windings);
139        begin
140            xcoord_local := xcoord(0 .. steps * windings);
141            ycoord_local := ycoord(0 .. steps * windings);
142            if fill then
143                plfill(xcoord_local, ycoord_local);
144            else
145                plline(xcoord_local, ycoord_local);
146            end if;
147        end;
148    end spiro;
149
150    procedure cycloid is
151    begin
152        null; -- TODO
153    end cycloid;
154
155    procedure arcs is
156       NSEG : constant Integer := 8;
157       theta : Long_Float;
158       dtheta : Long_Float;
159       a : Long_Float;
160       b : Long_Float;
161    begin
162
163       theta := 0.0;
164       dtheta := 360.0/Long_Float(NSEG);
165       plenv( -10.0, 10.0, -10.0, 10.0, 1, 0 );
166
167       -- Plot segments of circle in different colors
168       for i in 0 .. NSEG-1 loop
169          plcol0( (i mod 2) + 1 );
170          plarc(0.0, 0.0, 8.0, 8.0, theta, theta + dtheta, 0.0, False);
171          theta := theta + dtheta;
172       end loop;
173
174       -- Draw several filled ellipses inside the circle at different
175       -- angles.
176       a := 3.0;
177       b := a * tan( (dtheta/180.0*pi)/2.0 );
178       theta := dtheta/2.0;
179       for i in 0 .. NSEG-1 loop
180          plcol0( 2 - (i mod 2 ) );
181          plarc( a*cos(theta/180.0*pi), a*sin(theta/180.0*pi), a, b, 0.0, 360.0, theta, True);
182          theta := theta + dtheta;
183       end loop;
184
185    end arcs;
186
187begin
188    -- Parse and process command line arguments
189    plparseopts(PL_PARSE_FULL);
190
191    -- Initialize plplot
192    plinit;
193
194    -- Illustrate the construction of a cycloid
195    cycloid;
196
197    -- Loop over the various curves.
198    -- First an overview, then all curves one by one
199    plssub(3, 3) ; -- Three by three window
200
201    -- Overview
202    fill := False;
203    for i in params'range(1) loop
204        pladv(0);
205        plvpor(0.0, 1.0, 0.0, 1.0);
206        spiro(params, i, fill);
207    end loop;
208
209    -- Don't fill the curves.
210    pladv(0) ;
211    plssub(1, 1) ; -- One window per curve
212    for i in params'range(1) loop
213        pladv(0) ;
214        plvpor(0.0, 1.0, 0.0, 1.0);
215        spiro(params, i, fill);
216    end loop;
217
218    -- Fill the curves
219    fill := True;
220    pladv(0);
221    plssub(1, 1); -- One window per curve
222    for i in params'range(1) loop
223        pladv(0) ;
224        plvpor(0.0, 1.0, 0.0, 1.0);
225        spiro(params, i, fill);
226    end loop;
227
228    arcs;
229
230    -- Don't forget to call plend to finish off!
231    plend;
232end xtraditional27a;
233