chiark / gitweb /
X window size is a parameter; scale is a parameter; whether to draw original strip...
[moebius.git] / moebius.cc
1 /*
2  * Equation for a Moebius strip
3  */
4
5 #include <math.h>
6 #include <stdio.h>
7
8 #include "library.hh"
9 #include "moebius.hh"
10
11 Point MoebiusStrip::edgepoint(double t) {
12   double theta= (t-0.5)*4.0*M_PI;
13   double r= 1.0 - halfgap*(1.0 - sin(theta/2.0));
14   return Point(r*sin(theta),
15                -r*cos(theta),
16                halfbreadth*cos(theta/2.0));
17 }
18
19 Point MoebiusStrip::middlepoint(double t, double u) {
20   return edgepoint(t*0.5)*u + edgepoint((t+1)*0.5)*(1.0-u);
21 }
22
23 class Bezier {
24   double e,f,g,h;
25 public:
26   Bezier(double x0, double x1, double dx0, double dx1);
27   double operator()(double t) { return h + t*(g + t*(f + t*e)); }
28   void debug();
29 };
30
31 Bezier::Bezier(double x0, double x1, double dx0, double dx1) {
32 //  e= 2*x0 - 2*x1 + dx0 + dx1;
33 //  f= x1 - dx0 - x0 - e;
34   h= x0;
35   g= dx0;
36   e= g + 2*h + dx1 - 2*x1;
37   f= x1 - e - g - h;
38 }
39
40 void Bezier::debug() {
41   fprintf(stderr,"bz e %7.4f f %7.4f g %7.4f h %7.4f\n",e,f,g,h);
42 }
43
44 // The first end is at [sin(theta/2),-cos(theta/2),0]
45 // The second end is at [-theta/pi,0,sin(theta)]
46 // The first end is oriented towards [0,cos(theta),sin(theta)]
47 // The second end is oriented towards [0,-1,0]
48
49 Point MoebiusEnfoldment::edgepoint(double t) {
50   double theta= t*2.0*M_PI;
51   return Point(sin(theta),cos(theta),0);
52 }
53
54 Point MoebiusEnfoldment::middlepoint(double t, double u) {
55   if (t > bottomportion) {
56     t -= bottomportion;
57     t /= (1.0 - bottomportion);
58     double sizehere= sqrt(1-t*t);
59     return Point((u*2.0-1.0) * sizehere,
60                  t,
61                  sizehere * thickness * sin(u*2.0*M_PI));
62   } else {
63     t /= bottomportion;
64     double theta= (.5-u)*2*M_PI;
65     Bezier bx(sin(theta*.5), -theta/M_PI, 0, 0);
66     double ypushiness= (1-cos(theta))*2.5+1;
67 //        double fudge= (M_PI*sin(theta*.5)*cos(theta*.5))*(.5-u)*(.5-u)*4;
68     double fudge= (.5-u)*(.5-u)*4*cos(theta*.5);
69     Bezier by(-cos(theta*.5), 0,
70               cos(theta)*ypushiness + fudge*ypushiness,
71               ypushiness);
72 //by.debug();
73     Bezier bz(0, sin(theta), sin(theta), 0);
74     return Point( bx(t), by(t), thickness * bz(t) );
75   }
76 }