//*************************************************************
//**** 2D GEOMETRY CLASSES AND UTILITIES, Jarek Rossignac *****
//****            REVISED October 3, 2008                 *****
//*************************************************************

// Transforms the local CS
void push2D() {pushMatrix();}                           // saves current LCS (local coordinate system) by pushing matrix stack
void tra2D(vec2D V) {translate(V.x,V.y);}                 // translates view by vector V
void tra2D(pt2D P) {translate(P.x,P.y);}                  // translates view from origin to P
void traRev2D(pt2D P) {translate(-P.x,-P.y);}             // translates view back by vector from P to origin
void rot2D(float a) {rotate(a);}                        // rotates view around origin by angle a
void rot2D(float a, pt2D P) {tra2D(P); rot2D(a); traRev2D(P);}  // rotates view by a around point P
void sca2D(float s) {scale(s);}                         // scales view wrt origin
void sca2D(float s, pt2D P) {tra2D(P); sca2D(s); traRev2D(P);}  // scales view wrt point P
void pop2D() {popMatrix();}                             // restores previously saved CS by popping stack
// screen and mouse
pt2D screenCenter2D() {return(new pt2D(height/2,height/2));}                                           // returns point at center of screen
boolean mouseIsInWindow2D() {return(((mouseX>0)&&(mouseX<height)&&(mouseY>0)&&(mouseY<height)));}; // returns true if mouse is in window
pt2D Mouse2D() {return(new pt2D(mouseX,mouseY));};                                                     // returns point with current mouse location
pt2D pMouse2D() {return(new pt2D(pmouseX,pmouseY));};                                                  // returns point with previous mouse location
vec2D MouseDrag2D() {return new vec2D(mouseX-pmouseX,mouseY-pmouseY);};                                // vector representing recent mouse displacement
float angleTracedWithMouse2D () {vec2D V1=U2D(screenCenter2D(),pMouse2D()); vec2D V2=U2D(screenCenter2D(),Mouse2D()); return a2D(V1,V2);}
float scaleTracedWithMouse2D () {float d1=d2D(screenCenter2D(),pMouse2D()); float d2=d2D(screenCenter2D(),Mouse2D()); return d2/d1;}
// Make points
pt2D P2D(float x, float y) {return new pt2D(x,y); };                                                       // make point (x,y)
pt2D P2D() {return P2D(0,0); };                                                                            // make point (0,0)
pt2D P2D(pt2D P) {return P2D(P.x,P.y); };                                                                    // make copy of point P
// Weighted combinations of points
pt2D A2D(pt2D A, pt2D B) {return P2D((A.x+B.x)/2.0,(A.y+B.y)/2.0); };                                          // Average: (A+B)/2
pt2D A2D(pt2D A, pt2D B, pt2D C) {return P2D((A.x+B.x+C.x)/3.0,(A.y+B.y+C.y)/3.0); };                            // Average: (A+B+C)/3
pt2D L2D(pt2D A, float s, pt2D B) {return P2D(A.x+s*(B.x-A.x),A.y+s*(B.y-A.y)); };                             // Linear interpolation: A+sAB
pt2D S2D(float a, pt2D A, float b, pt2D B){return P2D(a*A.x+b*B.x,a*A.y+b*B.y); }                                                         // Weighted sum: aA+bB, make sure that a+b=1
pt2D S2D(float a, pt2D A, float b, pt2D B, float c, pt2D C){return P2D(a*A.x+b*B.x+c*C.x,a*A.y+b*B.y+c*C.y); }                              // Weighted sum: aA+bB+cC, make sure that a+b+c=1
pt2D S2D(float a, pt2D A, float b, pt2D B, float c, pt2D C, float d, pt2D D){return P2D(a*A.x+b*B.x+c*C.x+d*D.x,a*A.y+b*B.y+c*C.y+d*D.y); }   // Weighted sum: aA+bB+cC+dD, make sure that a+b+c+d=1
// transformed points
pt2D T2D(pt2D P, vec2D V) {return P2D(P.x + V.x, P.y + V.y); }                                                 // P+V (transaltes P by displacement vector V
pt2D T2D(pt2D P, float s, vec2D V) {return T2D(P,S2D(s,V)); }                                                    // P+sV (transaltes P by a fraction s of displacement vector V
pt2D T2D(pt2D P, float s, pt2D Q) {return T2D(P,s,V2D(P,Q)); }                                                   // P+sPQ/||PQ|| transaltes P by a fraction s of PQ (towards Q)
pt2D T2D(pt2D O, float x, vec2D I, float y, vec2D J) {return P2D(O.x + x*I.x+ y*J.x, O.y + x*I.y+ y*J.y); }      // O+xI+yJ  (change of coordinate systems)
pt2D R2D(pt2D Q, float a) {float c=cos(a), s=sin(a); return P2D(c*Q.x-s*Q.y, s*Q.x+c*Q.y); }                 // Rotated Q by a around origin
pt2D R2D(pt2D Q, float a, pt2D P) {return T2D(P,R2D(V2D(P,Q),a)); }                                                // Rotated Q by a around P
// distance
float d2D(pt2D P, pt2D Q) {return n2D(V2D(P,Q));};                                                // ||AB||
float d22D(pt2D P, pt2D Q) {return n22D(V2D(P,Q));};                                              // ||AB||^2
// Make vectors
vec2D V2D(float x, float y) {return new vec2D(x,y); };                                     // make V(x,y)
vec2D V2D() {return V2D(0,0);}                                                             // V(0,0)
vec2D V2D(vec2D V) {return new vec2D(V.x,V.y); };                                            // make copy of V
vec2D V2D(pt2D P, pt2D Q) {return new vec2D(Q.x-P.x,Q.y-P.y);};                                // PQ (vector from P to Q)
// scale, and combine vectors
vec2D S2D(float s,vec2D V) {return V2D(s*V.x,s*V.y);};                                                                                        // sV
vec2D S2D(vec2D U, vec2D V) {return V2D(U.x+V.x,U.y+V.y); };                                                                                    // U+V
vec2D S2D(float a, vec2D A, float b, vec2D B){return V2D(a*A.x+b*B.x,a*A.y+b*B.y); }                                                            // Weighted sum: aA+bB
vec2D S2D(float a, vec2D A, float b, vec2D B, float c, vec2D C){return V2D(a*A.x+b*B.x+c*C.x,a*A.y+b*B.y+c*C.y); }                                // Weighted sum: aA+bB+cC
vec2D S2D(float a, vec2D A, float b, vec2D B, float c, vec2D C, float d, vec2D D){return V2D(a*A.x+b*B.x+c*C.x+d*D.x,a*A.y+b*B.y+c*C.y+d*D.y); }    // Weighted sum: aA+bB+cC+dD
vec2D L2D(vec2D U,float s,vec2D V) {return V2D(U.x+s*(V.x-U.x),U.y+s*(V.y-U.y));};                                                              // (1-s)U+sV
vec2D A2D(vec2D U, vec2D V) {return L2D(U,0.5,V); };                                                                                            // (U+V)/2
// dot product, norm, normalization
float dot2D(vec2D U, vec2D V) {return U.x*V.x+U.y*V.y; };                                     // U*V  (dot product)
float n22D(vec2D V) {return sq(V.x)+sq(V.y);};                                              // V*V  (V^2)
float n2D(vec2D V) {return sqrt(n22D(V));};                                                   // ||V||
vec2D U2D(vec2D V) {float n = n2D(V); if (n==0) return null; else return S2D(1/n,V);};            // V/||V||  (normalized V, unit vector)
vec2D U2D(pt2D Q, pt2D P) {return U2D(V2D(Q,P)); };                                                 // Unit PQ
// angles and rotations
float a2D (vec2D U, vec2D V) {float a = atan2(dot2D(R2D(U),V),dot2D(U,V)); if(a>PI) return mPItoPIangle2D(a-2*PI); if(a<-PI) return mPItoPIangle2D(a+2*PI); return a; };  // angle(U,V) between -PI and PI
vec2D R2D(vec2D V) {return V2D(-V.y,V.x);};                                                         // V rotated 90 degrees (clockwise as seen on screen)
vec2D R2D(pt2D Q, pt2D P) {return R2D(V2D(Q,P)); };                                                     // vector QP rotated 90 degrees
vec2D R2D(vec2D U, float a) {float c=cos(a), s=sin(a); return V2D(c*U.x-s*U.y, s*U.x+c*U.y); };     // U rotated by a
vec2D R2D(vec2D U, float t, vec2D V) {return S2D(1+t*(n2D(V)-n2D(U))/n2D(U),R2D(U,t*angle2D(U,V)));  };         // interpolation (angle and length) between U and V
// three points
float a2D(pt2D A, pt2D B, pt2D C) {return  a2D(V2D(B,A),V2D(B,C)); }                                 // angle (BA,BC)
float turnAngle2D(pt2D A, pt2D B, pt2D C) {return  a2D(V2D(A,B),V2D(B,C)); }                         // angle (AB,BC)
boolean isRightTurn2D(pt2D A, pt2D B, pt2D C) {return dot2D( R2D(A,B),V2D(B,C)) > 0 ; };             // right turn (as seen on screen)
boolean isRightOf2D(pt2D A, pt2D Q, vec2D T) {return dot2D(R2D(T),V2D(Q,A)) > 0 ; };                 // A is on right of ray(Q,T) (as seen on screen)
boolean isInFrontOf2D(pt2D A, pt2D Q, vec2D T) {return dot2D(T,V2D(Q,A)) > 0 ; };                  // A is in frontof ray(Q,T) 
// display
void v2D(pt2D P) {vertex(P.x,P.y);};                                                     // next point when drawing polygons between beginShape(); and endShape();
void showCross2D(pt2D P, float r) {line(P.x-r,P.y,P.x+r,P.y); line(P.x,P.y-r,P.x,P.y+r);};   // shows P as cross of length r
void showCross2D(pt2D P) {showCross2D(P,2);};                                                      // shows P as small cross
void show2D(pt2D P, float r) {ellipse(P.x, P.y, 2*r, 2*r);};                             // draws circle of center r around point
void show2D(pt2D P) {ellipse(P.x, P.y, 4,4);};                                           // draws small circle around point
void show2D(pt2D P, pt2D Q) {line(P.x,P.y,Q.x,Q.y); };                                     // draws edge (P,Q)
void show2D(pt2D P, vec2D V) {line(P.x,P.y,P.x+V.x,P.y+V.y); }                             // show line from P along V
void show2D(pt2D P, float s, vec2D V) {show2D(P,S2D(s,V));}                                    // show line from P along sV
void arrow2D(pt2D P, pt2D Q) {arrow2D(P,V2D(P,Q)); }                                           // draws arrow from P to Q
void arrow2D(pt2D P, float s, vec2D V) {arrow2D(P,S2D(s,V));}                                   // show arrow from P along sV
void arrow2D(pt2D P, vec2D V) {show2D(P,V);  float n=n2D(V); float s=max(min(0.2,20./n),6./n);  // show arrow from P along V
  pt2D Q=T2D(P,V); vec2D U = S2D(-s,V); vec2D W = R2D(S2D(.3,U)); beginShape(); v2D(T2D(T2D(Q,U),W)); v2D(Q); v2D(T2D(T2D(Q,U),-1,W)); endShape(CLOSE);}; 

