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 MoebiusEnfoldmentAny::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 }
  77
  78 Point MoebiusEnfoldmentNew::middlepoint(double t, double u) {
  79   if (u <= 0.5) {
  80     if (t <= 0.5) {
  81       double sin_pi_t= sin(M_PI*t);
  82       double cos_pi_t= cos(M_PI*t);
  83       double sin_2pi_t= sin(2.0*M_PI*t);
  84       double cos_2pi_t= cos(2.0*M_PI*t);
  85
  86       double prop_circ= 0.5+0.5*cos(2.0*M_PI*u);
  87       //      double prop_circ= 1.0-2.0*u;
  88
  89       Point p_edge(-sin_pi_t,cos_pi_t,0);
  90
  91       double alpha_circ= (2.0*M_PI)*(2.0*t)*(2.0*u);
  92       Point v_edge1(cos_2pi_t*sin_pi_t,-cos_2pi_t*cos_pi_t,sin_2pi_t);
  93       Point v_edge2(sin_2pi_t*cos_pi_t,-sin_2pi_t*sin_pi_t,-cos_2pi_t);
  94       Point p_circ_cent= p_edge + v_edge2*(2.0*t);
  95       Point p_circ_pt= p_circ_cent + (v_edge1*sin(alpha_circ) + v_edge2*-cos(alpha_circ))*(2.0*t);
  96       p_circ_pt[1]= p_edge[1];
  97
  98       Point v_line(-cos_pi_t,0,sin_pi_t);
  99       Point p_line_start(0,1.0-2.0*t,0);
 100       Point p_line_pt= p_line_start + v_line*(1.0-2.0*u)*(M_PI*t);
 101
 102       return p_circ_pt*prop_circ + p_line_pt*(1.0-prop_circ);
 103     } else {
 104       return middlepoint(0.5,u) + Point(0,0.5-t,0);
 105     }
 106   } else {
 107     Point p_mirror= middlepoint(t,1.0-u);
 108     return Point(-p_mirror[0], p_mirror[1], -p_mirror[2]);
 109   }
 110 }