// done by Jarek Rossignac, with additions by Florian Hecht, Sept. 2006

//**************************************
// **** POINTS AND VECTOR CLASSES
class pt { float x,y; 
  pt (float px, float py) {x = px; y = py;};
  pt make() {return(new pt(x,y));};
  void show(int r) { ellipse(x, y, r, r); }; 
  void showLineTo (pt P) {line(x,y,P.x,P.y); }; 
  void setTo(pt P) { x = P.x; y = P.y; }; 
  void setToMouse() { x = mouseX; y = mouseY; if (!(x>=-width)&&(x<=2*width)) {println("***setToMouse x="+x);};}; 
  void write() {println("("+x+","+y+")");};
  void addVec(vec V) {x += V.x; y += V.y;};
  void addScaledVec(float s, vec V) {x += s*V.x; y += s*V.y;};
  void subVec(vec V) {x -= V.x; y -= V.y;};
  void vert() {vertex(x,y);};
  boolean isInWindow() {return(((x<0)||(x>600)||(y<0)||(y>600)));};
  void label(String s, vec D) {text(s, x+D.x,y+D.y);  };
  vec vecTo(pt P) {return(new vec(P.x-x,P.y-y)); };
  float disTo(pt P) {return(sqrt(sq(P.x-x)+sq(P.y-y))); };
  vec vecToMid(pt P, pt Q) {return(new vec((P.x+Q.x)/2.0-x,(P.y+Q.y)/2.0-y)); };
  vec vecToProp (pt B, pt D) {
      vec CB = this.vecTo(B); float LCB = CB.norm();
      vec CD = this.vecTo(D); float LCD = CD.norm();
      vec U = CB.make();
      vec V = CD.make(); V.sub(U); V.mul(LCB/(LCB+LCD));
      U.add(V);
      return(U);  
      };  
  void approach(pt P, float s) {x+=s*(P.x-x); y+=s*(P.y-y);};
  void interpolate(pt P, float s, pt Q) {x=P.x+s*(Q.x-P.x); y=P.y+s*(Q.y-P.y);};
  void addPt(pt P) {x+=P.x; y+=P.y;};
  void subPt(pt P) {x-=P.x; y-=P.y;};
  void mul(float f) {x*=f; y*=f;};
  void pers(float d) { y=d*y/(d+x); x=d*x/(d+x); };
  void pers2(float d) { y=d*y/(d+x); };
  void turnAround(pt P, float a) {float dx=x-P.x, dy=y-P.y, c=cos(a), s=sin(a); x=P.x+c*dx+s*dy; y=P.y-s*dx+c*dy; };
  pt movedBy(vec V) {return(new pt(x+V.x,y+V.y));};
  
  // new functions
  void setTo(float px, float py) {
    x = px; y = py;
  }
  
  vec toVec() {
    return new vec(x, y);
  }
} // end class pt

pt linter(pt a, float t, pt b) {
  pt p = new pt((1.0 - t)*a.x + t*b.x, (1.0 - t)*a.y + t*b.y);
  return p;
}

pt bezier(pt P1, pt P2, pt P3, pt P4, float s) {
  pt t1 = linter(P1, s, P2);
  pt t2 = linter(P2, s, P3);
  pt t3 = linter(P3, s, P4);
    
  pt t4 = linter(t1, s, t2);
  pt t5 = linter(t2, s, t3);
     
  return linter(t4, s, t5);
}

vec I=new vec (1,0);      // basis vector for angle calculations

float bezierAngle(pt P1, pt P2, pt P3, pt P4, float s) {
  float a;
  
  pt t1 = linter(P1, s, P2);
  pt t2 = linter(P2, s, P3);
  pt t3 = linter(P3, s, P4);
    
  pt t4 = linter(t1, s, t2);
  pt t5 = linter(t2, s, t3);
  
  vec tangent = t4.vecTo(t5);
  tangent.unit();
     
  if (tangent.y < 0)
    a = -acos(dot(I, tangent));
  else
    a = acos(dot(I, tangent));
    
  return a;
}

class vec { float x,y; 
  vec (float px, float py) {x = px; y = py;};
  vec make() {return(new vec(x,y));};
  void show (pt P) { ellipse(P.x, P.y, 3, 3); line(P.x,P.y,P.x+x,P.y+y); }; 
  void add(vec V) {x += V.x; y += V.y;};
  void addScaled(float m, vec V) {x += m*V.x; y += m*V.y;};
  void sub(vec V) {x -= V.x; y -= V.y;};
  void mul(float m) {x *= m; y *= m;};
  void div(float m) {x /= m; y /= m;};
  void write() {println("("+x+","+y+")");};
  float norm() {return(sqrt(sq(x)+sq(y)));}; 
  void unit() {float n=sqrt(sq(x)+sq(y)); if (n>0.0001) {x/=n; y/=n;};};
  vec left() {return(new vec(-y,x));};
  void back() {x= -x; y= -y;};
  vec turn(float a) {return(new vec(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a))); };
  void setTo(float px, float py) {x = px; y = py;}
  } 

pt mouse() {float x=mouseX, y=mouseY; x=max(x,0); y=max(y,0); x=min(x,width); y=min(y,height-55);  return(new pt(x,y));}; 
pt midPt(pt A, pt B) {return(new pt((A.x+B.x)/2,(A.y+B.y)/2)); };
vec midVec(vec U, vec V) {return(new vec((U.x+V.x)/2,(U.y+V.y)/2)); };
float dot(vec U, vec V) {return(U.x*V.x+U.y*V.y); };
float angle(vec U, vec V) {return(atan2(dot(U.left(),V),dot(U,V))); };

 void drawArc(pt A, pt B, pt C, int rec) {
   if (rec==0) { B.showLineTo(A); B.showLineTo(C); } 
   else { 
     float w = (A.vecTo(B).norm()+C.vecTo(B).norm())/2;
     float l = A.vecTo(C).norm()/2;
     float t = l/(w+l);
     pt L = A.make(); L.addScaledVec(t,A.vecTo(B)); 
     pt R = C.make(); R.addScaledVec(t,C.vecTo(B)); 
     pt M = midPt(L,R);
     drawArc(A,L, M,rec-1); drawArc(M,R, C,rec-1); 
     };
   };