//************************************************************************
//**** ANGLES
//************************************************************************
float angle2D(vec2D V) {return(atan2(V.y,V.x)); };
float angle2D(vec2D U, vec2D V) {return(atan2(dot2D(R2D(U),V),dot2D(U,V))); };
float mPItoPIangle2D(float a) { if(a>PI) return(mPItoPIangle2D(a-2*PI)); if(a<-PI) return(mPItoPIangle2D(a+2*PI)); return(a);};
float toDeg2D(float a) {return(a*180/PI);}
float toRad2D(float a) {return(a*PI/180);}
//************************************************************************
//**** POINTS
//************************************************************************
class pt2D { float x=0,y=0; 
  // CREATE
  pt2D () {}
  pt2D (float px, float py) {x = px; y = py;};
  pt2D (pt2D P) {x = P.x; y = P.y;};
  pt2D (pt2D P, vec2D V) {x = P.x+V.x; y = P.y+V.y;};
  pt2D (pt2D P, float s, vec2D V) {x = P.x+s*V.x; y = P.y+s*V.y;};
  pt2D (pt2D A, float s, pt2D B) {x = A.x+s*(B.x-A.x); y = A.y+s*(B.y-A.y);};

  // MODIFY
  void setTo(float px, float py) {x = px; y = py;};  
  void setTo(pt2D P) {x = P.x; y = P.y;}; 
  void setToMouse() { x = mouseX; y = mouseY; }; 
  void trackMouse() { x += mouseX-pmouseX; y += mouseY-pmouseY; }; 
  void translateToTrack(float s, pt2D P) {setTo(T2D(P,s,V2D(P,this)));};
  void track(float s, pt2D P) {setTo(T2D(P,s,V2D(P,this)));};
  pt2D scaleBy(float f) {x*=f; y*=f; return this;};
  void scaleBy(float u, float v) {x*=u; y*=v;};
  void add(vec2D V) {x += V.x; y += V.y;};   
  void add(float s, vec2D V) {x += s*V.x; y += s*V.y;};   
  void translateBy(vec2D V) {x += V.x; y += V.y;};   
  void translateBy(float s, vec2D V) {x += s*V.x; y += s*V.y;};   
  void translateBy(float u, float v) {x += u; y += v;};
  void translateTowards(float s, pt2D P) {x+=s*(P.x-x);  y+=s*(P.y-y); };
  void translateTowardsBy(float s, pt2D P) {vec2D V = this.makeVecTo(P); V.normalize(); this.translateBy(s,V); };
  void addPt(pt2D P) {x += P.x; y += P.y;};        // incorrect notation, but useful for computing weighted averages
  void addScaledPt(float s, pt2D P) {x += s*P.x; y += s*P.y;};        // incorrect notation, but useful for computing weighted averages
  void rotateBy(float a) {float dx=x, dy=y, c=cos(a), s=sin(a); x=c*dx-s*dy; y=+s*dx+c*dy; };     // around origin
  void rotateBy(float a, pt2D P) {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; };   // around point P
  void rotateBy(float s, float t, pt2D P) {float dx=x-P.x, dy=y-P.y; dx-=dy*t; dy+=dx*s; dx-=dy*t; x=P.x+dx; y=P.y+dy; };   // s=sin(a); t=tan(a/2);
  void clipToWindow() {x=max(x,0); y=max(y,0); x=min(x,height); y=min(y,height); }
  
