1#  Drawing "spirograph" curves - epitrochoids, cycolids, roulettes
2#
3#  Copyright (C) 2007  Arjen Markus
4#  Copyright (C) 2008  Andrew Ross
5#
6# This file is part of PLplot.
7#
8# PLplot is free software; you can redistribute it and/or modify
9# it under the terms of the GNU Library General Public License as published
10# by the Free Software Foundation; either version 2 of the License, or
11# (at your option) any later version.
12#
13# PLplot is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16# GNU Library General Public License for more details.
17#
18# You should have received a copy of the GNU Library General Public License
19# along with PLplot; if not, write to the Free Software
20# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21#
22#
23
24# --------------------------------------------------------------------------
25# main
26#
27# Generates two kinds of plots:
28#   - construction of a cycloid (animated)
29#   - series of epitrochoids and hypotrochoids
30# --------------------------------------------------------------------------
31
32proc x27 {{w loopback}} {
33
34  #     R, r, p, N
35  set params {
36    { 21.0   7.0   7.0   3.0  }
37    { 21.0   7.0  10.0   3.0  }
38    { 21.0  -7.0  10.0   3.0  }
39    { 20.0   3.0   7.0  20.0  }
40    { 20.0   3.0  10.0  20.0  }
41    { 20.0  -3.0  10.0  20.0  }
42    { 20.0  13.0   7.0  20.0  }
43    { 20.0  13.0  20.0  20.0  }
44    { 20.0 -13.0  20.0  20.0  }}
45
46  #  Illustrate the construction of a cycloid
47
48  cycloid $w
49
50  #  Loop over the various curves
51  #  First an overview, then all curves one by one
52  $w cmd pladv 0
53  $w cmd plssub 3 3
54
55  set fill 0
56
57  for { set i 0 } { $i < 9 } { incr i } {
58     $w cmd pladv 0
59     $w cmd plvpor 0.0  1.0  0.0  1.0
60     spiro $w [lindex $params $i] $fill
61  }
62  $w cmd pladv 0
63  $w cmd plssub 1 1
64
65  for { set i 0 } { $i < 9 } { incr i } {
66     $w cmd pladv 0
67     $w cmd plvpor 0.0  1.0  0.0  1.0
68     spiro $w [lindex $params $i] $fill
69  }
70
71  #  Fill the curves
72  set fill 1
73
74  $w cmd pladv 0
75  $w cmd plssub 1 1 ;# One window per curve
76
77  for { set i 0 } { $i < 9 } { incr i } {
78      $w cmd  pladv 0
79      $w cmd  plvpor 0.0 1.0 0.0 1.0
80      spiro $w [lindex $params $i] $fill
81  }
82
83  arcs $w
84}
85
86#--------------------------------------------------------------------------
87# Calculate greatest common divisor following pseudo-code for the
88# Euclidian algorithm at http://en.wikipedia.org/wiki/Euclidean_algorithm
89
90proc gcd {a b} {
91
92    set a [expr {int(abs($a))}]
93    set b [expr {int(abs($b))}]
94    while { $b != 0 } {
95        set t $b
96        set b [expr {$a % $b}]
97        set a $t
98    }
99    return $a
100}
101
102#  ===============================================================
103
104proc cycloid {w} {
105
106  #     TODO
107
108}
109
110#  ===============================================================
111
112proc spiro {w params fill} {
113
114  foreach {param1 param2 param3 param4} $params {break}
115
116  #     Fill the coordinates
117
118  set NPNT 2000
119
120  #     Proper termination of the angle loop very near the beginning
121  #     point, see
122  #     http://mathforum.org/mathimages/index.php/Hypotrochoid.
123  set windings [expr {int(abs($param2)/[gcd $param1 $param2])}]
124  set steps    [expr {int($NPNT/$windings)}]
125  set dphi     [expr {2.0*$::PLPLOT::PL_PI/double($steps)}]
126  # puts [ format "windings, steps, dphi = %d, %d, %f" $windings $steps $dphi ]
127
128  set n [expr {int($windings*$steps)+1}]
129
130  matrix xcoord f $n
131  matrix ycoord f $n
132
133  for { set i 0 } { $i < $n } { incr i } {
134     set phi  [expr {double($i) * $dphi}]
135     set phiw [expr {($param1-$param2)/$param2*$phi}]
136     xcoord $i = [expr {($param1-$param2)*cos($phi)+$param3*cos($phiw)}]
137     ycoord $i = [expr {($param1-$param2)*sin($phi)-$param3*sin($phiw)}]
138
139     if { $i == 0} {
140	set xmin [xcoord 0]
141	set xmax [xcoord 0]
142	set ymin [ycoord 0]
143	set ymax [ycoord 0]
144     }
145     if { $xmin > [xcoord $i] } { set xmin [xcoord $i] }
146     if { $xmax < [xcoord $i] } { set xmax [xcoord $i] }
147     if { $ymin > [ycoord $i] } { set ymin [ycoord $i] }
148     if { $ymax < [ycoord $i] } { set ymax [ycoord $i] }
149  }
150
151  set xrange_adjust [expr {0.15 * ($xmax - $xmin) }]
152  set xmin [expr {$xmin - $xrange_adjust }]
153  set xmax [expr {$xmax + $xrange_adjust }]
154  set yrange_adjust [expr {0.15 * ($ymax - $ymin) }]
155  set ymin [expr {$ymin - $yrange_adjust }]
156  set ymax [expr {$ymax + $yrange_adjust }]
157
158  $w cmd plwind $xmin $xmax $ymin $ymax
159
160  $w cmd plcol0 1
161
162  if { $fill } {
163      $w cmd plfill xcoord ycoord
164  } else {
165      $w cmd plline xcoord ycoord
166  }
167}
168
169proc arcs {w} {
170    set NSEG 8
171    set pi $::PLPLOT::PL_PI
172
173    set theta 0.0
174    set dtheta [expr {360.0 / $NSEG}]
175    $w cmd plenv -10.0 10.0 -10.0 10.0 1 0
176
177    # Plot segments of circle in different colors
178    for { set i 0 } { $i < $NSEG } {incr i} {
179        $w cmd plcol0 [expr {$i%2 + 1}]
180        $w cmd plarc 0.0 0.0 8.0 8.0 $theta [expr {$theta + $dtheta}] 0.0 0
181        set theta [expr {$theta + $dtheta}]
182    }
183
184    # Draw several filled ellipses inside the circle at different
185    # angles.
186    set a 3.0
187    set b [expr {$a * tan( ($dtheta/180.0*$pi)/2.0 )}]
188    set theta  [expr {$dtheta/2.0}]
189    for {set i 0} { $i < $NSEG } { incr i } {
190        $w cmd plcol0 [expr {2 - $i%2}]
191        $w cmd plarc [expr {$a*cos($theta/180.0*$pi)}] [expr {$a*sin($theta/180.0*$pi)}] $a $b 0.0 360.0 $theta 1
192	set theta [expr {$theta + $dtheta}]
193    }
194
195}
196