X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=moebius.cc;h=f169ab52520866f81c40703ce1c87dab3255504d;hb=HEAD;hp=0324e0d4a33a87a6934c04915d9ba77684d1b426;hpb=93bd402c99508c747f89c174a6d7e871def2a85e;p=moebius.git diff --git a/moebius.cc b/moebius.cc index 0324e0d..f169ab5 100644 --- a/moebius.cc +++ b/moebius.cc @@ -4,16 +4,12 @@ #include #include -#include #include "library.hh" #include "moebius.hh" -#include "parameter.hh" - -Parameter 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), @@ -50,6 +46,11 @@ void Bezier::debug() { // 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; @@ -57,13 +58,13 @@ Point MoebiusEnfoldment::middlepoint(double t, double u) { 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, @@ -74,86 +75,36 @@ Point MoebiusEnfoldment::middlepoint(double t, double u) { } } -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; }