  // OUTPUT POINT
  pt2D clone() {return new pt2D(x,y); };
  pt2D makeClone() {return new pt2D(x,y); };
  pt2D makeTranslatedBy(vec2D V) {return(new pt2D(x + V.x, y + V.y));};
  pt2D makeTranslatedBy(float s, vec2D V) {return(new pt2D(x + s*V.x, y + s*V.y));};
  pt2D makeTransaltedTowards(float s, pt2D P) {return(new pt2D(x + s*(P.x-x), y + s*(P.y-y)));};
  pt2D makeTranslatedBy(float u, float v) {return(new pt2D(x + u, y + v));};
  pt2D makeRotatedBy(float a, pt2D P) {float dx=x-P.x, dy=y-P.y, c=cos(a), s=sin(a); return(new pt2D(P.x+c*dx-s*dy, P.y+s*dx+c*dy)); };
  pt2D makeRotatedBy(float a) {float dx=x, dy=y, c=cos(a), s=sin(a); return(new pt2D(c*dx-s*dy, s*dx+c*dy)); };
  pt2D makeProjectionOnLine(pt2D P, pt2D Q) {float a=dot2D(P.makeVecTo(this),P.makeVecTo(Q)), b=dot2D(P.makeVecTo(Q),P.makeVecTo(Q)); return(P.makeTransaltedTowards(a/b,Q)); };
  pt2D makeOffset(pt2D P, pt2D Q, float r) {
    float a = angle2D(vecTo(P),vecTo(Q))/2;
    float h = r/tan(a); 
    vec2D T = vecTo(P); T.normalize(); vec2D N = T.left();
    pt2D R = new pt2D(x,y); R.translateBy(h,T); R.translateBy(r,N);
    return R; };

