+/*
+ * Equation for a Moebius strip
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <iostream.h>
+
+#include "library.hh"
+#include "moebius.hh"
+#include "parameter.hh"
+
+Parameter<double> stiff("stiff","stiffness of tangent specs.",0.4,0.1,0.0,1.0);
+
+Point MoebiusStrip::edgepoint(double t) {
+ double theta= (t-0.5)*4.0*PI;
+ double r= 1.0 - halfgap*(1.0 - sin(theta/2.0));
+ return Point(r*sin(theta),
+ -r*cos(theta),
+ halfbreadth*cos(theta/2.0));
+}
+
+Point MoebiusStrip::middlepoint(double t, double u) {
+ return edgepoint(t*0.5)*u + edgepoint((t+1)*0.5)*(1.0-u);
+}
+
+class Bezier {
+ double e,f,g,h;
+public:
+ Bezier(double x0, double x1, double dx0, double dx1);
+ double operator()(double t) { return h + t*(g + t*(f + t*e)); }
+ void debug();
+};
+
+Bezier::Bezier(double x0, double x1, double dx0, double dx1) {
+// e= 2*x0 - 2*x1 + dx0 + dx1;
+// f= x1 - dx0 - x0 - e;
+ h= x0;
+ g= dx0;
+ e= g + 2*h + dx1 - 2*x1;
+ f= x1 - e - g - h;
+}
+
+void Bezier::debug() {
+ fprintf(stderr,"bz e %7.4f f %7.4f g %7.4f h %7.4f\n",e,f,g,h);
+}
+
+// The first end is at [sin(theta/2),-cos(theta/2),0]
+// The second end is at [-theta/pi,0,sin(theta)]
+// The first end is oriented towards [0,cos(theta),sin(theta)]
+// The second end is oriented towards [0,-1,0]
+
+Point MoebiusEnfoldment::middlepoint(double t, double u) {
+ if (t > bottomportion) {
+ t -= bottomportion;
+ t /= (1.0 - bottomportion);
+ double sizehere= sqrt(1-t*t);
+ return Point((u*2.0-1.0) * sizehere,
+ t,
+ sizehere * thickness * sin(u*2.0*PI));
+ } else {
+ t /= bottomportion;
+ double theta= (.5-u)*2*PI;
+ Bezier bx(sin(theta*.5), -theta/PI, 0, 0);
+ double ypushiness= (1-cos(theta))*2.5+1;
+// double fudge= (PI*sin(theta*.5)*cos(theta*.5))*(.5-u)*(.5-u)*4;
+ double fudge= (.5-u)*(.5-u)*4*cos(theta*.5);
+ Bezier by(-cos(theta*.5), 0,
+ cos(theta)*ypushiness + fudge*ypushiness,
+ ypushiness);
+//by.debug();
+ Bezier bz(0, sin(theta), sin(theta), 0);
+ return Point( bx(t), by(t), thickness * bz(t) );
+ }
+}
+
+double smoothinterp(double t, const double values[], const double rates[]) {
+ const int ncells=4;
+ int n1= (int)floor(t*ncells);
+ int n2= n1+1;
+ double r= t*ncells-n1;
+ double v1= values[n1];
+ double v2= values[n2];
+ double r1= rates[n1]/ncells;
+ double r2= rates[n2]/ncells;
+ double dv= v2-v1;
+ double a= r1+r2-2*dv;
+ double b= 3*dv-2*r1-r2;
+ double c= r1;
+ double vr= a*r*r + b*r*r + c*r;
+ double v= vr + v1;
+// printf("smoothinterp %7.6f %7.6f %g(%g)..%g(%g) %g\n",t,r,v1,r1,v2,r2,v);
+ return v;
+}
+
+Point MoebiusNewEnfoldment::middlepoint(double t, double u) {
+ const double agammas[]= {
+ /* angle at rim for surface near edge, for values of t:
+ * 0.0 .25 0.5 0.75 1.0, not usu. used */
+ PI, 0.75*PI, 0.5*PI, -0.25*PI, -PI, -PI
+ };
+ const double agammarates[]= {
+ /* rates of change of gamma with t */
+ -PI, -PI, -PI, -3*PI, -3*PI, -PI
+ };
+ const double alambdas[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
+ const double alambdarates[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
+ const double bgammas[]= {
+ /* angle at middle (u=0.5) in direction of decreasing u, for values of t:
+ * 0.0 .25 0.5 0.75 1.0, not usu. used */
+ 0, 0.5*PI, 0.5*PI, 0.5*PI, PI, PI
+ };
+ const double bgammarates[]= {
+ /* rates of change of gamma with t */
+ 2*PI, 0, 0, 0, 2*PI, 2*PI
+ };
+ const double blambdas[]= {
+ /* shear along middle (u=0.5) for values of t:
+ * 0.0 .25 0.5 0.75 1.0, not usu. used */
+ 0.0, 0.5, 0.5, 0.5, 0, 0
+ };
+ const double blambdarates[]= { 0.0, 0.0, 0.0, 0.0, 0.0 };
+ if (u > 0.5) {
+ Point p= middlepoint(1.0-t,1.0-u);
+ return Point(-p[0],p[1],-p[2]);
+ }
+ double aangle= PI*t;
+ double asine= sin(aangle);
+ double acosine= cos(aangle);
+ double agamma= smoothinterp(t,agammas,agammarates);
+ double arotnl= smoothinterp(t,alambdas,alambdarates);
+ double aradial= cos(agamma);
+ double aaxial= sin(agamma);
+ Point a(asine, acosine, 0);
+ Point da(arotnl*acosine+aradial*asine,aradial*acosine-arotnl*asine,-aaxial);
+ double bangle= 2*PI*t;
+ double bsine= sin(bangle);
+ double bcosine= cos(bangle);
+ double bgamma= smoothinterp(t,bgammas,bgammarates);
+ double brotnl= smoothinterp(t,blambdas,blambdarates);
+ double bradial= cos(bgamma);
+ double baxial= sin(bgamma);
+ Point b(0, bcosine-1.0, -bsine);
+ Point db(baxial,bradial*bcosine-brotnl*bsine,brotnl*bcosine+bradial*bsine);
+ u*=2;
+ double ab= (a-b).magnitude();
+ double adscale= ab*stiff/da.magnitude();
+ double bdscale= ab*stiff/db.magnitude();
+ double facta= (1-u)*(1-u)*(1-u);
+ double factb= u*u*u;
+ double factan= 3*(1-u)*(1-u)*u;
+ double factbn= 3*(1-u)*u*u;
+ Point dan= a+da*adscale;
+ Point dbn= b+db*bdscale;
+ Point result= a*facta + dan*factan + dbn*factbn + b*factb;
+// printf("t=%7.6f u=%7.6f a=%-30s dan=%-30s dbn=%-30s b=%-30s\n",
+// t,u, a.printing(), dan.printing(), dbn.printing(), b.printing());
+ return result;
+}