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