   // OUTPUT VEC
  vec2D vecTo(pt2D P) {return(new vec2D(P.x-x,P.y-y)); };
  vec2D makeVecTo(pt2D P) {return(new vec2D(P.x-x,P.y-y)); };
  vec2D makeVecToCenter () {return(new vec2D(x-height/2.,y-height/2.)); };
  vec2D makeVecToAverage (pt2D P, pt2D Q) {return(new vec2D((P.x+Q.x)/2.0-x,(P.y+Q.y)/2.0-y)); };
  vec2D makeVecToAverage (pt2D P, pt2D Q, pt2D R) {return(new vec2D((P.x+Q.x+R.x)/3.0-x,(P.y+Q.y+R.x)/3.0-y)); };
  vec2D makeVecToMouse () {return(new vec2D(mouseX-x,mouseY-y)); };
  vec2D makeVecToBisectProjection (pt2D P, pt2D Q) {float a=disTo(P), b=disTo(Q);  return makeVecTo(L2D(P,a/(a+b),Q)); };
  vec2D makeVecToNormalProjection (pt2D P, pt2D Q) {float a=dot2D(V2D(P,this),V2D(P,Q)), b=d22D(P,Q); return V2D(this,L2D(P,a/b,Q)); };
 // vec makeVecTowards(pt P, float d) {vec V = makeVecTo(P); float n = V.norm(); V.normalize(); V.scaleBy(d-n); return V; };
 
