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_Standard; 26use 27 Ada.Numerics, 28 Ada.Numerics.Long_Elementary_Functions, 29 PLplot_Auxiliary, 30 PLplot_Standard; 31 32------------------------------------------------------------------------------ 33-- Generates two kinds of plots: 34-- - construction of a cycloid (animated) 35-- - series of epitrochoids and hypotrochoids 36------------------------------------------------------------------------------ 37 38procedure xstandard27a 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 Set_Viewport_World(xmin, xmax, ymin, ymax); 135 Set_Pen_Color(Red); 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 Fill_Polygon(xcoord_local, ycoord_local); 144 else 145 Draw_Curve(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 Set_Environment( -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 Set_Pen_Color( (i mod 2) + 1 ); 170 Draw_Arc(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 Set_Pen_Color( 2 - (i mod 2 ) ); 181 Draw_Arc( 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 Parse_Command_Line_Arguments(Parse_Full); 190 191 -- Initialize plplot 192 Initialize_PLplot; 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 Set_Number_Of_Subpages(3, 3) ; -- Three by three window 200 201 -- Overview 202 fill := False; 203 for i in params'range(1) loop 204 Advance_To_Subpage(Next_Subpage); 205 Set_Viewport_Normalized( 0.0, 1.0, 0.0, 1.0); 206 spiro(params, i, fill); 207 end loop; 208 209 -- Don't fill the curves. 210 Advance_To_Subpage(Next_Subpage); 211 Set_Number_Of_Subpages(1, 1) ; -- One window per curve 212 for i in params'range(1) loop 213 Advance_To_Subpage(Next_Subpage); 214 Set_Viewport_Normalized(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 Advance_To_Subpage(Next_Subpage); 221 Set_Number_Of_Subpages(1, 1); -- One window per curve 222 for i in params'range(1) loop 223 Advance_To_Subpage(Next_Subpage) ; 224 Set_Viewport_Normalized( 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 End_PLplot to finish off! 231 End_PLplot; 232end xstandard27a; 233