1(*
2        Drawing "spirograph" curves - epitrochoids, cycolids, roulettes
3
4   Copyright (C) 2007  Arjen Markus
5   Copyright (C) 2008  Hezekiah M. Carty
6
7  This file is part of PLplot.
8
9  PLplot is free software; you can redistribute it and/or modify
10  it under the terms of the GNU Library General Public License as published
11  by the Free Software Foundation; either version 2 of the License, or
12  (at your option) any later version.
13
14  PLplot is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  GNU Library General Public License for more details.
18
19  You should have received a copy of the GNU Library General Public License
20  along with PLplot; if not, write to the Free Software
21  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22
23*)
24
25open Plplot
26
27let pi = atan 1.0 *. 4.0
28
29let cycloid () = () (* TODO *)
30
31(* Calculate greatest common divisor following pseudo-code for the
32   Euclidian algorithm at http://en.wikipedia.org/wiki/Euclidean_algorithm *)
33let rec gcd a b =
34  let a = abs a in
35  let b = abs b in
36  match b with
37  | 0 -> a
38  | _ -> gcd b (a mod b)
39
40let spiro params fill =
41  let npnt = 2000 in
42  let xcoord = Array.make (npnt + 1) 0.0 in
43  let ycoord = Array.make (npnt + 1) 0.0 in
44
45  (* Fill the coordinates *)
46
47  (* Proper termination of the angle loop, very near the beginning
48     point, see
49     http://mathforum.org/mathimages/index.php/Hypotrochoid *)
50  let windings =
51    int_of_float (abs_float params.(1)) /
52      gcd (int_of_float params.(0)) (int_of_float params.(1))
53  in
54  let steps = npnt / windings in
55  let dphi = 2.0 *. pi /. float_of_int steps in
56
57  (* This initialisation appears to be necessary, but I (AWI) don't understand why!  *)
58  let xmin = ref 0.0 in
59  let xmax = ref 0.0 in
60  let ymin = ref 0.0 in
61  let ymax = ref 0.0 in
62
63
64  for i = 0 to windings * steps do
65    let phi = float_of_int i *. dphi in
66    let phiw = (params.(0) -. params.(1)) /. params.(1) *. phi in
67    xcoord.(i) <- (params.(0) -. params.(1)) *. cos phi +. params.(2) *. cos phiw;
68    ycoord.(i) <- (params.(0) -. params.(1)) *. sin phi -. params.(2) *. sin phiw;
69
70    if i = 0 then (
71      xmin := xcoord.(i);
72      xmax := xcoord.(i);
73      ymin := ycoord.(i);
74      ymax := ycoord.(i);
75    )
76    else (
77    );
78    if !xmin > xcoord.(i) then xmin := xcoord.(i) else ();
79    if !xmax < xcoord.(i) then xmax := xcoord.(i) else ();
80    if !ymin > ycoord.(i) then ymin := ycoord.(i) else ();
81    if !ymax < ycoord.(i) then ymax := ycoord.(i) else ();
82  done;
83
84  let xrange_adjust = 0.15 *. (!xmax -. !xmin) in
85  let xmin = !xmin -. xrange_adjust in
86  let xmax = !xmax +. xrange_adjust in
87  let yrange_adjust = 0.15 *. (!ymax -. !ymin) in
88  let ymin = !ymin -. yrange_adjust in
89  let ymax = !ymax +. yrange_adjust in
90
91  plwind xmin xmax ymin ymax;
92
93  plcol0 1;
94  let xcoord' = Array.sub xcoord 0 (1 + steps * windings) in
95  let ycoord' = Array.sub ycoord 0 (1 + steps * windings) in
96  if fill then (
97    plfill xcoord' ycoord';
98  )
99  else (
100    plline xcoord' ycoord';
101  )
102
103let deg_to_rad x = x *. pi /. 180.0
104
105let arcs () =
106  let nseg = 8 in
107
108  let theta = ref 0.0 in
109  let dtheta = 360.0 /. float_of_int nseg in
110  plenv ~-.10.0 10.0 ~-.10.0 10.0 1 0;
111
112  (* Plot segments of circle in different colors *)
113  for i = 0 to nseg - 1 do
114    plcol0 (i mod 2 + 1);
115    plarc 0.0 0.0 8.0 8.0 !theta (!theta +. dtheta) 0.0 false;
116    theta := !theta +. dtheta;
117  done;
118
119  (* Draw several filled ellipses inside the circle at different
120     angles. *)
121  let a = 3.0 in
122  let b = a *. tan (deg_to_rad dtheta /. 2.0) in
123  theta := dtheta /. 2.0;
124  for i = 0 to nseg - 1 do
125    plcol0 (2 - i mod 2);
126    plarc (a *. cos (deg_to_rad !theta)) (a *. sin (deg_to_rad !theta)) a b 0.0 360.0 !theta true;
127    theta := !theta +. dtheta;
128  done;
129  ()
130
131(*--------------------------------------------------------------------------*\
132 * Generates two kinds of plots:
133 *   - construction of a cycloid (animated)
134 *   - series of epitrochoids and hypotrochoids
135\*--------------------------------------------------------------------------*)
136let () =
137  (* R, r, p, N *)
138  let params =
139    [|
140      [|21.0;  7.0;  7.0;  3.0|];  (* Deltoid *)
141      [|21.0;  7.0; 10.0;  3.0|];
142      [|21.0; -7.0; 10.0;  3.0|];
143      [|20.0;  3.0;  7.0; 20.0|];
144      [|20.0;  3.0; 10.0; 20.0|];
145      [|20.0; -3.0; 10.0; 20.0|];
146      [|20.0; 13.0;  7.0; 20.0|];
147      [|20.0; 13.0; 20.0; 20.0|];
148      [|20.0;-13.0; 20.0; 20.0|];
149    |]
150  in
151
152  (* plplot initialization *)
153  (* Parse and process command line arguments *)
154  plparseopts Sys.argv [PL_PARSE_FULL];
155
156  (* Initialize plplot *)
157  plinit ();
158
159  (* Illustrate the construction of a cycloid *)
160  cycloid ();
161
162  (* Loop over the various curves
163     First an overview, then all curves one by one *)
164  (* Three by three window *)
165  plssub 3 3;
166
167  for i = 0 to 8 do
168    pladv 0;
169    plvpor 0.0 1.0 0.0 1.0;
170    spiro params.(i) false;
171  done;
172
173  pladv 0;
174  (* One window per curve *)
175  plssub 1 1;
176
177  for i = 0 to 8 do
178    pladv 0;
179    plvpor 0.0 1.0 0.0 1.0;
180    spiro params.(i) false;
181  done;
182
183  (* Fill the curves *)
184  pladv 0;
185  (* One window per curve *)
186  plssub 1 1;
187
188  for i = 0 to 8 do
189    pladv 0;
190    plvpor 0.0 1.0 0.0 1.0;
191    spiro params.(i) true;
192  done;
193
194  (* Finally, an example to test out plarc capabilities *)
195  arcs ();
196
197  (* Don't forget to call plend() to finish off! *)
198  plend ();
199  ()
200
201