  // OUTPUT TEST OR MEASURE
  float disTo(pt2D P) {return(sqrt(sq(P.x-x)+sq(P.y-y))); };
  float disToMouse() {return(sqrt(sq(x-mouseX)+sq(y-mouseY))); };
  boolean isInWindow() {return(((x>0)&&(x<height)&&(y>0)&&(y<height)));};
  boolean projectsBetween(pt2D P, pt2D Q) {float a=dot2D(P.makeVecTo(this),P.makeVecTo(Q)), b=dot2D(P.makeVecTo(Q),P.makeVecTo(Q)); return((0<a)&&(a<b)); };
  float ratioOfProjectionBetween(pt2D P, pt2D Q) {float a=dot2D(P.makeVecTo(this),P.makeVecTo(Q)), b=dot2D(P.makeVecTo(Q),P.makeVecTo(Q)); return(a/b); };
  float disToLine(pt2D P, pt2D Q) {float a=dot2D(P.makeVecTo(this),P.makeVecTo(Q).makeUnit().makeTurnedLeft()); return(abs(a)); };
  boolean isLeftOf(pt2D P, pt2D Q) {boolean l=dot2D(P.makeVecTo(this),P.makeVecTo(Q).makeTurnedLeft())>0; return(l);  };
  boolean isInTriangle(pt2D A, pt2D B, pt2D C) { boolean a = this.isLeftOf(B,C); boolean b = this.isLeftOf(C,A); boolean c = this.isLeftOf(A,B); return((a&&b&&c)||(!a&&!b&&!c));};
  boolean isInCircle(pt2D C, float r) {return d2D(this,C)<r; }  // returns true if point is in circle C,r
  float triArea(pt2D P, pt2D Q) {float a=dot2D(P.makeVecTo(this),P.makeVecTo(Q).makeTurnedLeft()); return(a/2.); };

  
  // DRAW , PRINT
  void show() {ellipse(x, y, height/200, height/200); }; // shows point as small dot
  void show(float r) {ellipse(x, y, 2*r, 2*r); }; // shows point as disk of radius r
  void showCross(float r) {line(x-r,y,x+r,y); line(x,y-r,x,y+r);}; 
  void v() {vertex(x,y);};  // used for drawing polygons between beginShape(); and endShape();
  void write() {println("("+x+","+y+")");};  // writes point coordinates in text window
  void showLabel(String s, vec2D D) {text(s, x+D.x,y+D.y);  };  // show string displaced by vector D from point
  void showLabel(String s) {text(s, x+5,y+4);  };
  void showLabel(int i) {text(str(i), x+5,y+4);  };  // shows integer number next to point
  void showLabel(String s, float u, float v) {text(s, x+u, y+v);  };
  void showSegmentTo (pt2D P) {line(x,y,P.x,P.y); }; // draws edge to another point
  void to (pt2D P) {line(x,y,P.x,P.y); }; // draws edge to another point

  } // end of pt class


