#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 theta= (t-0.5)*4.0*M_PI;
double r= 1.0 - halfgap*(1.0 - sin(theta/2.0));
return Point(r*sin(theta),
-r*cos(theta),
// The first end is oriented towards [0,cos(theta),sin(theta)]
// The second end is oriented towards [0,-1,0]
+Point MoebiusEnfoldmentAny::edgepoint(double t) {
+ double theta= t*2.0*M_PI;
+ return Point(sin(theta),cos(theta),0);
+}
+
Point MoebiusEnfoldment::middlepoint(double t, double u) {
if (t > bottomportion) {
t -= bottomportion;
double sizehere= sqrt(1-t*t);
return Point((u*2.0-1.0) * sizehere,
t,
- sizehere * thickness * sin(u*2.0*PI));
+ sizehere * thickness * sin(u*2.0*M_PI));
} else {
t /= bottomportion;
- double theta= (.5-u)*2*PI;
- Bezier bx(sin(theta*.5), -theta/PI, 0, 0);
+ double theta= (.5-u)*2*M_PI;
+ Bezier bx(sin(theta*.5), -theta/M_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= (M_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,
}
}
-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 MoebiusEnfoldmentNew::middlepoint(double t, double u) {
+ if (u <= 0.5) {
+ if (t <= 0.5) {
+ double sin_pi_t= sin(M_PI*t);
+ double cos_pi_t= cos(M_PI*t);
+ double sin_2pi_t= sin(2.0*M_PI*t);
+ double cos_2pi_t= cos(2.0*M_PI*t);
+
+ double prop_circ= 0.5+0.5*cos(2.0*M_PI*u);
+ // double prop_circ= 1.0-2.0*u;
-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]);
+ Point p_edge(-sin_pi_t,cos_pi_t,0);
+
+ double alpha_circ= (2.0*M_PI)*(2.0*t)*(2.0*u);
+ Point v_edge1(cos_2pi_t*sin_pi_t,-cos_2pi_t*cos_pi_t,sin_2pi_t);
+ Point v_edge2(sin_2pi_t*cos_pi_t,-sin_2pi_t*sin_pi_t,-cos_2pi_t);
+ Point p_circ_cent= p_edge + v_edge2*(2.0*t);
+ Point p_circ_pt= p_circ_cent + (v_edge1*sin(alpha_circ) + v_edge2*-cos(alpha_circ))*(2.0*t);
+ p_circ_pt[1]= p_edge[1];
+
+ Point v_line(-cos_pi_t,0,sin_pi_t);
+ Point p_line_start(0,1.0-2.0*t,0);
+ Point p_line_pt= p_line_start + v_line*(1.0-2.0*u)*(M_PI*t);
+
+ return p_circ_pt*prop_circ + p_line_pt*(1.0-prop_circ);
+ } else {
+ return middlepoint(0.5,u) + Point(0,0.5-t,0);
+ }
+ } else {
+ Point p_mirror= middlepoint(t,1.0-u);
+ return Point(-p_mirror[0], p_mirror[1], -p_mirror[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;
}
~ian / moebius.git / blobdiff