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*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::middlepoint(double t, double u) {
  50   if (t > bottomportion) {
  51     t -= bottomportion;
  52     t /= (1.0 - bottomportion);
  53     double sizehere= sqrt(1-t*t);
  54     return Point((u*2.0-1.0) * sizehere,
  55                  t,
  56                  sizehere * thickness * sin(u*2.0*PI));
  57   } else {
  58     t /= bottomportion;
  59     double theta= (.5-u)*2*PI;
  60     Bezier bx(sin(theta*.5), -theta/PI, 0, 0);
  61     double ypushiness= (1-cos(theta))*2.5+1;
  62 //        double fudge= (PI*sin(theta*.5)*cos(theta*.5))*(.5-u)*(.5-u)*4;
  63     double fudge= (.5-u)*(.5-u)*4*cos(theta*.5);
  64     Bezier by(-cos(theta*.5), 0,
  65               cos(theta)*ypushiness + fudge*ypushiness,
  66               ypushiness);
  67 //by.debug();
  68     Bezier bz(0, sin(theta), sin(theta), 0);
  69     return Point( bx(t), by(t), thickness * bz(t) );
  70   }
  71 }
  72
  73 double lininterp(double t, const double array[]) {
  74   const int ncells=4;
  75   int n1= (int)floor(t*ncells);
  76   int n2= n1+1;
  77   double r= t*ncells-n1;
  78   double v1= array[n1];
  79   double v2= array[n2];
  80   return v1*(1.0-r) + v2*r;
  81 }
  82
  83 Point MoebiusNewEnfoldment::middlepoint(double t, double u) {
  84   const double agammas[]= {
  85     /* angle at rim for surface near edge, for values of t:
  86      * 0.0   .25      0.5      0.75   1.0,  not usu. used */
  87        PI, 0.75*PI, 0.5*PI, -0.25*PI, -PI,    -PI
  88   };
  89   const double alambdas[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
  90   const double bgammas[]= {
  91     /* angle at middle (u=0.5) in direction of decreasing u, for values of t:
  92      * 0.0   .25     0.5     0.75    1.0,  not usu. used */
  93        0,   0.5*PI, 0.5*PI, 0.5*PI, 0.5*PI, 0.5*PI
  94   };
  95   const double blambdas[]= {
  96     /* shear along middle (u=0.5) for values of t:
  97      * 0.0   .25   0.5  0.75   1.0,  not usu. used */
  98        0.0,  0.5,  0.5,  0.5,  1.0,  1.0
  99   };
 100   if (u > 0.5) {
 101     Point p= middlepoint(1.0-t,1.0-u);
 102     return Point(-p[0],p[1],-p[2]);
 103   }
 104   double aangle= PI*t;
 105   double asine= sin(aangle);
 106   double acosine= cos(aangle);
 107   double agamma= lininterp(t,agammas);
 108   double arotnl= lininterp(t,alambdas);
 109   double aradial= cos(agamma);
 110   double aaxial= sin(agamma);
 111   Point a(asine,acosine,0);
 112   Point da(arotnl*acosine+aradial*asine,aradial*acosine-arotnl*asine,-aaxial);
 113   double bangle= 2*PI*t;
 114   double bsine= sin(bangle);
 115   double bcosine= cos(bangle);
 116   double bgamma= lininterp(t,bgammas);
 117   double brotnl= lininterp(t,blambdas);
 118   double bradial= cos(bgamma);
 119   double baxial= sin(bgamma);
 120   Point b(0,bcosine-1.0,bsine);
 121   Point db(baxial,bradial*bcosine-brotnl*bsine,brotnl*bcosine+bradial*bsine);
 122     Point tp;
 123     tp=a; a=b; b=tp;
 124     tp=da; da=db; db=tp;
 125   double ab= (a-b).magnitude();
 126   double adscale= ab*0.2/da.magnitude();
 127   double bdscale= ab*0.2/db.magnitude();
 128   u*=2;
 129   Bezier x(a[0],b[0],da[0]*adscale,db[0]*bdscale);
 130   Bezier y(a[1],b[1],da[1]*adscale,db[1]*bdscale);
 131   Bezier z(a[2],b[2],da[2]*adscale,db[2]*bdscale);
 132   return Point(x(u),y(u),z(u));
 133   return a*u+b*(1.0-u);
 134 }