/*
 * Equation for a Moebius strip
 */

#include <math.h>
#include <stdio.h>

#include "library.hh"
#include "moebius.hh"

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 lininterp(double t, const double array[]) {
  const int ncells=4;
  int n1= (int)floor(t*ncells);
  int n2= n1+1;
  double r= t*ncells-n1;
  double v1= array[n1];
  double v2= array[n2];
  return v1*(1.0-r) + v2*r;
}

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 alambdas[]= { 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, 0.5*PI, 0.5*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,  1.0,  1.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= lininterp(t,agammas);
  double arotnl= lininterp(t,alambdas);
  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= lininterp(t,bgammas);
  double brotnl= lininterp(t,blambdas);
  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);
    Point tp;
    tp=a; a=b; b=tp;
    tp=da; da=db; db=tp;
  double ab= (a-b).magnitude();
  double adscale= ab*0.2/da.magnitude();
  double bdscale= ab*0.2/db.magnitude();
  u*=2;
  Bezier x(a[0],b[0],da[0]*adscale,db[0]*bdscale);
  Bezier y(a[1],b[1],da[1]*adscale,db[1]*bdscale);
  Bezier z(a[2],b[2],da[2]*adscale,db[2]*bdscale);
  return Point(x(u),y(u),z(u));
  return a*u+b*(1.0-u);               
}