moebius.cc

   1 /*
   2  * Equation for a Moebius strip
   3  */
   4
   5 #include <math.h>
   6 #include <stdio.h>
   7 #include <iostream.h>
   8
   9 #include "library.hh"
  10 #include "moebius.hh"
  11 #include "parameter.hh"
  12
  13 Parameter<double> stiff("stiff","stiffness of tangent specs.",0.4,0.1,0.0,1.0);
  14
  15 Point MoebiusStrip::edgepoint(double t) {
  16   double theta= (t-0.5)*4.0*PI;
  17   double r= 1.0 - halfgap*(1.0 - sin(theta/2.0));
  18   return Point(r*sin(theta),
  19                -r*cos(theta),
  20                halfbreadth*cos(theta/2.0));
  21 }
  22
  23 Point MoebiusStrip::middlepoint(double t, double u) {
  24   return edgepoint(t*0.5)*u + edgepoint((t+1)*0.5)*(1.0-u);
  25 }
  26
  27 class Bezier {
  28   double e,f,g,h;
  29 public:
  30   Bezier(double x0, double x1, double dx0, double dx1);
  31   double operator()(double t) { return h + t*(g + t*(f + t*e)); }
  32   void debug();
  33 };
  34
  35 Bezier::Bezier(double x0, double x1, double dx0, double dx1) {
  36 //  e= 2*x0 - 2*x1 + dx0 + dx1;
  37 //  f= x1 - dx0 - x0 - e;
  38   h= x0;
  39   g= dx0;
  40   e= g + 2*h + dx1 - 2*x1;
  41   f= x1 - e - g - h;
  42 }
  43
  44 void Bezier::debug() {
  45   fprintf(stderr,"bz e %7.4f f %7.4f g %7.4f h %7.4f\n",e,f,g,h);
  46 }
  47
  48 // The first end is at [sin(theta/2),-cos(theta/2),0]
  49 // The second end is at [-theta/pi,0,sin(theta)]
  50 // The first end is oriented towards [0,cos(theta),sin(theta)]
  51 // The second end is oriented towards [0,-1,0]
  52
  53 Point MoebiusEnfoldment::middlepoint(double t, double u) {
  54   if (t > bottomportion) {
  55     t -= bottomportion;
  56     t /= (1.0 - bottomportion);
  57     double sizehere= sqrt(1-t*t);
  58     return Point((u*2.0-1.0) * sizehere,
  59                  t,
  60                  sizehere * thickness * sin(u*2.0*PI));
  61   } else {
  62     t /= bottomportion;
  63     double theta= (.5-u)*2*PI;
  64     Bezier bx(sin(theta*.5), -theta/PI, 0, 0);
  65     double ypushiness= (1-cos(theta))*2.5+1;
  66 //        double fudge= (PI*sin(theta*.5)*cos(theta*.5))*(.5-u)*(.5-u)*4;
  67     double fudge= (.5-u)*(.5-u)*4*cos(theta*.5);
  68     Bezier by(-cos(theta*.5), 0,
  69               cos(theta)*ypushiness + fudge*ypushiness,
  70               ypushiness);
  71 //by.debug();
  72     Bezier bz(0, sin(theta), sin(theta), 0);
  73     return Point( bx(t), by(t), thickness * bz(t) );
  74   }
  75 }
  76
  77 double smoothinterp(double t, const double values[], const double rates[]) {
  78   const int ncells=4;
  79   int n1= (int)floor(t*ncells);
  80   int n2= n1+1;
  81   double r= t*ncells-n1;
  82   double v1= values[n1];
  83   double v2= values[n2];
  84   double r1= rates[n1]/ncells;
  85   double r2= rates[n2]/ncells;
  86   double dv= v2-v1;
  87   double a= r1+r2-2*dv;
  88   double b= 3*dv-2*r1-r2;
  89   double c= r1;
  90   double vr= a*r*r + b*r*r + c*r;
  91   double v= vr + v1;
  92 //  printf("smoothinterp %7.6f  %7.6f %g(%g)..%g(%g) %g\n",t,r,v1,r1,v2,r2,v);
  93   return v;
  94 }
  95
  96 Point MoebiusNewEnfoldment::middlepoint(double t, double u) {
  97   const double agammas[]= {
  98     /* angle at rim for surface near edge, for values of t:
  99      * 0.0   .25      0.5      0.75   1.0,  not usu. used */
 100        PI, 0.75*PI, 0.5*PI, -0.25*PI, -PI,    -PI
 101   };
 102   const double agammarates[]= {
 103     /* rates of change of gamma with t */
 104        -PI,   -PI,   -PI,    -3*PI,  -3*PI,   -PI
 105   };
 106   const double alambdas[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
 107   const double alambdarates[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
 108   const double bgammas[]= {
 109     /* angle at middle (u=0.5) in direction of decreasing u, for values of t:
 110      * 0.0   .25     0.5     0.75    1.0,  not usu. used */
 111        0,   0.5*PI, 0.5*PI, 0.5*PI,   PI,   PI
 112   };
 113   const double bgammarates[]= {
 114     /* rates of change of gamma with t */
 115        2*PI,  0,      0,      0,     2*PI, 2*PI
 116   };
 117   const double blambdas[]= {
 118     /* shear along middle (u=0.5) for values of t:
 119      * 0.0   .25   0.5  0.75   1.0,  not usu. used */
 120        0.0,  0.5,  0.5,  0.5,   0,    0
 121   };
 122   const double blambdarates[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
 123   if (u > 0.5) {
 124     Point p= middlepoint(1.0-t,1.0-u);
 125     return Point(-p[0],p[1],-p[2]);
 126   }
 127   double aangle= PI*t;
 128   double asine= sin(aangle);
 129   double acosine= cos(aangle);
 130   double agamma= smoothinterp(t,agammas,agammarates);
 131   double arotnl= smoothinterp(t,alambdas,alambdarates);
 132   double aradial= cos(agamma);
 133   double aaxial= sin(agamma);
 134   Point a(asine, acosine, 0);
 135   Point da(arotnl*acosine+aradial*asine,aradial*acosine-arotnl*asine,-aaxial);
 136   double bangle= 2*PI*t;
 137   double bsine= sin(bangle);
 138   double bcosine= cos(bangle);
 139   double bgamma= smoothinterp(t,bgammas,bgammarates);
 140   double brotnl= smoothinterp(t,blambdas,blambdarates);
 141   double bradial= cos(bgamma);
 142   double baxial= sin(bgamma);
 143   Point b(0, bcosine-1.0, -bsine);
 144   Point db(baxial,bradial*bcosine-brotnl*bsine,brotnl*bcosine+bradial*bsine);
 145   u*=2;
 146   double ab= (a-b).magnitude();
 147   double adscale= ab*stiff/da.magnitude();
 148   double bdscale= ab*stiff/db.magnitude();
 149   double facta= (1-u)*(1-u)*(1-u);
 150   double factb= u*u*u;
 151   double factan= 3*(1-u)*(1-u)*u;
 152   double factbn=  3*(1-u)*u*u;
 153   Point dan= a+da*adscale;
 154   Point dbn= b+db*bdscale;
 155   Point result= a*facta + dan*factan + dbn*factbn + b*factb;
 156 //  printf("t=%7.6f u=%7.6f a=%-30s dan=%-30s dbn=%-30s b=%-30s\n",
 157 //         t,u, a.printing(), dan.printing(), dbn.printing(), b.printing());
 158   return result;
 159 }