//NEW CODE: Calculating barycentre of triangle (n-point) using code from triangles demo applet.
//for two points:
pt2D weightedSum(float a, pt2D A, float b, pt2D B) {pt2D P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B); return(P);}

//for three points:
pt2D weightedSum(float a, pt2D A, float b, pt2D B, float c, pt2D C) {pt2D P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B);  P.addScaledPt(c,C); return(P);}

//for four points:
pt2D weightedSum(float a, pt2D A, float b, pt2D B, float c, pt2D C, float d, pt2D D) 
        {pt2D P =  A.makeClone(); P.scaleBy(a); P.addScaledPt(b,B); P.addScaledPt(c,C); P.addScaledPt(d,D); return(P);}




//************************************************************************
//**** VECTORS
//************************************************************************
class vec2D { float x=0,y=0; 
 // CREATE
  vec2D () {};
  vec2D (vec2D V) {x = V.x; y = V.y;};
  vec2D (float s, vec2D V) {x = s*V.x; y = s*V.y;};
  vec2D (float px, float py) {x = px; y = py;};
  vec2D (pt2D P, pt2D Q) {x = Q.x-P.x; y = Q.y-P.y;};
 
 // MODIFY
  void setTo(float px, float py) {x = px; y = py;}; 
  void setTo(pt2D P, pt2D Q) {x = Q.x-P.x; y = Q.y-P.y;}; 
  void setTo(vec2D V) {x = V.x; y = V.y;}; 
  void scaleBy(float f) {x*=f; y*=f;};
  void back() {x=-x; y=-y;};
  void mul(float f) {x*=f; y*=f;};
  void div(float f) {x/=f; y/=f;};
  void scaleBy(float u, float v) {x*=u; y*=v;};
  void normalize() {float n=sqrt(sq(x)+sq(y)); if (n>0.000001) {x/=n; y/=n;};};
  void add(vec2D V) {x += V.x; y += V.y;};   
  void add(float s, vec2D V) {x += s*V.x; y += s*V.y;};   
  void add(float u, float v) {x += u; y += v;};
  void turnLeft() {float w=x; x=-y; y=w;};
  void rotateBy (float a) {float xx=x, yy=y; x=xx*cos(a)-yy*sin(a); y=xx*sin(a)+yy*cos(a); };
  
  // OUTPUT VEC
  vec2D makeClone() {return(new vec2D(x,y));}; 
  vec2D makeUnit() {float n=sqrt(sq(x)+sq(y)); if (n<0.000001) n=1; return(new vec2D(x/n,y/n));}; 
  vec2D unit() {float n=sqrt(sq(x)+sq(y)); if (n<0.000001) n=1; return(new vec2D(x/n,y/n));}; 
  vec2D makeScaledBy(float s) {return(new vec2D(x*s, y*s));};
  vec2D makeTurnedLeft() {return(new vec2D(-y,x));};
  vec2D left() {return(new vec2D(-y,x));};
  vec2D makeOffsetVec(vec2D V) {return(new vec2D(x + V.x, y + V.y));};
  vec2D makeOffsetVec(float s, vec2D V) {return(new vec2D(x + s*V.x, y + s*V.y));};
  vec2D makeOffsetVec(float u, float v) {return(new vec2D(x + u, y + v));};
  vec2D makeRotatedBy(float a) {return(new vec2D(x*cos(a)-y*sin(a),x*sin(a)+y*cos(a))); };
  vec2D makeReflectedVec(vec2D N) { return makeOffsetVec(-2.*dot2D(this,N),N);};

  // OUTPUT TEST MEASURE
  float norm() {return(sqrt(sq(x)+sq(y)));}
  boolean isNull() {return((abs(x)+abs(y)<0.000001));}
  float angle() {return(atan2(y,x)); }

  // DRAW, PRINT
  void write() {println("("+x+","+y+")");};
  void show (pt2D P) {line(P.x,P.y,P.x+x,P.y+y); }; 
  void showAt (pt2D P) {line(P.x,P.y,P.x+x,P.y+y); }; 
  void showArrowAt (pt2D P) {line(P.x,P.y,P.x+x,P.y+y); 
      float n=min(this.norm()/10.,height/50.); 
      pt2D Q=P.makeTranslatedBy(this); 
      vec2D U = this.makeUnit().makeScaledBy(-n);
      vec2D W = U.makeTurnedLeft().makeScaledBy(0.3);
      beginShape(); Q.makeTranslatedBy(U).makeTranslatedBy(W).v(); Q.v(); 
                    W.scaleBy(-1); Q.makeTranslatedBy(U).makeTranslatedBy(W).v(); endShape(CLOSE); }; 
  void showLabel(String s, pt2D P) {pt2D Q = P.makeTranslatedBy(0.5,this); 
           vec2D N = makeUnit().makeTurnedLeft(); Q.makeTranslatedBy(3,N).showLabel(s); };
  } // end vec class
 
 
 // fitting
 vec2D proVec2D (pt2D B, pt2D C, pt2D D) { return L2D(V2D(C,B), d2D(C,B)/(d2D(C,B)+d2D(C,D)),V2D(C,D)); }
 vec2D vecToCubic2D (pt2D A, pt2D B, pt2D C, pt2D D, pt2D E) {return V2D( (-A.x+4*B.x-6*C.x+4*D.x-E.x)/6, (-A.y+4*B.y-6*C.y+4*D.y-E.y)/6  );}
// CURVES
void drawCubicBezier2D(pt2D A, pt2D B, pt2D C, pt2D D) {bezier(A.x,A.y,B.x,B.y,C.x,C.y,D.x,D.y);}
pt2D B2D(pt2D A, pt2D B, pt2D C, float s) {return( L2D(L2D(B,s/4.,A),0.5,L2D(B,s/4.,C))); };                            // returns a tucked B towards its neighbors
pt2D F2D(pt2D A, pt2D B, pt2D C, pt2D D, float s) {return( L2D(L2D(A,1.+(1.-s)/8.,B) ,0.5, L2D(D,1.+(1.-s)/8.,C))); };    // returns a bulged mid-edge point 
pt2D limit2D(pt2D A, pt2D B, pt2D C, pt2D D, pt2D E, float s, int r) {
  if (r==0) return C.clone();
  else return limit2D(B2D(A,B,C,s),F2D(A,B,C,D,s),B2D(B,C,D,s),F2D(B,C,D,E,s),B2D(C,D,E,s),s,r-1);
  }
  
boolean edgesCross(pt2D A, pt2D B, pt2D C, pt2D D) {return (dot2D(U2D(R2D(V2D(A,B))),V2D(A,C))>0 != dot2D(U2D(R2D(V2D(A,B))),V2D(A,D))>0) // tests whether two edges cross
                                                && (dot2D(U2D(R2D(V2D(C,D))),V2D(C,A))>0 != dot2D(U2D(R2D(V2D(C,D))),V2D(C,B))>0) ; }

  
//************************************************************************
//**** CIRCLES
//************************************************************************
pt circumCenter (pt A3, pt B3, pt C3) {    // computes the center of a circumscirbing circle to triangle (A,B,C)
  pt2D A = new pt2D(A3.x,A3.y);
  pt2D B = new pt2D(B3.x,B3.y);
  pt2D C = new pt2D(C3.x,C3.y);
  vec2D AB =  A.makeVecTo(B);  float ab2 = dot2D(AB,AB);
  vec2D AC =  A.makeVecTo(C); AC.turnLeft();  float ac2 = dot2D(AC,AC);
  float d = 2*dot2D(AB,AC);
  AB.turnLeft();
  AB.scaleBy(-ac2); 
  AC.scaleBy(ab2);
  AB.add(AC);
  AB.scaleBy(1./d);
  pt2D X =  A.makeClone();
  X.translateBy(AB);
  return(new pt(X.x, X.y, 0));//returns a 3D point, but ignores the Z
  };

boolean ccw(pt A3, pt B3, pt C3) {
  pt2D A = new pt2D(A3.x,A3.y);
  pt2D B = new pt2D(B3.x,B3.y);
  pt2D C = new pt2D(C3.x,C3.y);
  return C.isLeftOf(A,B) ;
}