// geometry.asy

// Author: Philippe IVALDI 2007/09/01
// http://www.piprime.fr/

// This program is free software ; you can redistribute it and/or modify
// the Free Software Foundation ; either version 3 of the License, or
// (at your option) any later version.

// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY ; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License
// along with this program ; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

// COMMENTARY:
// An Asymptote geometry module.

// THANKS:
// Special thanks to Olivier Guibé for his help in mathematical issues.

// BUGS:

// CODE:

import math;
import markers;
// *=======================================================*

real epsgeo = 10 * sqrt(realEpsilon);/*Variable used in the approximate calculations.View index ordered by [name] - [type]*/

real lmargin = 0, real bmargin = 0,
real rmargin = lmargin, real tmargin = bmargin,
bool rigid = true, bool allObject = true)
{/*Add margins to 'pic' with respect to
the current bounding box of 'pic'.
If 'rigid' is false, margins are added iff an infinite curve will
be prolonged on the margin.
If 'allObject' is false, fixed - size objects (such as labels and
arrowheads) will be ignored.View index ordered by [name] - [type]*/
pair m = allObject ? truepoint(pic, SW) : point(pic, SW);
pair M = allObject ? truepoint(pic, NE) : point(pic, NE);
if(rigid) {
draw(m - inverse(pic.calculateTransform()) * (lmargin, bmargin), invisible);
draw(M + inverse(pic.calculateTransform()) * (rmargin, tmargin), invisible);
} else pic.addBox(m, M, -(lmargin, bmargin), (rmargin, tmargin));
}

real approximate(real t)
{
real ot = t;
if(abs(t - ceil(t)) < epsgeo) ot = ceil(t);
else if(abs(t - floor(t)) < epsgeo) ot = floor(t);
return ot;
}

real[] approximate(real[] T)
{
return map(approximate, T);
}

real binomial(real n, real k)
{/*Return n!/((n - k)!*k!)View index ordered by [name] - [type]*/
return gamma(n + 1)/(gamma(n - k + 1) * gamma(k + 1));
}

real rf(real x, real y, real z)
{/*Computes Carlson's elliptic integral of the first kind.
x, y, and z must be non negative, and at most one can be zero.View index ordered by [name] - [type]*/
real ERRTOL = 0.0025,
TINY = 1.5e-38,
BIG = 3e37,
THIRD = 1/3,
C1 = 1/24,
C2 = 0.1,
C3 = 3/44,
C4 = 1/14;
real alamb, ave, delx, dely, delz, e2, e3, sqrtx, sqrty, sqrtz, xt, yt, zt;
if(min(x, y, z) < 0 || min(x + y, x + z, y + z) < TINY ||
max(x, y, z) > BIG) abort("rf: invalid arguments.");
xt = x;
yt = y;
zt = z;
do {
sqrtx = sqrt(xt);
sqrty = sqrt(yt);
sqrtz = sqrt(zt);
alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz;
xt = 0.25 * (xt + alamb);
yt = 0.25 * (yt + alamb);
zt = 0.25 * (zt + alamb);
ave = THIRD * (xt + yt + zt);
delx = (ave - xt)/ave;
dely = (ave - yt)/ave;
delz = (ave - zt)/ave;
} while(max(fabs(delx), fabs(dely), fabs(delz)) > ERRTOL);
e2 = delx * dely - delz * delz;
e3 = delx * dely * delz;
return (1.0 + (C1 * e2 - C2 - C3 * e3) * e2 + C4 * e3)/sqrt(ave);
}

real rd(real x, real y, real z)
{/*Computes Carlson's elliptic integral of the second kind.
x and y must be positive, and at most one can be zero.
z must be non negative.View index ordered by [name] - [type]*/
real ERRTOL = 0.0015,
TINY = 1e-25,
BIG = 4.5 * 10.0^21,
C1 = (3/14),
C2 = (1/6),
C3 = (9/22),
C4 = (3/26),
C5 = (0.25 * C3),
C6 = (1.5 * C4);
real alamb, ave, delx, dely, delz, ea, eb, ec, ed, ee, fac, sqrtx, sqrty,
sqrtz, sum, xt, yt, zt;
if (min(x, y) < 0 || min(x + y, z) < TINY || max(x, y, z) > BIG)
abort("rd: invalid arguments");
xt = x;
yt = y;
zt = z;
sum = 0;
fac = 1;
do {
sqrtx = sqrt(xt);
sqrty = sqrt(yt);
sqrtz = sqrt(zt);
alamb = sqrtx * (sqrty + sqrtz) + sqrty * sqrtz;
sum += fac/(sqrtz * (zt + alamb));
fac = 0.25 * fac;
xt = 0.25 * (xt + alamb);
yt = 0.25 * (yt + alamb);
zt = 0.25 * (zt + alamb);
ave = 0.2 * (xt + yt + 3.0 * zt);
delx = (ave - xt)/ave;
dely = (ave - yt)/ave;
delz = (ave - zt)/ave;
} while (max(fabs(delx), fabs(dely), fabs(delz)) > ERRTOL);
ea = delx * dely;
eb = delz * delz;
ec = ea - eb;
ed = ea - 6 * eb;
ee = ed + ec + ec;
return 3 * sum + fac * (1.0 + ed * (-C1 + C5 * ed - C6 * delz * ee)
+delz * (C2 * ee + delz * (-C3 * ec + delz * C4 * ea)))/(ave * sqrt(ave));
}

real elle(real phi, real k)
{/*Legendre elliptic integral of the 2nd kind,
evaluated using Carlson's functions RD and RF.
The argument ranges are -infinity < phi < +infinity, 0 <= k * sin(phi) <= 1.View index ordered by [name] - [type]*/
real result;
if (phi >= 0 && phi <= pi/2) {
real cc, q, s;
s = sin(phi);
cc = cos(phi)^2;
q = (1 - s * k) * (1 + s * k);
result = s * (rf(cc, q, 1) - (s * k)^2 * rd(cc, q, 1)/3);
} else
if (phi <= pi && phi >= 0) {
result = 2 * elle(pi/2, k) - elle(pi - phi, k);
} else
if (phi <= 3 * pi/2 && phi >= 0) {
result = 2 * elle(pi/2, k) + elle(phi - pi, k);
} else
if (phi <= 2 * pi && phi >= 0) {
result = 4 * elle(pi/2, k) - elle(2 * pi - phi, k);
} else
if (phi >= 0) {
int nb = floor(0.5 * phi/pi);
result = nb * elle(2 * pi, k) + elle(phi%(2 * pi), k);
} else result = -elle(-phi, k);
return result;
}

pair[] intersectionpoints(pair A, pair B,
real a, real b, real c, real d, real f, real g)
{/*Intersection points with the line (AB) and the quadric curve
a * x^2 + b * x * y + c * y^2 + d * x + f * y + g = 0 given in the default coordinate systemView index ordered by [name] - [type]*/
pair[] op;
real ap = B.y - A.y,
bpp = A.x - B.x,
cp = A.y * B.x - A.x * B.y;
real sol[];
if (abs(ap) > epsgeo) {
real aa = ap * c + a * bpp^2/ap - b * bpp,
bb = ap * f - bpp * d + 2 * a * bpp * cp/ap - b * cp,
cc = ap * g - cp * d + a * cp^2/ap;
for (int i = 0; i < sol.length; ++i) {
op.push((-bpp * sol[i]/ap - cp/ap, sol[i]));
}
} else {
real aa = a * bpp,
bb = d * bpp - b * cp,
cc = g * bpp - cp * f + c * cp^2/bpp;
for (int i = 0; i < sol.length; ++i) {
op.push((sol[i], -cp/bpp));
}
}
return op;
}

pair[] intersectionpoints(pair A, pair B, real[] equation)
{/*Return the intersection points of the line AB with
the conic whose an equation is
equation[0] * x^2 + equation[1] * x * y + equation[2] * y^2 + equation[3] * x + equation[4] * y + equation[5] = 0View index ordered by [name] - [type]*/
if(equation.length != 6) abort("intersectionpoints: bad length of array for a conic equation.");
return intersectionpoints(A, B, equation[0], equation[1], equation[2],
equation[3], equation[4], equation[5]);
}
// *=======================================================*

// *=======================================================*
// *......................COORDINATES......................*

real EPS = sqrt(realEpsilon);

typedef pair convert(pair);/*Function type to convert pair in an other coordinate system.*/

typedef real abs(pair);/*Function type to calculate modulus of pair.*/

typedef real dot(pair, pair);/*Function type to calculate dot product.*/

typedef pair polar(real, real);/*Function type to calculate the coordinates from the polar coordinates.*/

struct coordsys
{/*This structure represents a coordinate system in the plane.*/

restricted convert relativetodefault = new pair(pair m){return m;};/*Convert a pair given relatively to this coordinate system to
the pair relatively to the default coordinate system.*/

restricted convert defaulttorelative = new pair(pair m){return m;};/*Convert a pair given relatively to the default coordinate system to
the pair relatively to this coordinate system.*/

restricted dot dot = new real(pair m, pair n){return dot(m, n);};/*Return the dot product of this coordinate system.*/

restricted abs abs = new real(pair m){return abs(m);};/*Return the modulus of a pair in this coordinate system.*/

restricted polar polar = new pair(real r, real a){return (r * cos(a), r * sin(a));};/*Polar coordinates routine of this coordinate system.*/

restricted pair O = (0, 0), i = (1, 0), j = (0, 1);/*Origin and units vector.*/

void init(convert rtd, convert dtr,
polar polar, dot dot)
{/*The default constructor of the coordinate system.*/
this.relativetodefault = rtd;
this.defaulttorelative = dtr;
this.polar = polar;
this.dot = dot;
this.abs = new real(pair m){return sqrt(dot(m, m));};;
this.O = rtd((0, 0));
this.i = rtd((1, 0)) - O;
this.j = rtd((0, 1)) - O;
}
}/*View index ordered by [name] - [type]*/

bool operator ==(coordsys c1, coordsys c2)
{/*Return true iff the coordinate system have the same origin and units vector.View index ordered by [name] - [type]*/
return c1.O == c2.O && c1.i == c2.i && c1.j == c2.j;
}

coordsys cartesiansystem(pair O = (0, 0), pair i, pair j)
{/*Return the Cartesian coordinate system (O, i, j).View index ordered by [name] - [type]*/
coordsys R;
real[][] P = {{0, 0}, {0, 0}};
real[][] iP;
P[0][0] = i.x;
P[0][1] = j.x;
P[1][0] = i.y;
P[1][1] = j.y;
iP = inverse(P);
real ni = abs(i);
real nj = abs(j);
real ij = angle(j) - angle(i);

pair rtd(pair m)
{
return O + (P[0][0] * m.x + P[0][1] * m.y, P[1][0] * m.x + P[1][1] * m.y);
}

pair dtr(pair m)
{
m-=O;
return (iP[0][0] * m.x + iP[0][1] * m.y, iP[1][0] * m.x + iP[1][1] * m.y);
}

pair polar(real r, real a)
{
real ca = sin(ij - a)/(ni * sin(ij));
real sa = sin(a)/(nj * sin(ij));
return r * (ca, sa);
}

real tdot(pair m, pair n)
{
return m.x * n.x * ni^2 + m.y * n.y * nj^2 + (m.x * n.y + n.x * m.y) * dot(i, j);
}

R.init(rtd, dtr, polar, tdot);
return R;
}

void show(picture pic = currentpicture, Label lo = "$O$",
Label li = "$\vec{\imath}$",
Label lj = "$\vec{\jmath}$",
coordsys R,
pen dotpen = currentpen, pen xpen = currentpen, pen ypen = xpen,
pen ipen = red,
pen jpen = ipen,
arrowbar arrow = Arrow)
{/*Draw the components (O, i, j, x - axis, y - axis) of 'R'.View index ordered by [name] - [type]*/
unravel R;
dot(pic, O, dotpen);
drawline(pic, O, O + i, xpen);
drawline(pic, O, O + j, ypen);
draw(pic, li, O--(O + i), ipen, arrow);
Label lj = lj.copy();
lj.align(lj.align, unit(I * j));
draw(pic, lj, O--(O + j), jpen, arrow);
draw(pic, lj, O--(O + j), jpen, arrow);
Label lo = lo.copy();
lo.align(lo.align, -2 * dir(O--O + i, O--O + j));
lo.p(dotpen);
label(pic, lo, O);
}

pair operator /(pair p, coordsys R)
{/*Return the xy - coordinates of 'p' relatively to
the coordinate system 'R'.
For example, if R = cartesiansystem((1, 2), (1, 0), (0, 1)), (0, 0)/R is (-1, -2).View index ordered by [name] - [type]*/
return R.defaulttorelative(p);
}

pair operator *(coordsys R, pair p)
{/*Return the coordinates of 'p' given in the
xy - coordinates 'R'.
For example, if R = cartesiansystem((1, 2), (1, 0), (0, 1)), R * (0, 0) is (1, 2).View index ordered by [name] - [type]*/
return R.relativetodefault(p);
}

path operator *(coordsys R, path g)
{/*Return the reconstructed path applying R * pair to each node, pre and post control point of 'g'.View index ordered by [name] - [type]*/
guide og = R * point(g, 0);
real l = length(g);
for(int i = 1; i <= l; ++i)
{
pair P = R * point(g, i);
pair post = R * postcontrol(g, i - 1);
pair pre = R * precontrol(g, i);
if(i == l && (cyclic(g)))
og = og..controls post and pre..cycle;
else
og = og..controls post and pre..P;
}
return og;
}

coordsys operator *(transform t,coordsys R)
{/*Provide transform * coordsys.
Note that shiftless(t) is applied to R.i and R.j.View index ordered by [name] - [type]*/
coordsys oc;
oc = cartesiansystem(t * R.O, shiftless(t) * R.i, shiftless(t) * R.j);
return oc;
}

restricted coordsys defaultcoordsys = cartesiansystem(0, (1, 0), (0, 1));/*One can always refer to the default coordinate system using this constant.View index ordered by [name] - [type]*/

coordsys currentcoordsys = defaultcoordsys;/*The coordinate system used by default.View index ordered by [name] - [type]*/

struct point
{/*This structure replaces the pair to embed its coordinate system.
For example, if 'P = point(cartesiansystem((1, 2), i, j), (0, 0))',
P is equal to the pair (1, 2).*/

coordsys coordsys;/*The coordinate system of this point.*/
restricted pair coordinates;/*The coordinates of this point relatively to the coordinate system 'coordsys'.*/
restricted real x, y;/*The xpart and the ypart of 'coordinates'.*/
/**/
real m = 1;/*Used to cast mass<->point.*/
void init(coordsys R, pair coordinates, real mass)
{/*The constructor.*/
this.coordsys = R;
this.coordinates = coordinates;
this.x = coordinates.x;
this.y = coordinates.y;
this.m = mass;
}
}/*View index ordered by [name] - [type]*/

point point(coordsys R, pair p, real m = 1)
{/*Return the point which has the coodinates 'p' in the
coordinate system 'R' and the mass 'm'.View index ordered by [name] - [type]*/
point op;
op.init(R, p, m);
return op;
}

point point(explicit pair p, real m)
{/*Return the point which has the coodinates 'p' in the current
coordinate system and the mass 'm'.View index ordered by [name] - [type]*/
point op;
op.init(currentcoordsys, p, m);
return op;
}

point point(coordsys R, explicit point M, real m = M.m)
{/*Return the point of 'R' which has the coordinates of 'M' and the mass 'm'.
Do not confuse this routine with the further routine 'changecoordsys'.View index ordered by [name] - [type]*/
point op;
op.init(R, M.coordinates, M.m);
return op;
}

point changecoordsys(coordsys R, point M)
{/*Return the point 'M' in the coordinate system 'coordsys'.
In other words, the returned point marks the same plot as 'M' does.View index ordered by [name] - [type]*/
point op;
coordsys mco = M.coordsys;
op.init(R, R.defaulttorelative(mco.relativetodefault(M.coordinates)), M.m);
return op;
}

pair coordinates(point M)
{/*Return the coordinates of 'M' in its coordinate system.View index ordered by [name] - [type]*/
return M.coordinates;
}

bool samecoordsys(bool warn = true ... point[] M)
{/*Return true iff all the points have the same coordinate system.
If 'warn' is true and the coordinate systems are different, a warning is sent.View index ordered by [name] - [type]*/
bool ret = true;
coordsys t = M[0].coordsys;
for (int i = 1; i < M.length; ++i) {
ret = (t == M[i].coordsys);
if(!ret) break;
t = M[i].coordsys;
}
if(warn && !ret)
warning("coodinatesystem",
"the coordinate system of two objects are not the same.
The operation will be done relative to the default coordinate system.");
return ret;
}

point[] standardizecoordsys(coordsys R = currentcoordsys,
bool warn = true ... point[] M)
{/*Return the points with the same coordinate system 'R'.
If 'warn' is true and the coordinate systems are different, a warning is sent.View index ordered by [name] - [type]*/
point[] op = new point[];
op = M;
if(!samecoordsys(warn ... M))
for (int i = 1; i < M.length; ++i)
op[i] = changecoordsys(R, M[i]);
return op;
}

pair operator cast(point P)
{/*Cast point to pair.View index ordered by [name] - [type]*/
return P.coordsys.relativetodefault(P.coordinates);
}

pair[] operator cast(point[] P)
{/*Cast point[] to pair[].View index ordered by [name] - [type]*/
pair[] op;
for (int i = 0; i < P.length; ++i) {
op.push((pair)P[i]);
}
return op;
}

point operator cast(pair p)
{/*Cast pair to point relatively to the current coordinate
system 'currentcoordsys'.View index ordered by [name] - [type]*/
return point(currentcoordsys, p);
}

point[] operator cast(pair[] p)
{/*Cast pair[] to point[] relatively to the current coordinate
system 'currentcoordsys'.View index ordered by [name] - [type]*/
pair[] op;
for (int i = 0; i < p.length; ++i) {
op.push((point)p[i]);
}
return op;
}

pair locate(point P)
{/*Return the coordinates of 'P' in the default coordinate system.View index ordered by [name] - [type]*/
return P.coordsys * P.coordinates;
}

point locate(pair p)
{/*Return the point in the current coordinate system 'currentcoordsys'.View index ordered by [name] - [type]*/
return p; //automatic casting 'pair to point'.
}

point operator *(real x, explicit point P)
{/*Multiply the coordinates (not the mass) of 'P' by 'x'.View index ordered by [name] - [type]*/
return point(P.coordsys, x * P.coordinates, P.m);
}

point operator /(explicit point P, real x)
{/*Divide the coordinates (not the mass) of 'P' by 'x'.View index ordered by [name] - [type]*/
return point(P.coordsys, P.coordinates/x, P.m);
}

point operator /(real x, explicit point P)
{/*View index ordered by [name] - [type]*/
return point(P.coordsys, x/P.coordinates, P.m);
}

point operator -(explicit point P)
{/*-P. The mass is inchanged.View index ordered by [name] - [type]*/
return point(P.coordsys, -P.coordinates, P.m);
}

point operator +(explicit point P1, explicit point P2)
{/*Provide 'point + point'.
If the two points haven't the same coordinate system, a warning is sent and the
returned point has the default coordinate system 'defaultcoordsys'.
The masses are added.View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(P1, P2);
coordsys R = P[0].coordsys;
return point(R, P[0].coordinates + P[1].coordinates, P1.m + P2.m);
}

point operator +(explicit point P1, explicit pair p2)
{/*Provide 'point + pair'.
The pair 'p2' is supposed to be coordinates relatively to the coordinates system of 'P1'.
The mass is not changed.View index ordered by [name] - [type]*/
coordsys R = currentcoordsys;
return point(R, P1.coordinates + point(R, p2).coordinates, P1.m);
}
point operator +(explicit pair p1, explicit point p2)
{
return p2 + p1;
}

point operator -(explicit point P1, explicit point P2)
{/*Provide 'point - point'.View index ordered by [name] - [type]*/
return P1 + (-P2);
}

point operator -(explicit point P1, explicit pair p2)
{/*Provide 'point - pair'.
The pair 'p2' is supposed to be coordinates relatively to the coordinates system of 'P1'.View index ordered by [name] - [type]*/
return P1 + (-p2);
}
point operator -(explicit pair p1, explicit point P2)
{
return p1 + (-P2);
}

point operator *(transform t, explicit point P)
{/*Provide 'transform * point'.
Note that the transforms scale, xscale, yscale and rotate are carried out relatively
the default coordinate system 'defaultcoordsys' which is not desired for point
defined in an other coordinate system.
On can use scale(real, point), xscale(real, point), yscale(real, point), rotate(real, point),
scaleO(real), xscaleO(real), yscaleO(real) and rotateO(real) (described further)
to change the coordinate system of reference.View index ordered by [name] - [type]*/
coordsys R = P.coordsys;
return point(R, (t * locate(P))/R, P.m);
}

point operator *(explicit point P1, explicit point P2)
{/*Provide 'point * point'.
The resulted mass is the mass of P2View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(P1, P2);
coordsys R = P[0].coordsys;
return point(R, P[0].coordinates * P[1].coordinates, P2.m);
}

point operator *(explicit point P1, explicit pair p2)
{/*Provide 'point * pair'.
The pair 'p2' is supposed to be the coordinates of
the point in the coordinates system of 'P1'.
'pair * point' is also defined.View index ordered by [name] - [type]*/
point P = point(P1.coordsys, p2, P1.m);
return P1 * P;
}
point operator *(explicit pair p1, explicit point p2)
{
return p2 * p1;
}

bool operator ==(explicit point M, explicit point N)
{/*Provide the test 'M == N' wish returns true iff MN < EPSView index ordered by [name] - [type]*/
return abs(locate(M) - locate(N)) < EPS;
}

bool operator !=(explicit point M, explicit point N)
{/*Provide the test 'M != N' wish return true iff MN >= EPSView index ordered by [name] - [type]*/
return !(M == N);
}

guide operator cast(point p)
{/*Cast point to guide.View index ordered by [name] - [type]*/
return locate(p);
}

path operator cast(point p)
{/*Cast point to path.View index ordered by [name] - [type]*/
return locate(p);
}

void dot(picture pic = currentpicture, Label L, explicit point Z,
align align = NoAlign,
string format = defaultformat, pen p = currentpen)
{/*View index ordered by [name] - [type]*/
Label L = L.copy();
L.position(locate(Z));
if(L.s == "") {
if(format == "") format = defaultformat;
L.s = "("+format(format, Z.x)+", "+format(format, Z.y)+")";
}
L.align(align, E);
L.p(p);
dot(pic, locate(Z), p);
}

real abs(coordsys R, pair m)
{/*Return the modulus |m| in the coordinate system 'R'.View index ordered by [name] - [type]*/
return R.abs(m);
}

real abs(explicit point M)
{/*Return the modulus |M| in its coordinate system.View index ordered by [name] - [type]*/
return M.coordsys.abs(M.coordinates);
}

real length(explicit point M)
{/*Return the modulus |M| in its coordinate system (same as 'abs').View index ordered by [name] - [type]*/
return M.coordsys.abs(M.coordinates);
}

point conj(explicit point M)
{/*Conjugate.View index ordered by [name] - [type]*/
return point(M.coordsys, conj(M.coordinates), M.m);
}

real degrees(explicit point M, coordsys R = M.coordsys, bool warn = true)
{/*Return the angle of M (in degrees) relatively to 'R'.View index ordered by [name] - [type]*/
return (degrees(locate(M) - R.O, warn) - degrees(R.i))%360;
}

real angle(explicit point M, coordsys R = M.coordsys, bool warn = true)
{/*Return the angle of M (in radians) relatively to 'R'.View index ordered by [name] - [type]*/
}

bool finite(explicit point p)
{/*Avoid to compute 'finite((pair)(infinite_point))'.View index ordered by [name] - [type]*/
return finite(p.coordinates);
}

real dot(point A, point B)
{/*Return the dot product in the coordinate system of 'A'.View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(A.coordsys, A, B);
return P[0].coordsys.dot(P[0].coordinates, P[1].coordinates);
}

real dot(point A, explicit pair B)
{/*Return the dot product in the default coordinate system.
dot(explicit pair, point) is also defined.View index ordered by [name] - [type]*/
return dot(locate(A), B);
}
real dot(explicit pair A, point B)
{
return dot(A, locate(B));
}

transform rotateO(real a)
{/*Rotation around the origin of the current coordinate system.View index ordered by [name] - [type]*/
return rotate(a, currentcoordsys.O);
};

transform projection(point A, point B)
{/*Return the orthogonal projection on the line (AB).View index ordered by [name] - [type]*/
pair dir = unit(locate(A) - locate(B));
pair a = locate(A);
real cof = dir.x * a.x + dir.y * a.y;
real tx = a.x - dir.x * cof;
real txx = dir.x^2;
real txy = dir.x * dir.y;
real ty = a.y - dir.y * cof;
real tyx = txy;
real tyy = dir.y^2;
transform t = (tx, ty, txx, txy, tyx, tyy);
return t;
}

transform projection(point A, point B, point C, point D, bool safe = false)
{/*Return the (CD) parallel projection on (AB).
If 'safe = true' and (AB)//(CD) return the identity.
If 'safe = false' and (AB)//(CD) return an infinity scaling.View index ordered by [name] - [type]*/
pair a = locate(A);
pair u = unit(locate(B) - locate(A));
pair v = unit(locate(D) - locate(C));
real c = u.x * a.y - u.y * a.x;
real d = (conj(u) * v).y;
if (abs(d) < epsgeo) {
return safe ? identity() : scale(infinity);
}
real tx = c * v.x/d;
real ty = c * v.y/d;
real txx = u.x * v.y/d;
real txy = -u.x * v.x/d;
real tyx = u.y * v.y/d;
real tyy = -u.y * v.x/d;
transform t = (tx, ty, txx, txy, tyx, tyy);
return t;
}

transform scale(real k, point M)
{/*Homothety.View index ordered by [name] - [type]*/
pair P = locate(M);
return shift(P) * scale(k) * shift(-P);
}

transform xscale(real k, point M)
{/*xscale from 'M' relatively to the x - axis of the coordinate system of 'M'.View index ordered by [name] - [type]*/
pair P = locate(M);
real a = degrees(M.coordsys.i);
return (shift(P) * rotate(a)) * xscale(k) * (rotate(-a) * shift(-P));
}

transform yscale(real k, point M)
{/*yscale from 'M' relatively to the y - axis of the coordinate system of 'M'.View index ordered by [name] - [type]*/
pair P = locate(M);
real a = degrees(M.coordsys.j) - 90;
return (shift(P) * rotate(a)) * yscale(k) * (rotate(-a) * shift(-P));
}

transform scale(real k, point A, point B, point C, point D, bool safe = false)
{/*http://fr.wikipedia.org/wiki/Affinit%C3%A9_%28math%C3%A9matiques%29
(help me for English translation...)
If 'safe = true' and (AB)//(CD) return the identity.
If 'safe = false' and (AB)//(CD) return a infinity scaling.View index ordered by [name] - [type]*/
pair a = locate(A);
pair u = unit(locate(B) - locate(A));
pair v = unit(locate(D) - locate(C));
real c = u.x * a.y - u.y * a.x;
real d = (conj(u) * v).y;
real d = (conj(u) * v).y;
if (abs(d) < epsgeo) {
return safe ? identity() : scale(infinity);
}
real tx = (1 - k) * c * v.x/d;
real ty = (1 - k) * c * v.y/d;
real txx = (1 - k) * u.x * v.y/d + k;
real txy = (k - 1) * u.x * v.x/d;
real tyx = (1 - k) * u.y * v.y/d;
real tyy = (k - 1) * u.y * v.x/d + k;
transform t = (tx, ty, txx, txy, tyx, tyy);
return t;
}

transform scaleO(real x)
{/*Homothety from the origin of the current coordinate system.View index ordered by [name] - [type]*/
return scale(x, (0, 0));
}

transform xscaleO(real x)
{/*xscale from the origin and relatively to the current coordinate system.View index ordered by [name] - [type]*/
return scale(x, (0, 0), (0, 1), (0, 0), (1, 0));
}

transform yscaleO(real x)
{/*yscale from the origin and relatively to the current coordinate system.View index ordered by [name] - [type]*/
return scale(x, (0, 0), (1, 0), (0, 0), (0, 1));
}

struct vector
{/*Like a point but casting to pair, adding etc does not take account
of the origin of the coordinate system.*/
point v;/*Coordinates as a point (embed coordinate system and pair).*/
}/*View index ordered by [name] - [type]*/

point operator cast(vector v)
{/*Cast vector 'v' to point 'M' so that OM = v.View index ordered by [name] - [type]*/
return v.v;
}

vector operator cast(pair v)
{/*Cast pair to vector relatively to the current coordinate
system 'currentcoordsys'.View index ordered by [name] - [type]*/
vector ov;
ov.v = point(currentcoordsys, v);
return ov;
}

vector operator cast(explicit point v)
{/*A point can be interpreted like a vector using the code
'(vector)a_point'.View index ordered by [name] - [type]*/
vector ov;
ov.v = v;
return ov;
}

pair operator cast(explicit vector v)
{/*Cast vector to pair (the coordinates of 'v' in the default coordinate system).View index ordered by [name] - [type]*/
return locate(v.v) - v.v.coordsys.O;
}

align operator cast(vector v)
{/*Cast vector to align.View index ordered by [name] - [type]*/
return (pair)v;
}

vector vector(coordsys R = currentcoordsys, pair v)
{/*Return the vector of 'R' which has the coordinates 'v'.View index ordered by [name] - [type]*/
vector ov;
ov.v = point(R, v);
return ov;
}

vector vector(point M)
{/*Return the vector OM, where O is the origin of the coordinate system of 'M'.
Useful to write 'vector(P - M);' instead of '(vector)(P - M)'.View index ordered by [name] - [type]*/
return M;
}

point point(explicit vector u)
{/*Return the point M so that OM = u, where O is the origin of the coordinate system of 'u'.View index ordered by [name] - [type]*/
return u.v;
}

pair locate(explicit vector v)
{/*Return the coordinates of 'v' in the default coordinate system (like casting vector to pair).View index ordered by [name] - [type]*/
return (pair)v;
}

void show(Label L, vector v, pen p = currentpen, arrowbar arrow = Arrow)
{/*Draw the vector v (from the origin of its coordinate system).View index ordered by [name] - [type]*/
coordsys R = v.v.coordsys;
draw(L, R.O--v.v, p, arrow);
}

vector changecoordsys(coordsys R, vector v)
{/*Return the vector 'v' relatively to coordinate system 'R'.View index ordered by [name] - [type]*/
vector ov;
ov.v = point(R, (locate(v) + R.O)/R);
return ov;
}

vector operator *(real x, explicit vector v)
{/*Provide real * vector.View index ordered by [name] - [type]*/
return x * v.v;
}

vector operator /(explicit vector v, real x)
{/*Provide vector/realView index ordered by [name] - [type]*/
return v.v/x;
}

vector operator *(transform t, explicit vector v)
{/*Provide transform * vector.View index ordered by [name] - [type]*/
return t * v.v;
}

vector operator *(explicit point M, explicit vector v)
{/*Provide point * vectorView index ordered by [name] - [type]*/
return M * v.v;
}

point operator +(point M, explicit vector v)
{/*Return 'M' shifted by 'v'.View index ordered by [name] - [type]*/
return shift(locate(v)) * M;
}

point operator -(point M, explicit vector v)
{/*Return 'M' shifted by '-v'.View index ordered by [name] - [type]*/
return shift(-locate(v)) * M;
}

vector operator -(explicit vector v)
{/*Provide -v.View index ordered by [name] - [type]*/
return -v.v;
}

point operator +(explicit pair m, explicit vector v)
{/*The pair 'm' is supposed to be the coordinates of
a point in the current coordinates system 'currentcoordsys'.
Return this point shifted by the vector 'v'.View index ordered by [name] - [type]*/
return locate(m) + v;
}

point operator -(explicit pair m, explicit vector v)
{/*The pair 'm' is supposed to be the coordinates of
a point in the current coordinates system 'currentcoordsys'.
Return this point shifted by the vector '-v'.View index ordered by [name] - [type]*/
return m + (-v);
}

vector operator +(explicit vector v1, explicit vector v2)
{/*Provide vector + vector.
If the two vector haven't the same coordinate system, the returned
vector is relative to the default coordinate system (without warning).View index ordered by [name] - [type]*/
coordsys R = v1.v.coordsys;
if(samecoordsys(false, v1, v2)){R = defaultcoordsys;}
return vector(R, (locate(v1) + locate(v2))/R);
}

vector operator -(explicit vector v1, explicit vector v2)
{/*Provide vector - vector.
If the two vector haven't the same coordinate system, the returned
vector is relative to the default coordinate system (without warning).View index ordered by [name] - [type]*/
return v1 + (-v2);
}

bool operator ==(explicit vector u, explicit vector v)
{/*Return true iff |u - v|<EPS.View index ordered by [name] - [type]*/
return abs(u - v) < EPS;
}

bool collinear(vector u, vector v)
{/*Return 'true' iff the vectors 'u' and 'v' are collinear.View index ordered by [name] - [type]*/
return abs(ypart((conj((pair)u) * (pair)v))) < EPS;
}

vector unit(point M)
{/*Return the unit vector according to the modulus of its coordinate system.View index ordered by [name] - [type]*/
return M/abs(M);
}

vector unit(vector u)
{/*Return the unit vector according to the modulus of its coordinate system.View index ordered by [name] - [type]*/
return u.v/abs(u.v);
}

real degrees(vector v,
coordsys R = v.v.coordsys,
bool warn = true)
{/*Return the angle of 'v' (in degrees) relatively to 'R'.View index ordered by [name] - [type]*/
return (degrees(locate(v), warn) - degrees(R.i))%360;
}

real angle(explicit vector v,
coordsys R = v.v.coordsys,
bool warn = true)
{/*Return the angle of 'v' (in radians) relatively to 'R'.View index ordered by [name] - [type]*/
}

vector conj(explicit vector u)
{/*Conjugate.View index ordered by [name] - [type]*/
return conj(u.v);
}

transform rotate(explicit vector dir)
{/*A rotation in the direction 'dir' limited to [-90, 90]
This is useful for rotating text along a line in the direction dir.
rotate(explicit point dir) is also defined.
View index ordered by [name] - [type]*/
return rotate(locate(dir));
}
transform rotate(explicit point dir){return rotate(locate(vector(dir)));}
// *......................COORDINATES......................*
// *=======================================================*

// *=======================================================*
// *.........................BASES.........................*

point origin = point(defaultcoordsys, (0, 0));/*The origin of the current coordinate system.View index ordered by [name] - [type]*/

point origin(coordsys R = currentcoordsys)
{/*Return the origin of the coordinate system 'R'.View index ordered by [name] - [type]*/
return point(R, (0, 0)); //use automatic casting;
}

real linemargin = 0;/*Margin used to draw lines.View index ordered by [name] - [type]*/

real linemargin()
{/*Return the margin used to draw lines.View index ordered by [name] - [type]*/
return linemargin;
}

pen addpenline = squarecap;/*Add this property to the drawing pen of "finish" lines.View index ordered by [name] - [type]*/
}

pen addpenarc = squarecap;/*Add this property to the drawing pen of arcs.View index ordered by [name] - [type]*/

string defaultmassformat = "$\left(%L;%.4g\right)$";/*Format used to construct the default label of masses.View index ordered by [name] - [type]*/

int sgnd(real x)
{/*Return the -1 if x < 0, 1 if x >= 0.View index ordered by [name] - [type]*/
return (x == 0) ? 1 : sgn(x);
}
int sgnd(int x)
{
return (x == 0) ? 1 : sgn(x);
}

bool defined(point P)
{/*Return true iff the coordinates of 'P' are finite.View index ordered by [name] - [type]*/
return finite(P.coordinates);
}

bool onpath(picture pic = currentpicture, path g, point M, pen p = currentpen)
{/*Return true iff 'M' is on the path drawn with the pen 'p' in 'pic'.View index ordered by [name] - [type]*/
transform t = inverse(pic.calculateTransform());
return intersect(g, shift(locate(M)) * scale(linewidth(p)/2) * t * unitcircle).length > 0;
}

bool sameside(point M, point N, point O)
{/*Return 'true' iff 'M' and 'N' are same side of the point 'O'.View index ordered by [name] - [type]*/
pair m = M, n = N, o = O;
return dot(m - o, n - o) >= -epsgeo;
}

bool between(point M, point O, point N)
{/*Return 'true' iff 'O' is between 'M' and 'N'.View index ordered by [name] - [type]*/
return (!sameside(N, M, O) || M == O || N == O);
}

typedef path pathModifier(path);
pathModifier NoModifier = new path(path g){return g;};

private void Drawline(picture pic = currentpicture, Label L = "", pair P, bool dirP = true, pair Q, bool dirQ = true,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None,
Label legend = "", marker marker = nomarker,
pathModifier pathModifier = NoModifier)
{/* Add the two parameters 'dirP' and 'dirQ' to the native routine
'drawline' of the module 'math'.
Segment [PQ] will be prolonged in direction of P if 'dirP = true', in
direction of Q if 'dirQ = true'.
If 'dirP = dirQ = true', the behavior is that of the native 'drawline'.
Add all the other parameters of 'Draw'.*/
pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
picture opic;
// Reduce the bounds by the size of the pen.
m -= min(p) - (linemargin(), linemargin()); M -= max(p) + (linemargin(), linemargin());

// Calculate the points and direction vector in the transformed space.
t = t * T;
pair z = t * P;
pair q = t * Q;
pair v = q - z;
// path g;
pair ptp, ptq;
real cp = dirP ? 1:0;
real cq = dirQ ? 1:0;
// Handle horizontal and vertical lines.
if(v.x == 0) {
if(m.x <= z.x && z.x <= M.x)
if (dot(v, m - z) < 0) {
ptp = (z.x, z.y + cp * (m.y - z.y));
ptq = (z.x, q.y + cq * (M.y - q.y));
} else {
ptq = (z.x, q.y + cq * (m.y - q.y));
ptp = (z.x, z.y + cp * (M.y - z.y));
}
} else if(v.y == 0) {
if (dot(v, m - z) < 0) {
ptp = (z.x + cp * (m.x - z.x), z.y);
ptq = (q.x + cq * (M.x - q.x), z.y);
} else {
ptq = (q.x + cq * (m.x - q.x), z.y);
ptp = (z.x + cp * (M.x - z.x), z.y);
}
} else {
// Calculate the maximum and minimum t values allowed for the
// parametric equation z + t * v
real mx = (m.x - z.x)/v.x, Mx = (M.x - z.x)/v.x;
real my = (m.y - z.y)/v.y, My = (M.y - z.y)/v.y;
real tmin = max(v.x > 0 ? mx : Mx, v.y > 0 ? my : My);
real tmax = min(v.x > 0 ? Mx : mx, v.y > 0 ? My : my);
pair pmin = z + tmin * v;
pair pmax = z + tmax * v;
if(tmin <= tmax) {
ptp = z + cp * tmin * v;
ptq = z + (cq == 0 ? v:tmax * v);
}
}
path g = ptp--ptq;
if (length(g)>0)
{
if(L.s != "") {
Label lL = L.copy();
if(L.defaultposition) lL.position(Relative(.9));
lL.p(p);
lL.out(opic, g);
}
g = pathModifier(g);
if(linetype(p).length == 0){
pair m = midpoint(g);
pen tp;
tp = dirP ? p : addpenline(p);
draw(opic, pathModifier(m--ptp), tp);
tp = dirQ ? p : addpenline(p);
draw(opic, pathModifier(m--ptq), tp);
} else {
draw(opic, g, p);
}
marker.markroutine(opic, marker.f, g);
arrow(opic, g, p, NoMargin);
}
});
}

void clipdraw(picture pic = currentpicture, Label L = "", path g,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
real xmargin = 0, real ymargin = xmargin,
Label legend = "", marker marker = nomarker)
{/*Draw the path 'g' on 'pic' clipped to the bounding box of 'pic'.View index ordered by [name] - [type]*/
if(L.s != "") {
picture tmp;
label(tmp, L, g, p);
}
pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
// Reduce the bounds by the size of the pen and the margins.
m += min(p) + (xmargin, ymargin); M -= max(p) + (xmargin, ymargin);
path bound = box(m, M);
picture tmp;
draw(tmp, "", t * T * g, align, p, arrow, bar, NoMargin, legend, marker);
clip(tmp, bound);
});
}

void distance(picture pic = currentpicture, Label L = "", point A, point B,
bool rotated = true, real offset = 3mm,
pen p = currentpen, pen joinpen = invisible,
arrowbar arrow = Arrows(NoFill))
{/*Draw arrow between A and B (from FAQ).View index ordered by [name] - [type]*/
pair A = A, B = B;
path g = A--B;
transform Tp = shift(-offset * unit(B - A) * I);
pic.add(new void(frame f, transform t) {
picture opic;
path G = Tp * t * g;
transform id = identity();
transform T = rotated ? rotate(B - A) : id;
Label L = L.copy();
L.align(L.align, Center);
if(abs(ypart((conj(A - B) * L.align.dir))) < epsgeo && L.filltype == NoFill)
L.filltype = UnFill(1);
draw(opic, T * L, G, p, arrow, Bars, PenMargins);
pair Ap = t * A, Bp = t * B;
draw(opic, (Ap--Tp * Ap)^^(Bp--Tp * Bp), joinpen);
}, true);
pic.addBox(min(g), max(g), Tp * min(p), Tp * max(p));
}

real perpfactor = 1;/*Factor for drawing perpendicular symbol.View index ordered by [name] - [type]*/

void perpendicularmark(picture pic = currentpicture, point z,
explicit pair align,
explicit pair dir = E, real size = 0,
pen p = currentpen,
margin margin = NoMargin,
filltype filltype = NoFill)
{/*Draw a perpendicular symbol at z aligned in the direction align
relative to the path z--z + dir.
dir(45 + n * 90), where n in N*, are common values for 'align'.View index ordered by [name] - [type]*/
p = squarecap + p;
if(size == 0) size = perpfactor * 3mm + sqrt(1 + linewidth(p)) - 1;
frame apic;
pair d1 = size * align * unit(dir) * dir(-45);
pair d2 = I * d1;
path g = d1--d1 + d2--d2;
g = margin(g, p).g;
draw(apic, g, p);
if(filltype != NoFill) filltype.fill(apic, (relpoint(g, 0) - relpoint(g, 0.5)+
relpoint(g, 1))--g--cycle, p + solid);
}

void perpendicularmark(picture pic = currentpicture, point z,
vector align,
vector dir = E, real size = 0,
pen p = currentpen,
margin margin = NoMargin,
filltype filltype = NoFill)
{/*Draw a perpendicular symbol at z aligned in the direction align
relative to the path z--z + dir.
dir(45 + n * 90), where n in N, are common values for 'align'.View index ordered by [name] - [type]*/
perpendicularmark(pic, z, (pair)align, (pair)dir, size,
p, margin, filltype);
}

void perpendicularmark(picture pic = currentpicture, point z, explicit pair align, path g,
real size = 0, pen p = currentpen,
margin margin = NoMargin,
filltype filltype = NoFill)
{/*Draw a perpendicular symbol at z aligned in the direction align
relative to the path z--z + dir(g, 0).
dir(45 + n * 90), where n in N, are common values for 'align'.View index ordered by [name] - [type]*/
perpendicularmark(pic, z, align, dir(g, 0), size, p, margin, filltype);
}

void perpendicularmark(picture pic = currentpicture, point z, vector align, path g,
real size = 0, pen p = currentpen,
margin margin = NoMargin,
filltype filltype = NoFill)
{/*Draw a perpendicular symbol at z aligned in the direction align
relative to the path z--z + dir(g, 0).
dir(45 + n * 90), where n in N, are common values for 'align'.View index ordered by [name] - [type]*/
perpendicularmark(pic, z, (pair)align, dir(g, 0), size, p, margin, filltype);
}

void markrightangle(picture pic = currentpicture, point A, point O,
point B, real size = 0, pen p = currentpen,
margin margin = NoMargin,
filltype filltype = NoFill)
{/*Mark the angle AOB with a perpendicular symbol.View index ordered by [name] - [type]*/
pair Ap = A, Bp = B, Op = O;
pair dir = Ap - Op;
real a1 = degrees(dir);
pair align = rotate(-a1) * unit(dir(Op--Ap, Op--Bp));
if (margin == NoMargin)
margin = TrueMargin(linewidth(currentpen)/2, linewidth(currentpen)/2);
perpendicularmark(pic = pic, z = O, align = align,
dir = dir, size = size, p = p,
margin = margin, filltype = filltype);
}

bool simeq(point A, point B, real fuzz = epsgeo)
{/*Return true iff abs(A - B) < fuzz.
This routine is used internally to know if two points are equal, in particular by the operator == in 'point == point'.View index ordered by [name] - [type]*/
return (abs(A - B) < fuzz);
}
bool simeq(point a, real b, real fuzz = epsgeo)
{
coordsys R = a.coordsys;
return (abs(a - point(R, ((pair)b)/R)) < fuzz);
}

pair attract(pair m, path g, real fuzz = 0)
{/*Return the nearest point (A PAIR) of 'm' which is on the path g.
'fuzz' is the argument 'fuzz' of 'intersect'.View index ordered by [name] - [type]*/
if(intersect(m, g, fuzz).length > 0) return m;
pair p;
real step = 1, r = 0;
real[] t;
static real eps = sqrt(realEpsilon);
do {// Find a radius for intersection
r += step;
t = intersect(shift(m) * scale(r) * unitcircle, g);
} while(t.length <= 0);
p = point(g, t[1]);
real rm = 0, rM = r;
while(rM - rm > eps) {
r = (rm + rM)/2;
t = intersect(shift(m) * scale(r) * unitcircle, g, fuzz);
if(t.length <= 0) {
rm = r;
} else {
rM = r;
p = point(g, t[1]);
}
}
return p;
}

point attract(point M, path g, real fuzz = 0)
{/*Return the nearest point (A POINT) of 'M' which is on the path g.
'fuzz' is the argument 'fuzz' of 'intersect'.View index ordered by [name] - [type]*/
return point(M.coordsys, attract(locate(M), g)/M.coordsys);
}

real[] intersect(path g, explicit pair p, real fuzz = 0)
{/*View index ordered by [name] - [type]*/
fuzz = fuzz <= 0 ? sqrt(realEpsilon) : fuzz;
real[] or;
real r = realEpsilon;
do{
or = intersect(g, shift(p) * scale(r) * unitcircle, fuzz);
r *= 2;
} while(or.length == 0);
return or;
}

real[] intersect(path g, explicit point P, real fuzz = epsgeo)
{/*View index ordered by [name] - [type]*/
return intersect(g, locate(P), fuzz);
}
// *.........................BASES.........................*
// *=======================================================*

// *=======================================================*
// *.........................LINES.........................*

struct line
{/*This structure provides the objects line, semi - line and segment oriented from A to B.
All the calculus with this structure will be as exact as Asymptote can do.
For a full precision, you must not cast 'line' to 'path' excepted for drawing routines.*/

restricted point A,B;/*Two line's points with same coordinate system.*/
bool extendA,extendB;/*If true,extend 'l' in direction of A (resp. B).*/
restricted vector u,v;/*u = unit(AB) = direction vector,v = normal vector.*/
restricted real a,b,c;/*Coefficients of the equation ax + by + c = 0 in the coordinate system of 'A'.*/
restricted real slope, origin;/*Slope and ordinate at the origin.*/

line copy()
{/*Copy a line in a new instance.*/
line l = new line;
l.A = A;
l.B = B;
l.a = a;
l.b = b;
l.c = c;
l.slope = slope;
l.origin = origin;
l.u = u;
l.v = v;
l.extendA = extendA;
l.extendB = extendB;
return l;
}

void init(point A, bool extendA = true, point B, bool extendB = true)
{/*Initialize line.
If 'extendA' is true, the "line" is infinite in the direction of A.*/
point[] P = standardizecoordsys(A, B);
this.A = P[0];
this.B = P[1];
this.a = B.y - A.y;
this.b = A.x - B.x;
this.c = A.y * B.x - A.x * B.y;
this.slope= (this.b == 0) ? infinity : -this.a/this.b;
this.origin = (this.b == 0) ? (this.c == 0) ? 0:infinity : -this.c/this.b;
this.u = unit(P[1]-P[0]);
//     int tmp = sgnd(this.slope);
//     this.u = (dot((pair)this.u, N) >= 0) ? tmp * this.u : -tmp * this.u;
this.v = rotate(90, point(P[0].coordsys, (0, 0))) * this.u;
this.extendA = extendA;
this.extendB = extendB;
}
}/*View index ordered by [name] - [type]*/

line line(point A, bool extendA = true, point B, bool extendB = true)
{/*Return the line passing through 'A' and 'B'.
If 'extendA' is true, the "line" is infinite in the direction of A.
A "line" can be half-line or segment.View index ordered by [name] - [type]*/
if (A == B) abort("line: the points must be distinct.");
line l;
l.init(A, extendA, B, extendB);
return l;
}

struct segment
{/*.*/
restricted point A, B;// Extremity.
restricted vector u, v;// u = direction vector, v = normal vector.
restricted real a, b, c;// Coefficients of the équation ax + by + c = 0
restricted real slope, origin;
segment copy()
{
segment s = new segment;
s.A = A;
s.B = B;
s.a = a;
s.b = b;
s.c = c;
s.slope = slope;
s.origin = origin;
s.u = u;
s.v = v;
return s;
}

void init(point A, point B)
{
line l;
l.init(A, B);
this.A = l.A; this.B = l.B;
this.a = l.a; this.b = l.b; this.c = l.c;
this.slope = l.slope; this.origin = l.origin;
this.u = l.u; this.v = l.v;
}
}/*View index ordered by [name] - [type]*/

segment segment(point A, point B)
{/*Return the segment whose the extremities are A and B.View index ordered by [name] - [type]*/
segment s;
s.init(A, B);
return s;
}

real length(segment s)
{/*Return the length of 's'.View index ordered by [name] - [type]*/
return abs(s.A - s.B);
}

line operator cast(segment s)
{/*A segment is casted to a "finite line".View index ordered by [name] - [type]*/
return line(s.A, false, s.B, false);
}

segment operator cast(line l)
{/*Cast line 'l' to segment [l.A l.B].View index ordered by [name] - [type]*/
return segment(l.A, l.B);
}

line operator *(transform t, line l)
{/*Provide transform * lineView index ordered by [name] - [type]*/
return line(t * l.A, l.extendA, t * l.B, l.extendB);
}

line operator /(line l, real x)
{/*Provide l/x.
Return the line passing through l.A/x and l.B/x.View index ordered by [name] - [type]*/
return line(l.A/x, l.extendA, l.B/x, l.extendB);
}
line operator /(line l, int x){return line(l.A/x, l.B/x);}

line operator *(real x, line l)
{/*Provide x * l.
Return the line passing through x * l.A and x * l.B.View index ordered by [name] - [type]*/
return line(x * l.A, l.extendA, x * l.B, l.extendB);
}
line operator *(int x, line l){return line(x * l.A, l.extendA, x * l.B, l.extendB);}

line operator *(point M, line l)
{/*Provide point * line.
Return the line passing through unit(M) * l.A and unit(M) * l.B.View index ordered by [name] - [type]*/
return line(unit(M) * l.A, l.extendA, unit(M) * l.B, l.extendB);
}

line operator +(line l, vector u)
{/*Provide line + vector (and so line + point).
Return the line 'l' shifted by 'u'.View index ordered by [name] - [type]*/
return line(l.A + u, l.extendA, l.B + u, l.extendB);
}

line operator -(line l, vector u)
{/*Provide line - vector (and so line - point).
Return the line 'l' shifted by '-u'.View index ordered by [name] - [type]*/
return line(l.A - u, l.extendA, l.B - u, l.extendB);
}

line[] operator ^^(line l1, line l2)
{/*Provide line^^line.
Return the line array {l1, l2}.View index ordered by [name] - [type]*/
line[] ol;
ol.push(l1); ol.push(l2);
return ol;
}

line[] operator ^^(line l1, line[] l2)
{/*Provide line^^line[].
Return the line array {l1, l2[0], l2[1]...}.
line[]^^line is also defined.View index ordered by [name] - [type]*/
line[] ol;
ol.push(l1);
for (int i = 0; i < l2.length; ++i) {
ol.push(l2[i]);
}
return ol;
}
line[] operator ^^(line[] l2, line l1)
{
line[] ol = l2;
ol.push(l1);
return ol;
}

line[] operator ^^(line l1[], line[] l2)
{/*Provide line[]^^line[].
Return the line array {l1[0], l1[1], ..., l2[0], l2[1], ...}.View index ordered by [name] - [type]*/
line[] ol = l1;
for (int i = 0; i < l2.length; ++i) {
ol.push(l2[i]);
}
return ol;
}

bool sameside(point M, point P, line l)
{/*Return 'true' iff 'M' and 'N' are same side of the line (or on the line) 'l'.View index ordered by [name] - [type]*/
pair A = l.A, B = l.B, m = M, p = P;
pair mil = (A + B)/2;
pair mA = rotate(90, mil) * A;
pair mB = rotate(-90, mil) * A;
return (abs(m - mA) <= abs(m - mB)) == (abs(p - mA) <= abs(p - mB));
// transform proj = projection(l.A, l.B);
// point Mp = proj * M;
// point Pp = proj * P;
// dot(Mp);dot(Pp);
// return dot(locate(Mp - M), locate(Pp - P)) >= 0;
}

line line(segment s)
{/*Return the line passing through 's.A'
and 's.B'.View index ordered by [name] - [type]*/
return line(s.A, s.B);
}

segment segment(line l)
{/*Return the segment whose extremities
are 'l.A' and 'l.B'.View index ordered by [name] - [type]*/
return segment(l.A, l.B);
}

point midpoint(segment s)
{/*Return the midpoint of 's'.View index ordered by [name] - [type]*/
return 0.5 * (s.A + s.B);
}

void write(explicit line l)
{/*Write some informations about 'l'.View index ordered by [name] - [type]*/
write("A = "+(string)((pair)l.A));
write("Extend A = "+(l.extendA ? "true" : "false"));
write("B = "+(string)((pair)l.B));
write("Extend B = "+(l.extendB ? "true" : "false"));
write("u = "+(string)((pair)l.u));
write("v = "+(string)((pair)l.v));
write("a = "+(string) l.a);
write("b = "+(string) l.b);
write("c = "+(string) l.c);
write("slope = "+(string) l.slope);
write("origin = "+(string) l.origin);
}

void write(explicit segment s)
{/*Write some informations about 's'.View index ordered by [name] - [type]*/
write("A = "+(string)((pair)s.A));
write("B = "+(string)((pair)s.B));
write("u = "+(string)((pair)s.u));
write("v = "+(string)((pair)s.v));
write("a = "+(string) s.a);
write("b = "+(string) s.b);
write("c = "+(string) s.c);
write("slope = "+(string) s.slope);
write("origin = "+(string) s.origin);
}

bool operator ==(line l1, line l2)
{/*Provide the test 'line == line'.View index ordered by [name] - [type]*/
return (collinear(l1.u, l2.u) &&
abs(ypart((locate(l1.A) - locate(l1.B))/(locate(l1.A) - locate(l2.B)))) < epsgeo &&
l1.extendA == l2.extendA && l1.extendB == l2.extendB);
}

bool operator !=(line l1, line l2)
{/*Provide the test 'line != line'.View index ordered by [name] - [type]*/
return !(l1 == l2);
}

bool operator @(point m, line l)
{/*Provide the test 'point @ line'.
Return true iff 'm' is on the 'l'.View index ordered by [name] - [type]*/
point M = changecoordsys(l.A.coordsys, m);
if (abs(l.a * M.x + l.b * M.y + l.c) >= epsgeo) return false;
if (l.extendA && l.extendB) return true;
if (!l.extendA && !l.extendB) return between(l.A, M, l.B);
if (l.extendA) return sameside(M, l.A, l.B);
return sameside(M, l.B, l.A);
}

coordsys coordsys(line l)
{/*Return the coordinate system in which 'l' is defined.View index ordered by [name] - [type]*/
return l.A.coordsys;
}

line reverse(line l)
{/*Permute the points 'A' and 'B' of 'l' and so its orientation.View index ordered by [name] - [type]*/
return line(l.B, l.extendB, l.A, l.extendA);
}

line extend(line l)
{/*Return the infinite line passing through 'l.A' and 'l.B'.View index ordered by [name] - [type]*/
line ol = l.copy();
ol.extendA = true;
ol.extendB = true;
return ol;
}

line complementary(explicit line l)
{/*Return the complementary of a half-line with respect of
the full line 'l'.View index ordered by [name] - [type]*/
if (l.extendA && l.extendB)
abort("complementary: the parameter is not a half-line.");
point origin = l.extendA ? l.B : l.A;
point ptdir = l.extendA ?
rotate(180, l.B) * l.A : rotate(180, l.A) * l.B;
return line(origin, false, ptdir);
}

line[] complementary(explicit segment s)
{/*Return the two half-lines of origin 's.A' and 's.B' respectively.View index ordered by [name] - [type]*/
line[] ol = new line[2];
ol[0] = complementary(line(s.A, false, s.B));
ol[1] = complementary(line(s.A, s.B, false));
return ol;
}

line Ox(coordsys R = currentcoordsys)
{/*Return the x-axis of 'R'.View index ordered by [name] - [type]*/
return line(point(R, (0, 0)), point(R, E));
}

restricted line Ox = Ox();/*the x-axis of
the default coordinate system.View index ordered by [name] - [type]*/

line Oy(coordsys R = currentcoordsys)
{/*Return the y-axis of 'R'.View index ordered by [name] - [type]*/
return line(point(R, (0, 0)), point(R, N));
}

restricted line Oy = Oy();/*the y-axis of
the default coordinate system.View index ordered by [name] - [type]*/

line line(real a, point A = point(currentcoordsys, (0, 0)))
{/*Return the line passing through 'A' with an
angle (in the coordinate system of A) 'a' in degrees.
line(point, real) is also defined.View index ordered by [name] - [type]*/
return line(A, A + point(A.coordsys, A.coordsys.polar(1, radians(a))));
}
line line(point A = point(currentcoordsys, (0, 0)), real a)
{
return line(a, A);
}
line line(int a, point A = point(currentcoordsys, (0, 0)))
{
return line((real)a, A);
}

line line(coordsys R = currentcoordsys, real slope, real origin)
{/*Return the line defined by slope and y-intercept relative to 'R'.View index ordered by [name] - [type]*/
if (slope == infinity || slope == -infinity)
abort("The slope is infinite. Please, use the routine 'vline'.");
return line(point(R, (0, origin)), point(R, (1, origin + slope)));
}

line line(coordsys R = currentcoordsys, real a, real b, real c)
{/*Retrun the line defined by equation relative to 'R'.View index ordered by [name] - [type]*/
if (a == 0 && b == 0) abort("line: inconsistent equation...");
pair M;
M = (a == 0) ? (0, -c/b) : (-c/a, 0);
return line(point(R, M), point(R, M + (-b, a)));
}

line vline(coordsys R = currentcoordsys)
{/*Return a vertical line in 'R' passing through the origin of 'R'.View index ordered by [name] - [type]*/
point P = point(R, (0, 0));
point PP = point(R, (R.O + N)/R);
return line(P, PP);
}

restricted line vline = vline();/*The vertical line in the current coordinate system passing
through the origin of this system.View index ordered by [name] - [type]*/

line hline(coordsys R = currentcoordsys)
{/*Return a horizontal line in 'R' passing through the origin of 'R'.View index ordered by [name] - [type]*/
point P = point(R, (0, 0));
point PP = point(R, (R.O + E)/R);
return line(P, PP);
}

line hline = hline();/*The horizontal line in the current coordinate system passing
through the origin of this system.View index ordered by [name] - [type]*/

line changecoordsys(coordsys R, line l)
{/*Return the line 'l' in the coordinate system 'R'.View index ordered by [name] - [type]*/
point A = changecoordsys(R, l.A);
point B = changecoordsys(R, l.B);
return line(A, B);
}

transform scale(real k, line l1, line l2, bool safe = false)
{/*Return the dilatation with respect to
'l1' in the direction of 'l2'.View index ordered by [name] - [type]*/
return scale(k, l1.A, l1.B, l2.A, l2.B, safe);
}

transform reflect(line l)
{/*Return the reflect about the line 'l'.View index ordered by [name] - [type]*/
return reflect((pair)l.A, (pair)l.B);
}

transform reflect(line l1, line l2, bool safe = false)
{/*Return the reflect about the line
'l1' in the direction of 'l2'.View index ordered by [name] - [type]*/
return scale(-1.0, l1, l2, safe);
}

point[] intersectionpoints(line l, path g)
{/*Return all points of intersection of the line 'l' with the path 'g'.View index ordered by [name] - [type]*/
// TODO utiliser la version 1.44 de intersections(path g, pair p, pair q)
// real [] t = intersections(g, l.A, l.B);
// coordsys R = coordsys(l);
// return sequence(new point(int n){return point(R, point(g, t[n])/R);}, t.length);
real [] t;
pair[] op;
pair A = l.A;
pair B = l.B;
real dy = B.y - A.y,
dx = A.x - B.x,
lg = length(g);

for (int i = 0; i < lg; ++i)
{
pair z0 = point(g, i),
z1 = point(g, i + 1),
c0 = postcontrol(g, i),
c1 = precontrol(g, i + 1),
t3 = z1 - z0 - 3 * c1 + 3 * c0,
t2 = 3 * z0 + 3 * c1 - 6 * c0,
t1 = 3 * c0 - 3z0;
real a = dy * t3.x + dx * t3.y,
b = dy * t2.x + dx * t2.y,
c = dy * t1.x + dx * t1.y,
d = dy * z0.x + dx * z0.y + A.y * B.x - A.x * B.y;

t = cubicroots(a, b, c, d);
for (int j = 0; j < t.length; ++j)
if (
t[j]>=0
&& (
t[j]<1
|| (
t[j] == 1
&& (i == lg - 1)
&& !cyclic(g)
)
)
) {
op.push(point(g, i + t[j]));
}
}

point[] opp;
for (int i = 0; i < op.length; ++i)
opp.push(point(coordsys(l), op[i]/coordsys(l)));
return opp;
}

point intersectionpoint(line l1, line l2)
{/*Return the point of intersection of line 'l1' with 'l2'.
If 'l1' and 'l2' have an infinity or none point of intersection,
this routine return (infinity, infinity).View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(l1.A, l1.B, l2.A, l2.B);
coordsys R = P[0].coordsys;
pair p = extension(P[0], P[1], P[2], P[3]);
if(finite(p)){
point p = point(R, p/R);
if (p @ l1 && p @ l2) return p;
}
return point(R, (infinity, infinity));
}

line parallel(point M, line l)
{/*Return the line parallel to 'l' passing through 'M'.View index ordered by [name] - [type]*/
point A, B;
if (M.coordsys != coordsys(l))
{
A = changecoordsys(M.coordsys, l.A);
B = changecoordsys(M.coordsys, l.B);
} else {A = l.A;B = l.B;}
return line(M, M - A + B);
}

line parallel(point M, explicit vector dir)
{/*Return the line of direction 'dir' and passing through 'M'.View index ordered by [name] - [type]*/
return line(M, M + locate(dir));
}

line parallel(point M, explicit pair dir)
{/*Return the line of direction 'dir' and passing through 'M'.View index ordered by [name] - [type]*/
return line(M, M + vector(currentcoordsys, dir));
}

bool parallel(line l1, line l2, bool strictly = false)
{/*Return 'true' if 'l1' and 'l2' are (strictly ?) parallel.View index ordered by [name] - [type]*/
bool coll = collinear(l1.u, l2.u);
return strictly ? coll && (l1 != l2) : coll;
}

bool concurrent(... line[] l)
{/*Returns true if all the lines 'l' are concurrent.View index ordered by [name] - [type]*/
if (l.length < 3) abort("'concurrent' needs at least for three lines ...");
pair point = intersectionpoint(l[0], l[1]);
bool conc;
for (int i = 2; i < l.length; ++i) {
pair pt = intersectionpoint(l[i - 1], l[i]);
conc = simeq(pt, point);
if (!conc) break;
}
return conc;
}

transform projection(line l)
{/*Return the orthogonal projection on 'l'.View index ordered by [name] - [type]*/
return projection(l.A, l.B);
}

transform projection(line l1, line l2, bool safe = false)
{/*Return the projection on (AB) in parallel of (CD).
If 'safe = true' and (l1)//(l2) return the identity.
If 'safe = false' and (l1)//(l2) return a infinity scaling.View index ordered by [name] - [type]*/
return projection(l1.A, l1.B, l2.A, l2.B, safe);
}

transform vprojection(line l, bool safe = false)
{/*Return the projection on 'l' in parallel of N--S.
If 'safe' is 'true' the projected point keeps the same place if 'l'
is vertical.View index ordered by [name] - [type]*/
coordsys R = defaultcoordsys;
return projection(l, line(point(R, N), point(R, S)), safe);
}

transform hprojection(line l, bool safe = false)
{/*Return the projection on 'l' in parallel of E--W.
If 'safe' is 'true' the projected point keeps the same place if 'l'
is horizontal.View index ordered by [name] - [type]*/
coordsys R = defaultcoordsys;
return projection(l, line(point(R, E), point(R, W)), safe);
}

line perpendicular(point M, line l)
{/*Return the perpendicular line of 'l' passing through 'M'.View index ordered by [name] - [type]*/
point Mp = projection(l) * M;
point A = Mp == l.A ? l.B : l.A;
return line(Mp, rotate(90, Mp) * A);
}

line perpendicular(point M, explicit vector normal)
{/*Return the line passing through 'M'
whose normal is \param{normal}.View index ordered by [name] - [type]*/
return perpendicular(M, line(M, M + locate(normal)));
}

line perpendicular(point M, explicit pair normal)
{/*Return the line passing through 'M'
whose normal is \param{normal} (given in the currentcoordsys).View index ordered by [name] - [type]*/
return perpendicular(M, line(M, M + vector(currentcoordsys, normal)));
}

bool perpendicular(line l1, line l2)
{/*Return 'true' if 'l1' and 'l2' are perpendicular.View index ordered by [name] - [type]*/
return abs(dot(locate(l1.u), locate(l2.u))) < epsgeo ;
}

real angle(line l, coordsys R = coordsys(l))
{/*Return the angle of the oriented line 'l',
in radian, in the interval ]-pi, pi] and relatively to 'R'.View index ordered by [name] - [type]*/
return angle(l.u, R, false);
}

real degrees(line l, coordsys R = coordsys(l))
{/*Returns the angle of the oriented line 'l' in degrees,
in the interval [0, 360[ and relatively to 'R'.View index ordered by [name] - [type]*/
return degrees(angle(l, R));
}

real sharpangle(line l1, line l2)
{/*Return the measure in radians of the sharp angle formed by 'l1' and 'l2'.View index ordered by [name] - [type]*/
vector u1 = l1.u;
vector u2 = (dot(l1.u, l2.u) < 0) ? -l2.u : l2.u;
real a12 = angle(locate(u2)) - angle(locate(u1));
a12 = a12%(sgnd(a12) * pi);
if (a12 <= -pi/2) {
a12 += pi;
} else if (a12 > pi/2) {
a12 -= pi;
}
return a12;
}

real angle(line l1, line l2)
{/*Return the measure in radians of oriented angle (l1.u, l2.u).View index ordered by [name] - [type]*/
return angle(locate(l2.u)) - angle(locate(l1.u));
}

real degrees(line l1, line l2)
{/*Return the measure in degrees of the
angle formed by the oriented lines 'l1' and 'l2'.View index ordered by [name] - [type]*/
return degrees(angle(l1, l2));
}

real sharpdegrees(line l1, line l2)
{/*Return the measure in degrees of the sharp angle formed by 'l1' and 'l2'.View index ordered by [name] - [type]*/
return degrees(sharpangle(l1, l2));
}

line bisector(line l1, line l2, real angle = 0, bool sharp = true)
{/*Return the bisector of the angle formed by 'l1' and 'l2'
rotated by the angle 'angle' (in degrees) around intersection point of 'l1' with 'l2'.
If 'sharp' is true (the default), this routine returns the bisector of the sharp angle.
Note that the returned line inherit of coordinate system of 'l1'.View index ordered by [name] - [type]*/
line ol;
if (l1 == l2) return l1;
point A = intersectionpoint(l1, l2);
if (finite(A)) {
if(sharp) ol = rotate(sharpdegrees(l1, l2)/2 + angle, A) * l1;
else {
coordsys R = coordsys(l1);
pair a = A, b = A + l1.u, c = A + l2.u;
pair pp = extension(a, a + dir(a--b, a--c), b, b + dir(b--a, b--c));
return rotate(angle, A) * line(A, point(R, pp/R));
}
} else {
ol = l1;
}
return ol;
}

line sector(int n = 2, int p = 1, line l1, line l2, real angle = 0, bool sharp = true)
{/*Return the p-th nth-sector of the angle
formed by the oriented line 'l1' and 'l2'
rotated by the angle 'angle' (in degrees) around the intersection point of 'l1' with 'l2'.
If 'sharp' is true (the default), this routine returns the bisector of the sharp angle.
Note that the returned line inherit of coordinate system of 'l1'.View index ordered by [name] - [type]*/
line ol;
if (l1 == l2) return l1;
point A = intersectionpoint(l1, l2);
if (finite(A)) {
if(sharp) ol = rotate(p * sharpdegrees(l1, l2)/n + angle, A) * l1;
else {
ol = rotate(p * degrees(l1, l2)/n + angle, A) * l1;
}
} else {
ol = l1;
}
return ol;
}

line bisector(point A, point B, point C, point D, real angle = 0, bool sharp = true)
{/*Return the bisector of the angle formed by the lines (AB) and (CD).
Look at bisector(line, line, real, bool).View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(A, B, C, D);
return bisector(line(P[0], P[1]), line(P[2], P[3]), angle, sharp);
}

line bisector(segment s, real angle = 0)
{/*Return the bisector of the segment line 's' rotated by 'angle' (in degrees) around the
midpoint of 's'.View index ordered by [name] - [type]*/
coordsys R = coordsys(s);
point m = midpoint(s);
vector dir = rotateO(90) * unit(s.A - m);
return rotate(angle, m) * line(m + dir, m - dir);
}

line bisector(point A, point B, real angle = 0)
{/*Return the bisector of the segment line [AB] rotated by 'angle' (in degrees) around the
midpoint of [AB].View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(A, B);
return bisector(segment(P[0], P[1]), angle);
}

real distance(point M, line l)
{/*Return the distance from 'M' to 'l'.
distance(line, point) is also defined.View index ordered by [name] - [type]*/
point A = changecoordsys(defaultcoordsys, l.A);
point B = changecoordsys(defaultcoordsys, l.B);
line ll = line(A, B);
pair m = locate(M);
return abs(ll.a * m.x + ll.b * m.y + ll.c)/sqrt(ll.a^2 + ll.b^2);
}

real distance(line l, point M)
{
return distance(M, l);
}

void draw(picture pic = currentpicture, Label L = "",
line l, bool dirA = l.extendA, bool dirB = l.extendB,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None,
Label legend = "", marker marker = nomarker,
pathModifier pathModifier = NoModifier)
{/*Draw the line 'l' without altering the size of picture pic.
The boolean parameters control the infinite section.
The global variable 'linemargin' (default value is 0) allows to modify
the bounding box in which the line must be drawn.View index ordered by [name] - [type]*/
if(!(dirA || dirB)) draw(l.A--l.B, invisible);// l is a segment.
Drawline(pic, L, l.A, dirP = dirA, l.B, dirQ = dirB,
align, p, arrow,
legend, marker, pathModifier);
}

void draw(picture pic = currentpicture, Label[] L = new Label[], line[] l,
align align = NoAlign, pen[] p = new pen[],
arrowbar arrow = None,
Label[] legend = new Label[], marker marker = nomarker,
pathModifier pathModifier = NoModifier)
{/*Draw each lines with the corresponding pen.View index ordered by [name] - [type]*/
for (int i = 0; i < l.length; ++i) {
draw(pic, L.length>0 ? L[i] : "", l[i],
align, p = p.length>0 ? p[i] : currentpen,
arrow, legend.length>0 ? legend[i] : "", marker,
pathModifier);
}
}

void draw(picture pic = currentpicture, Label[] L = new Label[], line[] l,
align align = NoAlign, pen p,
arrowbar arrow = None,
Label[] legend = new Label[], marker marker = nomarker,
pathModifier pathModifier = NoModifier)
{/*Draw each lines with the same pen 'p'.View index ordered by [name] - [type]*/
pen[] tp = sequence(new pen(int i){return p;}, l.length);
draw(pic, L, l, align, tp, arrow, legend, marker, pathModifier);
}

void show(picture pic = currentpicture, line l, pen p = red)
{/*Draw some informations of 'l'.View index ordered by [name] - [type]*/
dot("$A$", (pair)l.A, align = -locate(l.v), p);
dot("$B$", (pair)l.B, align = -locate(l.v), p);
draw(l, dotted);
draw("$\vec{u}$", locate(l.A)--locate(l.A + l.u), p, Arrow);
draw("$\vec{v}$", locate(l.A)--locate(l.A + l.v), p, Arrow);
}

point[] sameside(point M, line l1, line l2)
{/*Return two points on 'l1' and 'l2' respectively.
The first point is from the same side of M relatively to 'l2',
the second point is from the same side of M relatively to 'l1'.View index ordered by [name] - [type]*/
point[] op;
coordsys R1 = coordsys(l1);
coordsys R2 = coordsys(l2);
if (parallel(l1, l2)) {
op.push(projection(l1) * M);
op.push(projection(l2) * M);
} else {
point O = intersectionpoint(l1, l2);
if (M @ l2) op.push((sameside(M, O + l1.u, l2)) ? O + l1.u : rotate(180, O) * (O + l1.u));
else op.push(projection(l1, l2) * M);
if (M @ l1) op.push((sameside(M, O + l2.u, l1)) ? O + l2.u : rotate(180, O) * (O + l2.u));
else {op.push(projection(l2, l1) * M);}
}
return op;
}

void markangle(picture pic = currentpicture,
Label L = "", int n = 1, real radius = 0, real space = 0,
explicit line l1, explicit line l2, explicit pair align = dir(1),
arrowbar arrow = None, pen p = currentpen,
filltype filltype = NoFill,
margin margin = NoMargin, marker marker = nomarker)
{/*Mark the angle (l1, l2) aligned in the direction 'align' relative to 'l1'.
Commune values for 'align' are dir(real).View index ordered by [name] - [type]*/
if (parallel(l1, l2, true)) return;
real al = degrees(l1, defaultcoordsys);
pair O, A, B;
real d = degrees(locate(l1.u));
align = rotate(d) * align;
if (l1 == l2) {
O = midpoint(segment(l1.A, l1.B));
A = l1.A;B = l1.B;
if (sameside(rotate(sgn(angle(B-A)) * 45, O) * A, O + align, l1)) {radius = -radius;}
} else {
O = intersectionpoint(extend(l1), extend(l2));
pair R = O + align;
point [] ss = sameside(point(coordsys(l1), R/coordsys(l1)), l1, l2);
A = ss[0];
B = ss[1];
}
markangle(pic = pic, L = L, n = n, radius = radius, space = space,
O = O, A = A, B = B,
arrow = arrow, p = p, filltype = filltype,
margin = margin, marker = marker);
}

void markangle(picture pic = currentpicture,
Label L = "", int n = 1, real radius = 0, real space = 0,
explicit line l1, explicit line l2, explicit vector align,
arrowbar arrow = None, pen p = currentpen,
filltype filltype = NoFill,
margin margin = NoMargin, marker marker = nomarker)
{/*Mark the angle (l1, l2) in the direction 'dir' given relatively to 'l1'.View index ordered by [name] - [type]*/
markangle(pic, L, n, radius, space, l1, l2, (pair)align, arrow,
p, filltype, margin, marker);
}

// void markangle(picture pic = currentpicture,
//                Label L = "", int n = 1, real radius = 0, real space = 0,
//                explicit line l1, explicit line l2,
//                arrowbar arrow = None, pen p = currentpen,
//                filltype filltype = NoFill,
//                margin margin = NoMargin, marker marker = nomarker)
// {/*Mark the oriented angle (l1, l2).View index ordered by [name] - [type]*/
//   if (parallel(l1, l2, true)) return;
//   real al = degrees(l1, defaultcoordsys);
//   pair O, A, B;
//   real d = degrees(locate(l1.u));
//   if (l1 == l2) {
//     O = midpoint(segment(l1.A, l1.B));
//   } else {
//     O = intersectionpoint(extend(l1), extend(l2));
//   }
//   A = O + locate(l1.u);
//   B = O + locate(l2.u);
//   markangle(pic = pic, L = L, n = n, radius = radius, space = space,
//             O = O, A = A, B = B,
//             arrow = arrow, p = p, filltype = filltype,
//             margin = margin, marker = marker);
// }

void perpendicularmark(picture pic = currentpicture, line l1, line l2,
real size = 0, pen p = currentpen, int quarter = 1,
margin margin = NoMargin, filltype filltype = NoFill)
{/*Draw a right angle at the intersection point of lines and
aligned in the 'quarter' nth quarter of circle formed by 'l1.u' and
'l2.u'.View index ordered by [name] - [type]*/
point P = intersectionpoint(l1, l2);
pair align = rotate(90 * (quarter - 1)) * dir(45);
perpendicularmark(P, align, locate(l1.u), size, p, margin, filltype);
}
// *.........................LINES.........................*
// *=======================================================*

// *=======================================================*
// *........................CONICS.........................*

struct bqe

real[] a;/*a[0] * x^2 + a[1] * x * y + a[2] * y^2 + a[3] * x + a[4] * y + a[5] = 0*/
coordsys coordsys;
}/*View index ordered by [name] - [type]*/

bqe bqe(coordsys R = currentcoordsys,
real a, real b, real c, real d, real e, real f)
a[0] * x^2 + a[1] * x * y + a[2] * y^2 + a[3] * x + a[4] * y + a[5] = 0
relatively to the coordinate system R.View index ordered by [name] - [type]*/
bqe obqe;
obqe.coordsys = R;
obqe.a = new real[] {a, b, c, d, e, f};
return obqe;
}

bqe changecoordsys(coordsys R, bqe bqe)
{/*Returns the bivariate quadratic equation relatively to 'R'.View index ordered by [name] - [type]*/
pair i = coordinates(changecoordsys(R, vector(defaultcoordsys,
bqe.coordsys.i)));
pair j = coordinates(changecoordsys(R, vector(defaultcoordsys,
bqe.coordsys.j)));
pair O = coordinates(changecoordsys(R, point(defaultcoordsys,
bqe.coordsys.O)));
real a = bqe.a[0], b = bqe.a[1], c = bqe.a[2], d = bqe.a[3], f = bqe.a[4], g = bqe.a[5];
real ux = i.x, uy = i.y;
real vx = j.x, vy = j.y;
real ox = O.x, oy = O.y;
real D = ux * vy - uy * vx;
real ap = (a * vy^2 - b * uy * vy + c * uy^2)/D^2;
real bpp = (-2 * a * vx * vy + b * ux * vy + b * uy * vx - 2 * c * ux * uy)/D^2;
real cp = (a * vx^2 - b * ux * vx + c * ux^2)/D^2;
real dp = (-2a * ox * vy^2 + 2a * oy * vx * vy + 2b * ox * uy * vy-
b * oy * ux * vy - b * oy * uy * vx - 2c * ox * uy^2 + 2c * oy * uy * ux)/D^2+
(d * vy - f * uy)/D;
real fp = (2a * ox * vx * vy - b * ox * ux * vy - 2a * oy * vx^2-
b * ox * uy * vx + 2 * b * oy * ux * vx + 2c * ox * ux * uy - 2c * oy * ux^2)/D^2+
(f * ux - d * vx)/D;
g = (a * ox^2 * vy^2 - 2a * ox * oy * vx * vy - b * ox^2 * uy * vy + b * ox * oy * ux * vy+
a * oy^2 * vx^2 + b * ox * oy * uy * vx - b * oy^2 * ux * vx + c * ox^2 * uy^2-
2 * c * ox * oy * ux * uy + c * oy^2 * ux^2)/D^2+
(d * oy * vx + f * ox * uy - d * ox * vy - f * oy * ux)/D + g;
bqe obqe;
obqe.a = approximate(new real[] {ap, bpp, cp, dp, fp, g});
obqe.coordsys = R;
return obqe;
}

bqe bqe(point M1, point M2, point M3, point M4, point M5)
{/*Return the bqe of conic passing through the five points (if possible).View index ordered by [name] - [type]*/
coordsys R;
pair[] pts;
if (samecoordsys(M1, M2, M3, M4, M5)) {
R = M1.coordsys;
pts= new pair[] {M1.coordinates, M2.coordinates, M3.coordinates, M4.coordinates, M5.coordinates};
} else {
R = defaultcoordsys;
pts= new pair[] {M1, M2, M3, M4, M5};
}
real[][] M;
real[] x;
bqe bqe;
bqe.coordsys = R;
for (int i = 0; i < 5; ++i) {// Try a = -1
M[i] = new real[] {pts[i].x * pts[i].y, pts[i].y^2, pts[i].x, pts[i].y, 1};
x[i] = pts[i].x^2;
}
if(abs(determinant(M)) < 1e-5) {// Try c = -1
for (int i = 0; i < 5; ++i) {
M[i] = new real[] {pts[i].x^2, pts[i].x * pts[i].y, pts[i].x, pts[i].y, 1};
x[i] = pts[i].y^2;
}
real[] coef = solve(M, x);
bqe.a = new real[] {coef[0], coef[1], -1, coef[2], coef[3], coef[4]};
} else {
real[] coef = solve(M, x);
bqe.a = new real[] {-1, coef[0], coef[1], coef[2], coef[3], coef[4]};
}
bqe.a = approximate(bqe.a);
return bqe;
}

bool samecoordsys(bool warn = true ... bqe[] bqes)
{/*Return true if all the bivariate quadratic equations have the same coordinate system.View index ordered by [name] - [type]*/
bool ret = true;
coordsys t = bqes[0].coordsys;
for (int i = 1; i < bqes.length; ++i) {
ret = (t == bqes[i].coordsys);
if(!ret) break;
t = bqes[i].coordsys;
}
if(warn && !ret)
warning("coodinatesystem",
"the coordinate system of two bivariate quadratic equations are not
the same. The operation will be done relatively to the default coordinate
system.");
return ret;
}

real[] realquarticroots(real a, real b, real c, real d, real e)
{/*Return the real roots of the quartic equation ax^4 + b^x3 + cx^2 + dx = 0.View index ordered by [name] - [type]*/
static real Fuzz = sqrt(realEpsilon);
pair[] zroots = quarticroots(a, b, c, d, e);
real[] roots;
real p(real x){return a * x^4 + b * x^3 + c * x^2 + d * x + e;}
real prime(real x){return 4 * a * x^3 + 3 * b * x^2 + 2 * c * x + d;}
real x;
bool search = true;
int n;
{
bool exist = false;
for (int i = 0; i < roots.length; ++i) {
if(abs(roots[i]-x) < 1e-5) {exist = true; break;}
}
if(!exist) roots.push(x);
}
for(int i = 0; i < zroots.length; ++i) {
if(zroots[i].y == 0 || abs(p(zroots[i].x)) < Fuzz) addroot(zroots[i].x);
else {
if(abs(zroots[i].y) < 1e-3) {
x = zroots[i].x;
search = true;
n = 200;
while(search) {
real tx = abs(p(x)) < Fuzz ? x : newton(iterations = n, p, prime, x);
if(tx < realMax) {
if(abs(p(tx)) < Fuzz) {
search = false;
} else if(n < 200) n *=2;
else {
search = false;
}
} else search = false; //It's not a real root.
}
}
}
}
return roots;
}

point[] intersectionpoints(bqe bqe1, bqe bqe2)
{/*Return the interscetion of the two conic sections whose equations are 'bqe1' and 'bqe2'.View index ordered by [name] - [type]*/
coordsys R = bqe1.coordsys;
bqe lbqe1, lbqe2;
real[] a, b;
if(R != bqe2.coordsys) {
R = currentcoordsys;
a = changecoordsys(R, bqe1).a;
b = changecoordsys(R, bqe2).a;
} else {
a = bqe1.a;
b = bqe2.a;
}
static real e = 100 * sqrt(realEpsilon);
real[] x, y, c;
point[] P;
if(abs(a[0]-b[0]) > e || abs(a[1]-b[1]) > e || abs(a[2]-b[2]) > e) {
c = new real[] {-2 * a[0]*a[2]*b[0]*b[2]+a[0]*a[2]*b[1]^2 - a[0]*a[1]*b[2]*b[1]+a[1]^2 * b[0]*b[2]-
a[2]*a[1]*b[0]*b[1]+a[0]^2 * b[2]^2 + a[2]^2 * b[0]^2,
-a[2]*a[1]*b[0]*b[4]-a[2]*a[4]*b[0]*b[1]-a[1]*a[3]*b[2]*b[1]+2 * a[0]*a[2]*b[1]*b[4]-
a[0]*a[1]*b[2]*b[4]+a[1]^2 * b[2]*b[3]-2 * a[2]*a[3]*b[0]*b[2]-2 * a[0]*a[2]*b[2]*b[3]+
a[2]*a[3]*b[1]^2 - a[2]*a[1]*b[1]*b[3]+2 * a[1]*a[4]*b[0]*b[2]+2 * a[2]^2 * b[0]*b[3]-
a[0]*a[4]*b[2]*b[1]+2 * a[0]*a[3]*b[2]^2,
-a[3]*a[4]*b[2]*b[1]+a[2]*a[5]*b[1]^2 - a[1]*a[5]*b[2]*b[1]-a[1]*a[3]*b[2]*b[4]+
a[1]^2 * b[2]*b[5]-2 * a[2]*a[3]*b[2]*b[3]+2 * a[2]^2 * b[0]*b[5]+2 * a[0]*a[5]*b[2]^2 + a[3]^2 * b[2]^2-
2 * a[2]*a[5]*b[0]*b[2]+2 * a[1]*a[4]*b[2]*b[3]-a[2]*a[4]*b[1]*b[3]-2 * a[0]*a[2]*b[2]*b[5]+
a[2]^2 * b[3]^2 + 2 * a[2]*a[3]*b[1]*b[4]-a[2]*a[4]*b[0]*b[4]+a[4]^2 * b[0]*b[2]-a[2]*a[1]*b[3]*b[4]-
a[2]*a[1]*b[1]*b[5]-a[0]*a[4]*b[2]*b[4]+a[0]*a[2]*b[4]^2,
-a[4]*a[5]*b[2]*b[1]+a[2]*a[3]*b[4]^2 + 2 * a[3]*a[5]*b[2]^2 - a[2]*a[1]*b[4]*b[5]-
a[2]*a[4]*b[3]*b[4]+2 * a[2]^2 * b[3]*b[5]-2 * a[2]*a[3]*b[2]*b[5]-a[3]*a[4]*b[2]*b[4]-
2 * a[2]*a[5]*b[2]*b[3]-a[2]*a[4]*b[1]*b[5]+2 * a[1]*a[4]*b[2]*b[5]-a[1]*a[5]*b[2]*b[4]+
a[4]^2 * b[2]*b[3]+2 * a[2]*a[5]*b[1]*b[4],
-2 * a[2]*a[5]*b[2]*b[5]+a[4]^2 * b[2]*b[5]+a[5]^2 * b[2]^2 - a[4]*a[5]*b[2]*b[4]+a[2]*a[5]*b[4]^2+
a[2]^2 * b[5]^2 - a[2]*a[4]*b[4]*b[5]};
x = realquarticroots(c[0], c[1], c[2], c[3], c[4]);
} else {
if(abs(b[4]-a[4]) > e){
real D = (b[4]-a[4])^2;
c = new real[] {(a[0]*b[4]^2 + (-a[1]*b[3]-2 * a[0]*a[4]+a[1]*a[3]) * b[4]+a[2]*b[3]^2+
(a[1]*a[4]-2 * a[2]*a[3]) * b[3]+a[0]*a[4]^2 - a[1]*a[3]*a[4]+a[2]*a[3]^2)/D,
-((a[1]*b[4]-2 * a[2]*b[3]-a[1]*a[4]+2 * a[2]*a[3]) * b[5]-a[3]*b[4]^2 + (a[4]*b[3]-a[1]*a[5]+a[3]*a[4]) * b[4]+(2 * a[2]*a[5]-a[4]^2) * b[3]+(a[1]*a[4]-2 * a[2]*a[3]) * a[5])/D,
a[2]*(a[5]-b[5])^2/D + a[4]*(a[5]-b[5])/(b[4]-a[4]) + a[5]};
} else {
if(abs(a[3]-b[3]) > e) {
real D = b[3]-a[3];
c = new real[] {a[2], (-a[1]*b[5] + a[4]*b[3] + a[1]*a[5] - a[3]*a[4])/D,
a[0]*(a[5]-b[5])^2/D^2 + a[3]*(a[5]-b[5])/D + a[5]};
for (int i = 0; i < y.length; ++i) {
c = new real[] {a[0], a[1]*y[i]+a[3], a[2]*y[i]^2 + a[4]*y[i]+a[5]};
for (int j = 0; j < x.length; ++j) {
if(abs(b[0]*x[j]^2 + b[1]*x[j]*y[i]+b[2]*y[i]^2 + b[3]*x[j]+b[4]*y[i]+b[5]) < 1e-5)
P.push(point(R, (x[j], y[i])));
}
}
return P;
} else {
if(abs(a[5]-b[5]) < e) abort("intersectionpoints: intersection of identical conics.");
}
}
}
for (int i = 0; i < x.length; ++i) {
c = new real[] {a[2], a[1]*x[i]+a[4], a[0]*x[i]^2 + a[3]*x[i]+a[5]};
for (int j = 0; j < y.length; ++j) {
if(abs(b[0]*x[i]^2 + b[1]*x[i]*y[j]+b[2]*y[j]^2 + b[3]*x[i]+b[4]*y[j]+b[5]) < 1e-5)
P.push(point(R, (x[i], y[j])));
}
}
return P;
}

struct conic
{/**/
real e, p, h;/*BE CAREFUL: h = distance(F, D) and p = h * e (http://en.wikipedia.org/wiki/Ellipse)
While http://mathworld.wolfram.com/ takes p = distance(F,D).*/
point F;/*Focus.*/
line D;/*Directrix.*/
line[] l;/*Case of degenerated conic (not yet implemented !).*/
}/*View index ordered by [name] - [type]*/

bool degenerate(conic c)
{
return !finite(c.p) || !finite(c.h);
}

/*ANCconic conic(point, line, real)ANC*/
conic conic(point F, line l, real e)
{/*DOC
The conic section define by the eccentricity 'e', the focus 'F'
and the directrix 'l'.
Note that an eccentricity equal to 0 defines a circle centered at F,
with a radius equal at the distance from 'F' to 'l'.
If the coordinate system of 'F' and 'l' are not identical, the conic is
attached to 'defaultcoordsys'.
DOC*/
if(e < 0) abort("conic: 'e' can't be negative.");
conic oc;
point[] P = standardizecoordsys(F, l.A, l.B);
line ll;
ll = line(P[1], P[2]);
oc.e = e < epsgeo ? 0 : e; // Handle case of circle.
oc.F = P[0];
oc.D = ll;
oc.h = distance(P[0], ll);
oc.p = abs(e) < epsgeo ? oc.h : e * oc.h;
return oc;
}

struct circle
{/*All the calculus with this structure will be as exact as Asymptote can do.
For a full precision, you must not cast 'circle' to 'path' excepted for drawing routines.*/

point C;/*Center*/
line l;/*If the radius is infinite, this line is used instead of circle.*/
}/*View index ordered by [name] - [type]*/

bool degenerate(circle c)
{
return !finite(c.r);
}

line line(circle c){
if(finite(c.r)) abort("Circle can not be casted to line here.");
return c.l;
}

struct ellipse
{/*Look at http://mathworld.wolfram.com/Ellipse.html*/

restricted point F1,F2,C;/*Foci and center.*/
restricted real a,b,c,e,p;
restricted real angle;/*Value is degrees(F1 - F2).*/
restricted line D1,D2;/*Directrices.*/
line l;/*If one axis is infinite, this line is used instead of ellipse.*/

void init(point f1, point f2, real a)
{/*Ellipse given by foci and semimajor axis*/
point[] P = standardizecoordsys(f1, f2);
this.F1 = P[0];
this.F2 = P[1];
this.angle = abs(P[1]-P[0]) < 10 * epsgeo ? 0 : degrees(P[1]-P[0]);
this.C = (P[0] + P[1])/2;
this.a = a;
if(!finite(a)) {
this.l = line(P[0], P[1]);
this.b = infinity;
this.e = 0;
this.c = 0;
} else {
this.c = abs(C - P[0]);
this.b = this.c < epsgeo ? a : sqrt(a^2 - c^2); // Handle case of circle.
this.e = this.c < epsgeo ? 0 : this.c/a; // Handle case of circle.
if(this.e >= 1) abort("ellipse.init: wrong parameter: e >= 1.");
this.p = a * (1 - this.e^2);
if (this.c != 0) {// directrix is not set for a circle.
point A = this.C + (a^2/this.c) * unit(P[0]-this.C);
this.D1 = line(A, A + rotateO(90) * unit(A - this.C));
this.D2 = reverse(rotate(180, C) * D1);
}
}
}
}/*View index ordered by [name] - [type]*/

bool degenerate(ellipse el)
{
return (!finite(el.a) || !finite(el.b));
}

struct parabola
{/*Look at http://mathworld.wolfram.com/Parabola.html*/
restricted point F,V;/*Focus and vertex*/
restricted real a,p,e = 1;
restricted real angle;/*Angle, in degrees, of the line (FV).*/
restricted line D;/*Directrix*/
pair bmin, bmax;/*The (left, bottom) and (right, top) coordinates of region bounding box for drawing the parabola.
If unset the current picture bounding box is used instead.*/

void init(point F, line directrix)
{/*Parabola given by focus and directrix.*/
point[] P = standardizecoordsys(F, directrix.A, directrix.B);
line l = line(P[1], P[2]);
this.F = P[0];
this.D = l;
this.a = distance(P[0], l)/2;
this.p = 2 * a;
this.V = 0.5 * (F + projection(D) * P[0]);
this.angle = degrees(F - V);
}
}/*View index ordered by [name] - [type]*/

struct hyperbola
{/*Look at http://mathworld.wolfram.com/Hyperbola.html*/
restricted point F1,F2;/*Foci.*/
restricted point C,V1,V2;/*Center and vertices.*/
restricted real a,b,c,e,p;/**/
restricted real angle;/*Angle,in degrees,of the line (F1F2).*/
restricted line D1,D2,A1,A2;/*Directrices and asymptotes.*/
pair bmin, bmax; /*The (left, bottom) and (right, top) coordinates of region bounding box for drawing the hyperbola.
If unset the current picture bounding box is used instead.*/

void init(point f1, point f2, real a)
{/*Hyperbola given by foci and semimajor axis.*/
point[] P = standardizecoordsys(f1, f2);
this.F1 = P[0];
this.F2 = P[1];
this.angle = degrees(F2 - F1);
this.a = a;
this.C = (P[0] + P[1])/2;
this.c = abs(C - P[0]);
this.e = this.c/a;
if(this.e <= 1) abort("hyperbola.init: wrong parameter: e <= 1.");
this.b = a * sqrt(this.e^2 - 1);
this.p = a * (this.e^2 - 1);
point A = this.C + (a^2/this.c) * unit(P[0]-this.C);
this.D1 = line(A, A + rotateO(90) * unit(A - this.C));
this.D2 = reverse(rotate(180, C) * D1);
this.V1 = C + a * unit(F1 - C);
this.V2 = C + a * unit(F2 - C);
this.A1 = line(C, V1 + b * unit(rotateO(-90) * (C - V1)));
this.A2 = line(C, V1 + b * unit(rotateO(90) * (C - V1)));
}
}/*View index ordered by [name] - [type]*/

int conicnodesfactor = 1;/*Factor for the node number of all conics.View index ordered by [name] - [type]*/

int circlenodesnumberfactor = 100;/*Factor for the node number of circles.View index ordered by [name] - [type]*/

int circlenodesnumber(real r)
{/*Return the number of nodes for drawing a circle of radius 'r'.View index ordered by [name] - [type]*/
if (circlenodesnumberfactor < 100)
warning("circlenodesnumberfactor",
"variable 'circlenodesnumberfactor' may be too small.");
int oi = ceil(circlenodesnumberfactor * abs(r)^0.1);
oi = 45 * floor(oi/45);
return oi == 0 ? 4 : conicnodesfactor * oi;
}

int circlenodesnumber(real r, real angle1, real angle2)
{/*Return the number of nodes to draw a circle arc.View index ordered by [name] - [type]*/
return (r > 0) ?
ceil(circlenodesnumber(r) * abs(angle1 - angle2)/360) :
ceil(circlenodesnumber(r) * abs((1 - abs(angle1 - angle2)/360)));
}

int ellipsenodesnumberfactor = 250;/*Factor for the node number of ellispe (non-circle).View index ordered by [name] - [type]*/

int ellipsenodesnumber(real a, real b)
{/*Return the number of nodes to draw a ellipse of axis 'a' and 'b'.View index ordered by [name] - [type]*/
if (ellipsenodesnumberfactor < 250)
write("ellipsenodesnumberfactor",
"variable 'ellipsenodesnumberfactor' maybe too small.");
int tmp = circlenodesnumberfactor;
circlenodesnumberfactor = ellipsenodesnumberfactor;
int oi = circlenodesnumber(max(abs(a), abs(b))/min(abs(a), abs(b)));
circlenodesnumberfactor = tmp;
return conicnodesfactor * oi;
}

int ellipsenodesnumber(real a, real b, real angle1, real angle2, bool dir)
{/*Return the number of nodes to draw an ellipse arc.View index ordered by [name] - [type]*/
real d;
real da = angle2 - angle1;
if(dir) {
d = angle1 < angle2 ? da : 360 + da;
} else {
d = angle1 < angle2 ? -360 + da : da;
}
int n = floor(ellipsenodesnumber(a, b) * abs(d)/360);
return n < 5 ? 5 : n;
}

int parabolanodesnumberfactor = 100;/*Factor for the number of nodes of parabolas.View index ordered by [name] - [type]*/

int parabolanodesnumber(parabola p, real angle1, real angle2)
{/*Return the number of nodes for drawing a parabola.View index ordered by [name] - [type]*/
return conicnodesfactor * floor(0.01 * parabolanodesnumberfactor * abs(angle1 - angle2));
}

int hyperbolanodesnumberfactor = 100;/*Factor for the number of nodes of hyperbolas.View index ordered by [name] - [type]*/

int hyperbolanodesnumber(hyperbola h, real angle1, real angle2)
{/*Return the number of nodes for drawing an hyperbola.View index ordered by [name] - [type]*/
return conicnodesfactor * floor(0.01 * hyperbolanodesnumberfactor * abs(angle1 - angle2)/h.e);
}

conic operator +(conic c, explicit point M)
{/*View index ordered by [name] - [type]*/
return conic(c.F + M, c.D + M, c.e);
}

conic operator -(conic c, explicit point M)
{/*View index ordered by [name] - [type]*/
return conic(c.F - M, c.D - M, c.e);
}

conic operator +(conic c, explicit pair m)
{/*View index ordered by [name] - [type]*/
point M = point(c.F.coordsys, m);
return conic(c.F + M, c.D + M, c.e);
}

conic operator -(conic c, explicit pair m)
{/*View index ordered by [name] - [type]*/
point M = point(c.F.coordsys, m);
return conic(c.F - M, c.D - M, c.e);
}

conic operator +(conic c, vector v)
{/*View index ordered by [name] - [type]*/
return conic(c.F + v, c.D + v, c.e);
}

conic operator -(conic c, vector v)
{/*View index ordered by [name] - [type]*/
return conic(c.F - v, c.D - v, c.e);
}

coordsys coordsys(conic co)
{/*Return the coordinate system of 'co'.View index ordered by [name] - [type]*/
return co.F.coordsys;
}

conic changecoordsys(coordsys R, conic co)
{/*Change the coordinate system of 'co' to 'R'View index ordered by [name] - [type]*/
line l = changecoordsys(R, co.D);
point F = changecoordsys(R, co.F);
return conic(F, l, co.e);
}

typedef path polarconicroutine(conic co, real angle1, real angle2, int n, bool direction);/*Routine type used to draw conics from 'angle1' to 'angle2'*/

path arcfromfocus(conic co, real angle1, real angle2, int n = 400, bool direction = CCW)
{/*Return the path of the conic section 'co' from angle1 to angle2 in degrees,
drawing in the given direction, with n nodes.View index ordered by [name] - [type]*/
guide op;
if (n < 1) return op;
if (angle1 > angle2) {
path g = arcfromfocus(co, angle2, angle1, n, !direction);
return g == nullpath ? g : reverse(g);
}
point O = projection(co.D) * co.F;
pair i = unit(locate(co.F) - locate(O));
pair j = rotate(90) * i;
coordsys Rp = cartesiansystem(co.F, i, j);
real a2 = direction ? radians(angle2) : radians(angle1) + 2 * pi;
real step = n == 1 ? 0 : (a2 - a1)/(n - 1);
real a, r;
for (int i = 0; i < n; ++i) {
a = a1 + i * step;
if(co.e >= 1) {
r = 1 - co.e * cos(a);
if(r > epsgeo) {
r = co.p/r;
op = op--Rp * Rp.polar(r, a);
}
} else {
r = co.p/(1 - co.e * cos(a));
op = op..Rp * Rp.polar(r, a);
}
}
if(co.e < 1 && abs(abs(a2 - a1) - 2 * pi) < epsgeo) op = (path)op..cycle;

return (direction ? op : op == nullpath ? op :reverse(op));
}

polarconicroutine currentpolarconicroutine = arcfromfocus;/*Default routine used to cast conic section to path.View index ordered by [name] - [type]*/

point angpoint(conic co, real angle)
{/*Return the point of 'co' whose the angular (in degrees)
coordinate is 'angle' (mesured from the focus of 'co', relatively
to its 'natural coordinate system').View index ordered by [name] - [type]*/
coordsys R = coordsys(co);
return point(R, point(arcfromfocus(co, angle, angle, 1, CCW), 0)/R);
}

bool operator @(point M, conic co)
{/*Return true iff 'M' on 'co'.View index ordered by [name] - [type]*/
if(co.e == 0) return abs(abs(co.F - M) - co.p) < 10 * epsgeo;
return abs(co.e * distance(M, co.D) - abs(co.F - M)) < 10 * epsgeo;
}

coordsys coordsys(ellipse el)
{/*Return the coordinate system of 'el'.View index ordered by [name] - [type]*/
return el.F1.coordsys;
}

coordsys canonicalcartesiansystem(ellipse el)
{/*Return the canonical cartesian system of the ellipse 'el'.View index ordered by [name] - [type]*/
if(degenerate(el)) return cartesiansystem(el.l.A, el.l.u, el.l.v);
pair O = locate(el.C);
pair i = el.e == 0 ? el.C.coordsys.i : unit(locate(el.F1) - O);
pair j = rotate(90) * i;
return cartesiansystem(O, i, j);
}

coordsys canonicalcartesiansystem(parabola p)
{/*Return the canonical cartesian system of a parabola,
so that Origin = vertex of 'p' and directrix: x = -a.View index ordered by [name] - [type]*/
point A = projection(p.D) * p.F;
pair O = locate((A + p.F)/2);
pair i = unit(locate(p.F) - O);
pair j = rotate(90) * i;
return cartesiansystem(O, i, j);
}

coordsys canonicalcartesiansystem(hyperbola h)
{/*Return the canonical cartesian system of an hyperbola.View index ordered by [name] - [type]*/
pair O = locate(h.C);
pair i = unit(locate(h.F2) - O);
pair j = rotate(90) * i;
return cartesiansystem(O, i, j);
}

ellipse ellipse(point F1, point F2, real a)
{/*Return the ellipse whose the foci are 'F1' and 'F2'
and the semimajor axis is 'a'.View index ordered by [name] - [type]*/
ellipse oe;
oe.init(F1, F2, a);
return oe;
}

restricted bool byfoci = true, byvertices = false;/*Constants useful for the routine 'hyperbola(point P1, point P2, real ae, bool byfoci = byfoci)'View index ordered by [name] - [type]*/

hyperbola hyperbola(point P1, point P2, real ae, bool byfoci = byfoci)
{/*if 'byfoci = true':
return the hyperbola whose the foci are 'P1' and 'P2'
and the semimajor axis is 'ae'.
else return the hyperbola whose vertexes are 'P1' and 'P2' with eccentricity 'ae'.View index ordered by [name] - [type]*/
hyperbola oh;
point[] P = standardizecoordsys(P1, P2);
if(byfoci) {
oh.init(P[0], P[1], ae);
} else {
real a = abs(P[0]-P[1])/2;
vector V = unit(P[0]-P[1]);
point F1 = P[0] + a * (ae - 1) * V;
point F2 = P[1]-a * (ae - 1) * V;
oh.init(F1, F2, a);
}
return oh;
}

ellipse ellipse(point F1, point F2, point M)
{/*Return the ellipse passing through 'M' whose the foci are 'F1' and 'F2'.View index ordered by [name] - [type]*/
point P[] = standardizecoordsys(false, F1, F2, M);
real a = abs(F1 - M) + abs(F2 - M);
return ellipse(F1, F2, finite(a) ? a/2 : a);
}

ellipse ellipse(point C, real a, real b, real angle = 0)
{/*Return the ellipse centered at 'C' with semimajor axis 'a' along C--C + dir(angle),
semiminor axis 'b' along the perpendicular.View index ordered by [name] - [type]*/
ellipse oe;
coordsys R = C.coordsys;
angle += degrees(R.i);
if(a < b) {angle += 90; real tmp = a; a = b; b = tmp;}
if(finite(a) && finite(b)) {
real c = sqrt(abs(a^2 - b^2));
point f1, f2;
if(abs(a - b) < epsgeo) {
f1 = C; f2 = C;
} else {
f1 = point(R, (locate(C) + rotate(angle) * (-c, 0))/R);
f2 = point(R, (locate(C) + rotate(angle) * (c, 0))/R);
}
oe.init(f1, f2, a);
} else {
if(finite(b) || !finite(a)) oe.init(C, C + R.polar(1, angle), infinity);
else oe.init(C, C + R.polar(1, 90 + angle), infinity);
}
return oe;
}

ellipse ellipse(bqe bqe)
{/*Return the ellipse a[0] * x^2 + a[1] * xy + a[2] * y^2 + a[3] * x + a[4] * y + a[5] = 0
given in the coordinate system of 'bqe' with a[i] = bque.a[i].
http://mathworld.wolfram.com/Ellipse.html.View index ordered by [name] - [type]*/
bqe lbqe = changecoordsys(defaultcoordsys, bqe);
real a = lbqe.a[0], b = lbqe.a[1]/2, c = lbqe.a[2], d = lbqe.a[3]/2, f = lbqe.a[4]/2, g = lbqe.a[5];
coordsys R = bqe.coordsys;
string message = "ellipse: the given equation is not an equation of an ellipse.";
real u = b^2 * g + d^2 * c + f^2 * a;
real delta = a * c * g + b * f * d + d * b * f - u;
if(abs(delta) < epsgeo) abort(message);
real j = b^2 - a * c;
real i = a + c;
real dd = j * (sgnd(c - a) * sqrt((a - c)^2 + 4 * (b^2)) - c-a);
real ddd = j * (-sgnd(c - a) * sqrt((a - c)^2 + 4 * (b^2)) - c-a);

if(abs(ddd) < epsgeo || abs(dd) < epsgeo ||
j >= -epsgeo || delta/sgnd(i) > 0) abort(message);

real x = (c * d - b * f)/j, y = (a * f - b * d)/j;
// real dir = abs(b) < epsgeo ? 0 : pi/2-0.5 * acot(0.5 * (c-a)/b);
real dir = abs(b) < epsgeo ? 0 : 0.5 * acot(0.5 * (c - a)/b);
if(dir * (c - a) * b < 0) dir = dir < 0 ? dir + pi/2 : dir - pi/2;
real cd = cos(dir), sd = sin(dir);
real t = a * cd^2 - 2 * b * cd * sd + c * sd^2;
real tt = a * sd^2 + 2 * b * cd * sd + c * cd^2;
real gg = -g + ((d * cd - f * sd)^2)/t + ((d * sd + f * cd)^2)/tt;
t = t/gg; tt = tt/gg;
// The equation of the ellipse is t * (x - center.x)^2 + tt * (y - center.y)^2 = 1;
real aa, bb;
aa = sqrt(2 * (u - 2 * b * d * f - a * c * g)/dd);
bb = sqrt(2 * (u - 2 * b * d * f - a * c * g)/ddd);
a = t > tt ? max(aa, bb) : min(aa, bb);
b = t > tt ? min(aa, bb) : max(aa, bb);
return ellipse(point(R, (x, y)/R),
a, b, degrees(pi/2 - dir - angle(R.i)));
}

ellipse ellipse(point M1, point M2, point M3, point M4, point M5)
{/*Return the ellipse passing through the five points (if possible)View index ordered by [name] - [type]*/
return ellipse(bqe(M1, M2, M3, M4, M5));
}

bool inside(ellipse el, point M)
{/*Return 'true' iff 'M' is inside 'el'.View index ordered by [name] - [type]*/
return abs(el.F1 - M) + abs(el.F2 - M) - 2 * el.a < -epsgeo;
}

bool inside(parabola p, point M)
{/*Return 'true' if 'M' is inside 'p'.View index ordered by [name] - [type]*/
return distance(p.D, M) - abs(p.F - M) > epsgeo;
}

parabola parabola(point F, line l)
{/*Return the parabola whose focus is 'F' and directrix is 'l'.View index ordered by [name] - [type]*/
parabola op;
op.init(F, l);
return op;
}

parabola parabola(point F, point vertex)
{/*Return the parabola whose focus is 'F' and vertex is 'vertex'.View index ordered by [name] - [type]*/
parabola op;
point[] P = standardizecoordsys(F, vertex);
point A = rotate(180, P[1]) * P[0];
point B = A + rotateO(90) * unit(P[1]-A);
op.init(P[0], line(A, B));
return op;
}

parabola parabola(point F, real a, real angle)
{/*Return the parabola whose focus is F, latus rectum is 4a and
the angle of the axis of symmetry (in the coordinate system of F) is 'angle'.View index ordered by [name] - [type]*/
parabola op;
coordsys R = F.coordsys;
point A = F - point(R, R.polar(2a, radians(angle)));
point B = A + point(R, R.polar(1, radians(90 + angle)));
op.init(F, line(A, B));
return op;
}

bool isparabola(bqe bqe)
{/*Return true iff 'bqe' is the equation of a parabola.View index ordered by [name] - [type]*/
bqe lbqe = changecoordsys(defaultcoordsys, bqe);
real a = lbqe.a[0], b = lbqe.a[1]/2, c = lbqe.a[2], d = lbqe.a[3]/2, f = lbqe.a[4]/2, g = lbqe.a[5];
real delta = a * c * g + b * f * d + d * b * f - (b^2 * g + d^2 * c + f^2 * a);
return (abs(delta) > epsgeo && abs(b^2 - a * c) < epsgeo);
}

parabola parabola(bqe bqe)
{/*Return the parabola a[0]x^2 + a[1]xy + a[2]y^2 + a[3]x + a[4]y + a[5]] = 0 (a[n] means bqe.a[n]).
http://mathworld.wolfram.com/Parabola.htmlView index ordered by [name] - [type]*/
bqe lbqe = changecoordsys(defaultcoordsys, bqe);
real a = lbqe.a[0], b = lbqe.a[1]/2, c = lbqe.a[2], d = lbqe.a[3]/2, f = lbqe.a[4]/2, g = lbqe.a[5];
string message = "parabola: the given equation is not an equation of a parabola.";
real delta = a * c * g + b * f * d + d * b * f - (b^2 * g + d^2 * c + f^2 * a);
if(abs(delta) < 10 * epsgeo || abs(b^2 - a * c) > 10 * epsgeo) abort(message);
real dir = abs(b) < epsgeo ? 0 : 0.5 * acot(0.5 * (c - a)/b);
if(dir * (c - a) * b < 0) dir = dir < 0 ? dir + pi/2 : dir - pi/2;
real cd = cos(dir), sd = sin(dir);
real ap = a * cd^2 - 2 * b * cd * sd + c * sd^2;
real cp = a * sd^2 + 2 * b * cd * sd + c * cd^2;
real dp = d * cd - f * sd;
real fp = d * sd + f * cd;
real gp = g;
parabola op;
coordsys R = bqe.coordsys;
// The equation of the parabola is ap * x'^2 + cp * y'^2 + 2dp * x'+2fp * y'+gp = 0
if (abs(ap) < epsgeo) {/* directrix parallel to the rotated(dir) y-axis
equation: (y-vertex.y)^2 = 4 * a * (x-vertex)
*/
pair pvertex = rotate(degrees(-dir)) * (0.5(-gp + fp^2/cp)/dp, -fp/cp);
real a = -0.5 * dp/cp;
point vertex = point(R, pvertex/R);
point focus = point(R, (pvertex + a * expi(-dir))/R);
op = parabola(focus, vertex);

} else {/* directrix parallel to the rotated(dir) x-axis
equation: (x-vertex)^2 = 4 * a * (y-vertex.y)
*/
pair pvertex = rotate(degrees(-dir)) * (-dp/ap, 0.5 * (-gp + dp^2/ap)/fp);
real a = -0.5 * fp/ap;
point vertex = point(R, pvertex/R);
point focus = point(R, (pvertex + a * expi(pi/2 - dir))/R);
op = parabola(focus, vertex);
}
return op;
}

parabola parabola(point M1, point M2, point M3, line l)
{/*Return the parabola passing through the three points with its directix
parallel to the line 'l'.View index ordered by [name] - [type]*/
coordsys R;
pair[] pts;
if (samecoordsys(M1, M2, M3)) {
R = M1.coordsys;
} else {
R = defaultcoordsys;
}
real gle = degrees(l);
coordsys Rp = cartesiansystem(R.O, rotate(gle) * R.i, rotate(gle) * R.j);
pts = new pair[] {coordinates(changecoordsys(Rp, M1)),
coordinates(changecoordsys(Rp, M2)),
coordinates(changecoordsys(Rp, M3))};
real[][] M;
real[] x;
for (int i = 0; i < 3; ++i) {
M[i] = new real[] {pts[i].x, pts[i].y, 1};
x[i] = -pts[i].x^2;
}
real[] coef = solve(M, x);
return parabola(changecoordsys(R, bqe(Rp, 1, 0, 0, coef[0], coef[1], coef[2])));
}

parabola parabola(point M1, point M2, point M3, point M4, point M5)
{/*Return the parabola passing through the five points.View index ordered by [name] - [type]*/
return parabola(bqe(M1, M2, M3, M4, M5));
}

hyperbola hyperbola(point C, real a, real b, real angle = 0)
{/*Return the hyperbola centered at 'C' with semimajor axis 'a' along C--C + dir(angle),
semiminor axis 'b' along the perpendicular.View index ordered by [name] - [type]*/
hyperbola oh;
coordsys R = C.coordsys;
angle += degrees(R.i);
real c = sqrt(a^2 + b^2);
point f1 = point(R, (locate(C) + rotate(angle) * (-c, 0))/R);
point f2 = point(R, (locate(C) + rotate(angle) * (c, 0))/R);
oh.init(f1, f2, a);
return oh;
}

hyperbola hyperbola(bqe bqe)
{/*Return the hyperbola a[0]x^2 + a[1]xy + a[2]y^2 + a[3]x + a[4]y + a[5]] = 0 (a[n] means bqe.a[n]).
http://mathworld.wolfram.com/Hyperbola.htmlView index ordered by [name] - [type]*/
bqe lbqe = changecoordsys(defaultcoordsys, bqe);
real a = lbqe.a[0], b = lbqe.a[1]/2, c = lbqe.a[2], d = lbqe.a[3]/2, f = lbqe.a[4]/2, g = lbqe.a[5];
string message = "hyperbola: the given equation is not an equation of a hyperbola.";
real delta = a * c * g + b * f * d + d * b * f - (b^2 * g + d^2 * c + f^2 * a);
if(abs(delta) < 10 * epsgeo || abs(b^2 - a * c) < 0) abort(message);
real dir = abs(b) < epsgeo ? 0 : 0.5 * acot(0.5 * (c - a)/b);
real cd = cos(dir), sd = sin(dir);
real ap = a * cd^2 - 2 * b * cd * sd + c * sd^2;
real cp = a * sd^2 + 2 * b * cd * sd + c * cd^2;
real dp = d * cd - f * sd;
real fp = d * sd + f * cd;
real gp = -g + dp^2/ap + fp^2/cp;
hyperbola op;
coordsys R = bqe.coordsys;
real j = b^2 - a * c;
point C = point(R, ((c * d - b * f)/j, (a * f - b * d)/j)/R);
real aa = gp/ap, bb = gp/cp;
real a = sqrt(abs(aa)), b = sqrt(abs(bb));
if(aa < 0) {dir -= pi/2; aa = a; a = b; b = aa;}
return hyperbola(C, a, b, degrees(-dir - angle(R.i)));
}

hyperbola hyperbola(point M1, point M2, point M3, point M4, point M5)
{/*Return the hyperbola passing through the five points (if possible).View index ordered by [name] - [type]*/
return hyperbola(bqe(M1, M2, M3, M4, M5));
}

hyperbola conj(hyperbola h)
{/*Conjugate.View index ordered by [name] - [type]*/
return hyperbola(h.C, h.b, h.a, 90 + h.angle);
}

circle circle(explicit point C, real r)
{/*Circle given by center and radius.View index ordered by [name] - [type]*/
circle oc = new circle;
oc.C = C;
oc.r = r;
if(!finite(r)) oc.l = line(C, C + vector(C.coordsys, (1, 0)));
return oc;
}

circle circle(point A, point B)
{/*Return the circle of diameter AB.View index ordered by [name] - [type]*/
real r;
circle oc;
real a = abs(A), b = abs(B);
if(finite(a) && finite(b)) {
oc = circle((A + B)/2, abs(A - B)/2);
} else {
oc.r = infinity;
if(finite(abs(A))) oc.l = line(A, A + unit(B));
else {
if(finite(abs(B))) oc.l = line(B, B + unit(A));
else if(finite(abs(A - B)/2)) oc = circle((A + B)/2, abs(A - B)/2); else
oc.l = line(A, B);
}
}
return oc;
}

circle circle(segment s)
{/*Return the circle of diameter 's'.View index ordered by [name] - [type]*/
return circle(s.A, s.B);
}

point circumcenter(point A, point B, point C)
{/*Return the circumcenter of triangle ABC.View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(A, B, C);
coordsys R = P[0].coordsys;
pair a = A, b = B, c = C;
pair mAB = (a + b)/2;
pair mAC = (a + c)/2;
pair pp = extension(mAB, rotate(90, mAB) * a, mAC, rotate(90, mAC) * c);
return point(R, pp/R);
}

circle circle(point A, point B, point C)
{/*Return the circumcircle of the triangle ABC.View index ordered by [name] - [type]*/
if(collinear(A - B, A - C)) {
circle oc;
oc.r = infinity;
oc.C = (A + B + C)/3;
oc.l = line(oc.C, oc.C == A ? B : A);
return oc;
}
point c = circumcenter(A, B, C);
return circle(c, abs(c - A));
}

circle circumcircle(point A, point B, point C)
{/*Return the circumcircle of the triangle ABC.View index ordered by [name] - [type]*/
return circle(A, B, C);
}

circle operator *(real x, explicit circle c)
{/*Multiply the radius of 'c'.View index ordered by [name] - [type]*/
return finite(c.r) ? circle(c.C, x * c.r) : c;
}
circle operator *(int x, explicit circle c)
{
return finite(c.r) ? circle(c.C, x * c.r) : c;
}

circle operator /(explicit circle c, real x)
{/*Divide the radius of 'c'View index ordered by [name] - [type]*/
return finite(c.r) ? circle(c.C, c.r/x) : c;
}
circle operator /(explicit circle c, int x)
{
return finite(c.r) ? circle(c.C, c.r/x) : c;
}

circle operator +(explicit circle c, explicit point M)
{/*Translation of 'c'.View index ordered by [name] - [type]*/
return circle(c.C + M, c.r);
}

circle operator -(explicit circle c, explicit point M)
{/*Translation of 'c'.View index ordered by [name] - [type]*/
return circle(c.C - M, c.r);
}

circle operator +(explicit circle c, pair m)
{/*Translation of 'c'.
'm' represent coordinates in the coordinate system where 'c' is defined.View index ordered by [name] - [type]*/
return circle(c.C + m, c.r);
}

circle operator -(explicit circle c, pair m)
{/*Translation of 'c'.
'm' represent coordinates in the coordinate system where 'c' is defined.View index ordered by [name] - [type]*/
return circle(c.C - m, c.r);
}

circle operator +(explicit circle c, vector m)
{/*Translation of 'c'.View index ordered by [name] - [type]*/
return circle(c.C + m, c.r);
}

circle operator -(explicit circle c, vector m)
{/*Translation of 'c'.View index ordered by [name] - [type]*/
return circle(c.C - m, c.r);
}

real operator ^(point M, explicit circle c)
{/*The power of 'M' with respect to the circle 'c'View index ordered by [name] - [type]*/
return xpart((abs(locate(M) - locate(c.C)), c.r)^2);
}

bool operator @(point M, explicit circle c)
{/*Return true iff 'M' is on the circle 'c'.View index ordered by [name] - [type]*/
return finite(c.r) ?
abs(abs(locate(M) - locate(c.C)) - abs(c.r)) <= 10 * epsgeo :
M @ c.l;
}

ellipse operator cast(circle c)
{/*View index ordered by [name] - [type]*/
return finite(c.r) ? ellipse(c.C, c.r, c.r, 0) : ellipse(c.l.A, c.l.B, infinity);
}

circle operator cast(ellipse el)
{/*View index ordered by [name] - [type]*/
circle oc;
bool infb = (!finite(el.a) || !finite(el.b));
if(!infb && abs(el.a - el.b) > epsgeo)
abort("Can not cast ellipse with different axis values to circle");
oc = circle(el.C, infb ? infinity : el.a);
oc.l = el.l.copy();
return oc;
}

ellipse operator cast(conic co)
{/*Cast a conic to an ellipse (can be a circle).View index ordered by [name] - [type]*/
if(degenerate(co) && co.e < 1) return ellipse(co.l[0].A, co.l[0].B, infinity);
ellipse oe;
if(co.e < 1) {
real a = co.p/(1 - co.e^2);
real c = co.e * a;
vector v = co.D.v;
if(!sameside(co.D.A + v, co.F, co.D)) v = -v;
point f2 = co.F + 2 * c * v;
f2 = changecoordsys(co.F.coordsys, f2);
oe = a == 0 ? ellipse(co.F, co.p, co.p, 0) : ellipse(co.F, f2, a);
} else
abort("casting: The conic section is not an ellipse.");
return oe;
}

parabola operator cast(conic co)
{/*Cast a conic to a parabola.View index ordered by [name] - [type]*/
parabola op;
if(abs(co.e - 1) > epsgeo) abort("casting: The conic section is not a parabola.");
op.init(co.F, co.D);
return op;
}

conic operator cast(parabola p)
{/*Cast a parabola to a conic section.View index ordered by [name] - [type]*/
return conic(p.F, p.D, 1);
}

hyperbola operator cast(conic co)
{/*Cast a conic section to an hyperbola.View index ordered by [name] - [type]*/
hyperbola oh;
if(co.e > 1) {
real a = co.p/(co.e^2 - 1);
real c = co.e * a;
vector v = co.D.v;
if(sameside(co.D.A + v, co.F, co.D)) v = -v;
point f2 = co.F + 2 * c * v;
f2 = changecoordsys(co.F.coordsys, f2);
oh = hyperbola(co.F, f2, a);
} else
abort("casting: The conic section is not an hyperbola.");
return oh;
}

conic operator cast(hyperbola h)
{/*Hyperbola to conic section.View index ordered by [name] - [type]*/
return conic(h.F1, h.D1, h.e);
}

conic operator cast(ellipse el)
{/*Ellipse to conic section.View index ordered by [name] - [type]*/
conic oc;
if(abs(el.c) > epsgeo) {
real x = el.a^2/el.c;
point O = (el.F1 + el.F2)/2;
point A = O + x * unit(el.F1 - el.F2);
oc = conic(el.F1, perpendicular(A, line(el.F1, el.F2)), el.e);
} else {//The ellipse is a circle
coordsys R = coordsys(el);
point M = el.F1 + point(R, R.polar(el.a, 0));
line l = line(rotate(90, M) * el.F1, M);
oc = conic(el.F1, l, 0);
}
if(degenerate(el)) {
oc.p = infinity;
oc.h = infinity;
oc.l = new line[]{el.l};
}
return oc;
}

conic operator cast(circle c)
{/*Circle to conic section.View index ordered by [name] - [type]*/
return (conic)((ellipse)c);
}

circle operator cast(conic c)
{/*Conic section to circle.View index ordered by [name] - [type]*/
ellipse el = (ellipse)c;
circle oc;
if(abs(el.a - el.b) < epsgeo) {
oc = circle(el.C, el.a);
if(degenerate(c)) oc.l = c.l[0];
}
else abort("Can not cast this conic to a circle");
return oc;
}

ellipse operator *(transform t, ellipse el)
{/*Provide transform * ellipse.View index ordered by [name] - [type]*/
if(!degenerate(el)) {
point[] ep;
for (int i = 0; i < 360; i += 72) {
ep.push(t * angpoint(el, i));
}
ellipse oe = ellipse(ep[0], ep[1], ep[2], ep[3], ep[4]);
if(angpoint(oe, 0) != ep[0]) return ellipse(oe.F2, oe.F1, oe.a);
return oe;
}
return ellipse(t * el.l.A, t * el.l.B, infinity);
}

parabola operator *(transform t, parabola p)
{/*Provide transform * parabola.View index ordered by [name] - [type]*/
point[] P;
P.push(t * angpoint(p, 45));
P.push(t * angpoint(p, -45));
P.push(t * angpoint(p, 180));
parabola op = parabola(P[0], P[1], P[2], t * p.D);
op.bmin = p.bmin;
op.bmax = p.bmax;

return op;
}

ellipse operator *(transform t, circle c)
{/*Provide transform * circle.
For example, 'circle C = scale(2) * circle' and 'ellipse E = xscale(2) * circle' are valid
but 'circle C = xscale(2) * circle' is invalid.View index ordered by [name] - [type]*/
return t * ((ellipse)c);
}

hyperbola operator *(transform t, hyperbola h)
{/*Provide transform * hyperbola.View index ordered by [name] - [type]*/
if (t == identity()) {
return h;
}

point[] ep;
for (int i = 90; i <= 270; i += 45) {
ep.push(t * angpoint(h, i));
}

hyperbola oe = hyperbola(ep[0], ep[1], ep[2], ep[3], ep[4]);
if(angpoint(oe, 90) != ep[0]) {
oe = hyperbola(oe.F2, oe.F1, oe.a);
}

oe.bmin = h.bmin;
oe.bmax = h.bmax;

return oe;
}

conic operator *(transform t, conic co)
{/*Provide transform * conic.View index ordered by [name] - [type]*/
if(co.e < 1) return (t * ((ellipse)co));
if(co.e == 1) return (t * ((parabola)co));
return (t * ((hyperbola)co));
}

ellipse operator *(real x, ellipse el)
{/*Identical but more efficient (rapid) than 'scale(x, el.C) * el'.View index ordered by [name] - [type]*/
return degenerate(el) ? el : ellipse(el.C, x * el.a, x * el.b, el.angle);
}

ellipse operator /(ellipse el, real x)
{/*Identical but more efficient (rapid) than 'scale(1/x, el.C) * el'.View index ordered by [name] - [type]*/
return degenerate(el) ? el : ellipse(el.C, el.a/x, el.b/x, el.angle);
}

path arcfromcenter(ellipse el, real angle1, real angle2,
bool direction=CCW,
int n=ellipsenodesnumber(el.a,el.b,angle1,angle2,direction))
{/*Return the path of the ellipse 'el' from angle1 to angle2 in degrees,
drawing in the given direction, with n nodes.
The angles are mesured relatively to the  axis (C,x-axis) where C is
the center of the ellipse.View index ordered by [name] - [type]*/
if(degenerate(el)) abort("arcfromcenter: can not convert degenerated ellipse to path.");
if (angle1 > angle2)
return reverse(arcfromcenter(el, angle2, angle1, !direction, n));

guide op;
coordsys Rp=coordsys(el);
if (n < 1) return op;

interpolate join = operator ..;
real stretch = max(el.a/el.b, el.b/el.a);

if (stretch > 10) {
n *= floor(stretch/5);
join = operator --;
}

real a2 = direction ? radians(angle2) : radians(angle1) + 2 * pi;
real step=(a2 - a1)/(n != 1 ? n-1 : 1);
real a, r;

for (int i=0; i < n; ++i) {
a = a1 + i * step;
r = el.b/sqrt(1 - (el.e * cos(a))^2);
op = join(op, Rp*Rp.polar(r, da + a));
}

return shift(el.C.x*Rp.i + el.C.y*Rp.j) * (direction ? op : reverse(op));
}

path arcfromcenter(hyperbola h, real angle1, real angle2,
int n = hyperbolanodesnumber(h, angle1, angle2),
bool direction = CCW)
{/*Return the path of the hyperbola 'h' from angle1 to angle2 in degrees,
drawing in the given direction, with n nodes.
The angles are mesured relatively to the axis (C, x-axis) where C is
the center of the hyperbola.View index ordered by [name] - [type]*/
guide op;
coordsys Rp = coordsys(h);
if (n < 1) return op;
if (angle1 > angle2) {
path g = reverse(arcfromcenter(h, angle2, angle1, n, !direction));
return g == nullpath ? g : reverse(g);
}
real a2 = direction ? radians(angle2) : radians(angle1) + 2 * pi;
real step = (a2 - a1)/(n != 1 ? n - 1 : 1);
real a, r;
typedef guide interpolate(... guide[]);
interpolate join = operator ..;
for (int i = 0; i < n; ++i) {
a = a1 + i * step;
r = (h.b * cos(a))^2 - (h.a * sin(a))^2;
if(r > epsgeo) {
r = sqrt(h.a^2 * h.b^2/r);
op = join(op, Rp * Rp.polar(r, a + da));
join = operator ..;
} else join = operator --;
}
return shift(h.C.x * Rp.i + h.C.y * Rp.j)*
(direction ? op : op == nullpath ? op : reverse(op));
}

path arcfromcenter(explicit conic co, real angle1, real angle2,
int n, bool direction = CCW)
{/*Use arcfromcenter(ellipse, ...) or arcfromcenter(hyperbola, ...) depending of
the eccentricity of 'co'.View index ordered by [name] - [type]*/
path g;
if(co.e < 1)
g = arcfromcenter((ellipse)co, angle1,
angle2, direction, n);
else if(co.e > 1)
g = arcfromcenter((hyperbola)co, angle1,
angle2, n, direction);
else abort("arcfromcenter: does not exist for a parabola.");
return g;
}

restricted polarconicroutine fromCenter = arcfromcenter;/*View index ordered by [name] - [type]*/

restricted polarconicroutine fromFocus = arcfromfocus;/*View index ordered by [name] - [type]*/

bqe equation(ellipse el)
{/*Return the coefficients of the equation of the ellipse in its coordinate system:
bqe.a[0] * x^2 + bqe.a[1] * x * y + bqe.a[2] * y^2 + bqe.a[3] * x + bqe.a[4] * y + bqe.a[5] = 0.
One can change the coordinate system of 'bqe' using the routine 'changecoordsys'.View index ordered by [name] - [type]*/
pair[] pts;
for (int i = 0; i < 360; i += 72)
pts.push(locate(angpoint(el, i)));

real[][] M;
real[] x;
for (int i = 0; i < 5; ++i) {
M[i] = new real[] {pts[i].x * pts[i].y, pts[i].y^2, pts[i].x, pts[i].y, 1};
x[i] = -pts[i].x^2;
}
real[] coef = solve(M, x);
bqe bqe = changecoordsys(coordsys(el),
bqe(defaultcoordsys,
1, coef[0], coef[1], coef[2], coef[3], coef[4]));
bqe.a = approximate(bqe.a);
return bqe;
}

bqe equation(parabola p)
{/*Return the coefficients of the equation of the parabola in its coordinate system.
bqe.a[0] * x^2 + bqe.a[1] * x * y + bqe.a[2] * y^2 + bqe.a[3] * x + bqe.a[4] * y + bqe.a[5] = 0
One can change the coordinate system of 'bqe' using the routine 'changecoordsys'.View index ordered by [name] - [type]*/
coordsys R = canonicalcartesiansystem(p);
parabola tp = changecoordsys(R, p);
point A = projection(tp.D) * point(R, (0, 0));
real a = abs(A);
return changecoordsys(coordsys(p),
bqe(R, 0, 0, 1, -4 * a, 0, 0));
}

bqe equation(hyperbola h)
{/*Return the coefficients of the equation of the hyperbola in its coordinate system.
bqe.a[0] * x^2 + bqe.a[1] * x * y + bqe.a[2] * y^2 + bqe.a[3] * x + bqe.a[4] * y + bqe.a[5] = 0
One can change the coordinate system of 'bqe' using the routine 'changecoordsys'.View index ordered by [name] - [type]*/
coordsys R = canonicalcartesiansystem(h);
return changecoordsys(coordsys(h),
bqe(R, 1/h.a^2, 0, -1/h.b^2, 0, 0, -1));
}

path operator cast(ellipse el)
{/*Cast ellipse to path.View index ordered by [name] - [type]*/
if(degenerate(el))
abort("Casting degenerated ellipse to path is not possible.");
int n = el.e == 0 ? circlenodesnumber(el.a) : ellipsenodesnumber(el.a, el.b);
return arcfromcenter(el, 0.0, 360, CCW, n)&cycle;
}

path operator cast(circle c)
{/*Cast circle to path.View index ordered by [name] - [type]*/
return (path)((ellipse)c);
}

real[] bangles(picture pic = currentpicture, parabola p)
{/*Return the array {ma, Ma} where 'ma' and 'Ma' are respectively
the smaller and the larger angles for which the parabola 'p' is included
in the bounding box of the picture 'pic'.View index ordered by [name] - [type]*/
pair bmin, bmax;
pair[] b;
if (p.bmin == p.bmax) {
bmin = pic.userMin();
bmax = pic.userMax();
} else {
bmin = p.bmin;bmax = p.bmax;
}
if(bmin.x == bmax.x || bmin.y == bmax.y || !finite(abs(bmin)) || !finite(abs(bmax)))
return new real[] {0, 0};
b[0] = bmin;
b[1] = (bmax.x, bmin.y);
b[2] = bmax;
b[3] = (bmin.x, bmax.y);
real[] eq = changecoordsys(defaultcoordsys, equation(p)).a;
pair[] inter;
for (int i = 0; i < 4; ++i) {
pair[] tmp = intersectionpoints(b[i], b[(i + 1)%4], eq);
for (int j = 0; j < tmp.length; ++j) {
if(dot(b[i]-tmp[j], b[(i + 1)%4]-tmp[j]) <= epsgeo)
inter.push(tmp[j]);
}
}
pair F = p.F, V = p.V;
real d = degrees(F - V);
real[] a = sequence(new real(int n){
return (360 - d + degrees(inter[n]-F))%360;
}, inter.length);
real ma = a.length != 0 ? min(a) : 0, Ma= a.length != 0 ? max(a) : 0;
return new real[] {ma, Ma};
}

real[][] bangles(picture pic = currentpicture, hyperbola h)
{/*Return the array {{ma1, Ma1}, {ma2, Ma2}} where 'maX' and 'MaX' are respectively
the smaller and the bigger angles (from h.FX) for which the hyperbola 'h' is included
in the bounding box of the picture 'pic'.View index ordered by [name] - [type]*/
pair bmin, bmax;
pair[] b;
if (h.bmin == h.bmax) {
bmin = pic.userMin();
bmax = pic.userMax();
} else {
bmin = h.bmin;bmax = h.bmax;
}
if(bmin.x == bmax.x || bmin.y == bmax.y || !finite(abs(bmin)) || !finite(abs(bmax)))
return new real[][] {{0, 0}, {0, 0}};
b[0] = bmin;
b[1] = (bmax.x, bmin.y);
b[2] = bmax;
b[3] = (bmin.x, bmax.y);
real[] eq = changecoordsys(defaultcoordsys, equation(h)).a;
pair[] inter0, inter1;
pair C = locate(h.C);
pair F1 = h.F1;
for (int i = 0; i < 4; ++i) {
pair[] tmp = intersectionpoints(b[i], b[(i + 1)%4], eq);
for (int j = 0; j < tmp.length; ++j) {
if(dot(b[i]-tmp[j], b[(i + 1)%4]-tmp[j]) <= epsgeo) {
if(dot(F1 - C, tmp[j]-C) > 0) inter0.push(tmp[j]);
else inter1.push(tmp[j]);
}
}
}
real d = degrees(F1 - C);
real[] ma, Ma;
pair[][] inter = new pair[][] {inter0, inter1};
for (int i = 0; i < 2; ++i) {
real[] a = sequence(new real(int n){
return (360 - d + degrees(inter[i][n]-F1))%360;
}, inter[i].length);
ma[i] = a.length != 0 ? min(a) : 0;
Ma[i] = a.length != 0 ? max(a) : 0;
}
return new real[][] {{ma[0], Ma[0]}, {ma[1], Ma[1]}};
}

path operator cast(parabola p)
{/*Cast parabola to path.
If possible, the returned path is restricted to the actual bounding box
of the current picture if the variables 'p.bmin' and 'p.bmax' are not set else
the bounding box of box(p.bmin, p.bmax) is used instead.View index ordered by [name] - [type]*/
real[] bangles = bangles(p);
int n = parabolanodesnumber(p, bangles[0], bangles[1]);
return arcfromfocus(p, bangles[0], bangles[1], n, CCW);
}

void draw(picture pic = currentpicture, Label L = "", circle c,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
margin margin = NoMargin, Label legend = "", marker marker = nomarker)
{/*View index ordered by [name] - [type]*/
if(degenerate(c)) draw(pic, L, c.l, align, p, arrow, legend, marker);
else draw(pic, L, (path)c, align, p, arrow, bar, margin, legend, marker);
}

void draw(picture pic = currentpicture, Label L = "", ellipse el,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
margin margin = NoMargin, Label legend = "", marker marker = nomarker)
{/*Draw the ellipse 'el' if it is not degenerated else draw 'el.l'.View index ordered by [name] - [type]*/
if(degenerate(el)) draw(pic, L, el.l, align, p, arrow, legend, marker);
else draw(pic, L, (path)el, align, p, arrow, bar, margin, legend, marker);
}

void draw(picture pic = currentpicture, Label L = "", parabola parabola,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
margin margin = NoMargin, Label legend = "", marker marker = nomarker)
{/*Draw the parabola 'p' on 'pic' without (if possible) altering the
size of picture pic.View index ordered by [name] - [type]*/
pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
// Reduce the bounds by the size of the pen and the margins.
m -= min(p); M -= max(p);
parabola.bmin = inverse(t) * m;
parabola.bmax = inverse(t) * M;
picture tmp;
path pp = t * ((path) (T * parabola));

if (pp != nullpath) {
draw(tmp, L, pp, align, p, arrow, bar, NoMargin, legend, marker);
}
}, true);

pair m = pic.userMin(), M = pic.userMax();
if(m != M) {
}
}

path operator cast(hyperbola h)
{/*Cast hyperbola to path.
If possible, the returned path is restricted to the actual bounding box
of the current picture unless the variables 'h.bmin' and 'h.bmax'
are set; in this case the bounding box of box(h.bmin, h.bmax) is used instead.
Only the branch on the side of 'h.F1' is considered.View index ordered by [name] - [type]*/
real[][] bangles = bangles(h);
int n = hyperbolanodesnumber(h, bangles[0][0], bangles[0][1]);
return arcfromfocus(h, bangles[0][0], bangles[0][1], n, CCW);
}

void draw(picture pic = currentpicture, Label L = "", hyperbola h,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
margin margin = NoMargin, Label legend = "", marker marker = nomarker)
{/*Draw the hyperbola 'h' on 'pic' without (if possible) altering the
size of the picture pic.View index ordered by [name] - [type]*/
pic.add(new void (frame f, transform t, transform T, pair m, pair M) {
// Reduce the bounds by the size of the pen and the margins.
m -= min(p); M -= max(p);
h.bmin = inverse(t) * m;
h.bmax = inverse(t) * M;
path hp;

picture tmp;
hp = t * ((path) (T * h));
if (hp != nullpath) {
draw(tmp, L, hp, align, p, arrow, bar, NoMargin, legend, marker);
}

hyperbola ht = hyperbola(h.F2, h.F1, h.a);
ht.bmin = h.bmin;
ht.bmax = h.bmax;

hp = t * ((path) (T * ht));
if (hp != nullpath) {
draw(tmp, "", hp, align, p, arrow, bar, NoMargin, marker);
}

}, true);

pair m = pic.userMin(), M = pic.userMax();
if(m != M)
}

void draw(picture pic = currentpicture, Label L = "", explicit conic co,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None,
margin margin = NoMargin, Label legend = "", marker marker = nomarker)
{/*Use one of the routine 'draw(ellipse, ...)',
'draw(parabola, ...)' or 'draw(hyperbola, ...)' depending of the value of eccentricity of 'co'.View index ordered by [name] - [type]*/
if(co.e == 0)
draw(pic, L, (circle)co, align, p, arrow, bar, margin, legend, marker);
else
if(co.e < 1) draw(pic, L, (ellipse)co, align, p, arrow, bar, margin, legend, marker);
else
if(co.e == 1) draw(pic, L, (parabola)co, align, p, arrow, bar, margin, legend, marker);
else
if(co.e > 1) draw(pic, L, (hyperbola)co, align, p, arrow, bar, margin, legend, marker);
else abort("draw: unknown conic.");
}

int conicnodesnumber(conic co, real angle1, real angle2, bool dir = CCW)
{/*Return the number of node to draw a conic arc.View index ordered by [name] - [type]*/
int oi;
if(co.e == 0) {
circle c = (circle)co;
oi = circlenodesnumber(c.r, angle1, angle2);
} else if(co.e < 1) {
ellipse el = (ellipse)co;
oi = ellipsenodesnumber(el.a, el.b, angle1, angle2, dir);
} else if(co.e == 1) {
parabola p = (parabola)co;
oi = parabolanodesnumber(p, angle1, angle2);
} else {
hyperbola h = (hyperbola)co;
oi = hyperbolanodesnumber(h, angle1, angle2);
}
return oi;
}

path operator cast(conic co)
{/*Cast conic section to path.View index ordered by [name] - [type]*/
if(co.e < 1) return (path)((ellipse)co);
if(co.e == 1) return (path)((parabola)co);
return (path)((hyperbola)co);
}

bqe equation(explicit conic co)
{/*Return the coefficients of the equation of conic section in its coordinate system:
bqe.a[0] * x^2 + bqe.a[1] * x * y + bqe.a[2] * y^2 + bqe.a[3] * x + bqe.a[4] * y + bqe.a[5] = 0.
One can change the coordinate system of 'bqe' using the routine 'changecoordsys'.View index ordered by [name] - [type]*/
bqe obqe;
if(co.e == 0)
obqe = equation((circle)co);
else
if(co.e < 1) obqe = equation((ellipse)co);
else
if(co.e == 1) obqe = equation((parabola)co);
else
if(co.e > 1) obqe = equation((hyperbola)co);
else abort("draw: unknown conic.");
return obqe;
}

string conictype(bqe bqe)
{/*Returned values are "ellipse" or "parabola" or "hyperbola"
depending of the conic section represented by 'bqe'.View index ordered by [name] - [type]*/
bqe lbqe = changecoordsys(defaultcoordsys, bqe);
string os = "degenerated";
real a = lbqe.a[0], b = lbqe.a[1]/2, c = lbqe.a[2], d = lbqe.a[3]/2, f = lbqe.a[4]/2, g = lbqe.a[5];
real delta = a * c * g + b * f * d + d * b * f - (b^2 * g + d^2 * c + f^2 * a);
if(abs(delta) < 10 * epsgeo) return os;
real J = a * c - b^2;
real I = a + c;
if(J > epsgeo) {
if(delta/I < -epsgeo);
os = "ellipse";
} else {
if(abs(J) < epsgeo) os = "parabola"; else os = "hyperbola";
}
return os;
}

conic conic(point M1, point M2, point M3, point M4, point M5)
{/*Return the conic passing through 'M1', 'M2', 'M3', 'M4' and 'M5' if the conic is not degenerated.View index ordered by [name] - [type]*/
bqe bqe = bqe(M1, M2, M3, M4, M5);
string ct = conictype(bqe);
if(ct == "degenerated") abort("conic: degenerated conic passing through five points.");
if(ct == "ellipse") return ellipse(bqe);
if(ct == "parabola") return parabola(bqe);
return hyperbola(bqe);
}

coordsys canonicalcartesiansystem(explicit conic co)
{/*Return the canonical cartesian system of the conic 'co'.View index ordered by [name] - [type]*/
if(co.e < 1) return canonicalcartesiansystem((ellipse)co);
else if(co.e == 1) return canonicalcartesiansystem((parabola)co);
return canonicalcartesiansystem((hyperbola)co);
}

bqe canonical(bqe bqe)
{/*Return the bivariate quadratic equation relative to the
canonical coordinate system of the conic section represented by 'bqe'.View index ordered by [name] - [type]*/
string type = conictype(bqe);
if(type == "") abort("canonical: the equation can not be performed.");
bqe obqe;
if(type == "ellipse") {
ellipse el = ellipse(bqe);
obqe = changecoordsys(canonicalcartesiansystem(el), equation(el));
} else {
if(type == "parabola") {
parabola p = parabola(bqe);
obqe = changecoordsys(canonicalcartesiansystem(p), equation(p));
} else {
hyperbola h = hyperbola(bqe);
obqe = changecoordsys(canonicalcartesiansystem(h), equation(h));
}
}
return obqe;
}

conic conic(bqe bqe)
{/*Return the conic section represented by the bivariate quartic equation 'bqe'.View index ordered by [name] - [type]*/
string type = conictype(bqe);
if(type == "") abort("canonical: the equation can not be performed.");
conic oc;
if(type == "ellipse") {
oc = ellipse(bqe);
} else {
if(type == "parabola") oc = parabola(bqe); else oc = hyperbola(bqe);
}
return oc;
}

real arclength(circle c)
{/*View index ordered by [name] - [type]*/
return c.r * 2 * pi;
}

real focusToCenter(ellipse el, real a)
{/*Return the angle relatively to the center of 'el' for the angle 'a'
given relatively to the focus of 'el'.View index ordered by [name] - [type]*/
pair p = point(fromFocus(el, a, a, 1, CCW), 0);
pair c = locate(el.C);
real d = degrees(p - c) - el.angle;
d = abs(d) < epsgeo ? 0 : d; // Avoid -1e-15
return d%(sgnd(a) * 360);
}

real centerToFocus(ellipse el, real a)
{/*Return the angle relatively to the focus of 'el' for the angle 'a'
given relatively to the center of 'el'.View index ordered by [name] - [type]*/
pair P = point(fromCenter(el, a, a, 1, CCW), 0);
pair F1 = locate(el.F1);
pair F2 = locate(el.F2);
real d = degrees(P - F1) - degrees(F2 - F1);
d = abs(d) < epsgeo ? 0 : d; // Avoid -1e-15
return d%(sgnd(a) * 360);
}

real arclength(ellipse el)
{/*View index ordered by [name] - [type]*/
return degenerate(el) ? infinity : 4 * el.a * elle(pi/2, el.e);
}

real arclength(ellipse el, real angle1, real angle2,
bool direction = CCW,
polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return the length of the arc of the ellipse between 'angle1'
and 'angle2'.
'angle1' and 'angle2' must be in the interval ]-360;+oo[ if polarconicroutine = fromFocus,
]-oo;+oo[ if polarconicroutine = fromCenter.View index ordered by [name] - [type]*/
if(degenerate(el)) return infinity;
if(angle1 > angle2) return arclength(el, angle2, angle1, !direction, polarconicroutine);
//   path g;int n = 1000;
//   if(el.e == 0) g = arcfromcenter(el, angle1, angle2, n, direction);
//   if(el.e != 1) g = polarconicroutine(el, angle1, angle2, n, direction);
//   write("with path = ", arclength(g));
if(polarconicroutine == fromFocus) {
//   dot(point(fromFocus(el, angle1, angle1, 1, CCW), 0), 2mm + blue);
//   dot(point(fromFocus(el, angle2, angle2, 1, CCW), 0), 2mm + blue);
//   write("fromfocus1 = ", angle1);
//   write("fromfocus2 = ", angle2);
real gle1 = focusToCenter(el, angle1);
real gle2 = focusToCenter(el, angle2);
if((gle1 - gle2) * (angle1 - angle2) > 0) {
angle1 = gle1; angle2 = gle2;
} else {
angle1 = gle2; angle2 = gle1;
}
//   dot(point(fromCenter(el, angle1, angle1, 1, CCW), 0), 1mm + red);
//   dot(point(fromCenter(el, angle2, angle2, 1, CCW), 0), 1mm + red);
//   write("fromcenter1 = ", angle1);
//   write("fromcenter2 = ", angle2);
}
if(angle1 < 0 || angle2 < 0) return arclength(el, 180 + angle1, 180 + angle2, direction, fromCenter);
real a1 = direction ? angle1 : angle2;
real a2 = direction ? angle2 : angle1 + 360;
real elleq = el.a * elle(pi/2, el.e);
real S(real a)
{//Return the arclength from 0 to the angle 'a' (in degrees)
// given form the center of the ellipse.
real gle = atan(el.a * Tan(a)/el.b)+
pi * (((a%90 == 0 && a != 0) ? floor(a/90) - 1 : floor(a/90)) -
((a%180 == 0) ? 0 : floor(a/180)) -
(a%360 == 0 ? floor(a/(360)) : 0));
/* // Uncomment to visualize the used branches
unitsize(2cm, 1cm);
import graph;

real xmin = 0, xmax = 3pi;

xlimits( xmin, xmax);
ylimits( 0, 10);
yaxis( "y" , LeftRight(), RightTicks(pTick=.8red, ptick = lightgrey, extend = true));
xaxis( "x - value", BottomTop(), Ticks(Label("$%.2f$", red), Step = pi/2, step = pi/4, pTick=.8red, ptick = lightgrey, extend = true));

real p2 = pi/2;
real f(real t)
{
return atan(0.6 * tan(t))+
pi * ((t%p2 == 0 && t != 0) ? floor(t/p2) - 1 : floor(t/p2)) -
((t%pi == 0) ? 0 : pi * floor(t/pi)) - (t%(2pi) == 0 ? pi * floor(t/(2 * pi)) : 0);
}

draw(graph(f, xmin, xmax, 100));
write(degrees(f(pi/2)));
write(degrees(f(pi)));
write(degrees(f(3pi/2)));
write(degrees(f(2pi)));
draw(graph(new real(real t){return t;}, xmin, xmax, 3));
*/
return elleq - el.a * elle(pi/2 - gle, el.e);
}
return S(a2) - S(a1);
}

real arclength(parabola p, real angle)
{/*Return the arclength from 180 to 'angle' given from focus in the
canonical coordinate system of 'p'.View index ordered by [name] - [type]*/
real a = p.a; /* In canonicalcartesiansystem(p) the equation of p
is x = y^2/(4a) */
// integrate(sqrt(1 + (x/(2 * a))^2), x);
real S(real t){return 0.5 * t * sqrt(1 + t^2/(4 * a^2)) + a * asinh(t/(2 * a));}
real R(real gle){return 2 * a/(1 - Cos(gle));}
real t = Sin(angle) * R(angle);
return S(t);
}

real arclength(parabola p, real angle1, real angle2)
{/*Return the arclength from 'angle1' to 'angle2' given from
focus in the canonical coordinate system of 'p'View index ordered by [name] - [type]*/
return arclength(p, angle1) - arclength(p, angle2);
}

real arclength(parabola p)
{/*Return the length of the arc of the parabola bounded to the bounding
box of the current picture.View index ordered by [name] - [type]*/
real[] b = bangles(p);
return arclength(p, b[0], b[1]);
}
// *........................CONICS.........................*
// *=======================================================*

// *=======================================================*
// *.......................ABSCISSA........................*

struct abscissa
{/*Provide abscissa structure on a curve used in the routine-like 'point(object, abscissa)'
where object can be 'line','segment','ellipse','circle','conic'...*/
real x;/*The abscissa value.*/
int system;/*0 = relativesystem; 1 = curvilinearsystem; 2 = angularsystem; 3 = nodesystem*/
polarconicroutine polarconicroutine = fromCenter;/*The routine used with angular system and two foci conic section.
Possible values are 'formCenter' and 'formFocus'.*/

abscissa copy()
{/*Return a copy of this abscissa.*/
abscissa oa = new abscissa;
oa.x = this.x;
oa.system = this.system;
oa.polarconicroutine = this.polarconicroutine;
return oa;
}
}/*View index ordered by [name] - [type]*/

restricted int relativesystem = 0, curvilinearsystem = 1, angularsystem = 2, nodesystem = 3;/*Constant used to set the abscissa system.View index ordered by [name] - [type]*/

abscissa operator cast(explicit position position)
{/*Cast position to abscissa.
If 'position' is relative, the abscissa is relative else it's a curvilinear abscissa.View index ordered by [name] - [type]*/
abscissa oarcc;
oarcc.x = position.position.x;
oarcc.system = position.relative ? relativesystem : curvilinearsystem;
return oarcc;
}

abscissa operator +(real x, explicit abscissa a)
{/*Provide 'real + abscissa'.
Return abscissa b so that b.x = a.x + x.
+(explicit abscissa, real), -(real, explicit abscissa) and -(explicit abscissa, real) are also defined.View index ordered by [name] - [type]*/
abscissa oa = a.copy();
oa.x = a.x + x;
return oa;
}

abscissa operator +(explicit abscissa a, real x)
{
return x + a;
}
abscissa operator +(int x, explicit abscissa a)
{
return ((real)x) + a;
}

abscissa operator -(explicit abscissa a)
{/*Return the abscissa b so that b.x = -a.x.View index ordered by [name] - [type]*/
abscissa oa;
oa.system = a.system;
oa.x = -a.x;
return oa;
}

abscissa operator -(real x, explicit abscissa a)
{
abscissa oa;
oa.system = a.system;
oa.x = x - a.x;
return oa;
}
abscissa operator -(explicit abscissa a, real x)
{
abscissa oa;
oa.system = a.system;
oa.x = a.x - x;
return oa;
}
abscissa operator -(int x, explicit abscissa a)
{
return ((real)x) - a;
}

abscissa operator *(real x, explicit abscissa a)
{/*Provide 'real * abscissa'.
Return abscissa b so that b.x = x * a.x.
*(explicit abscissa, real), /(real, explicit abscissa) and /(explicit abscissa, real) are also defined.View index ordered by [name] - [type]*/
abscissa oa;
oa.system = a.system;
oa.x = a.x * x;
return oa;
}
abscissa operator *(explicit abscissa a, real x)
{
return x * a;
}

abscissa operator /(real x, explicit abscissa a)
{
abscissa oa;
oa.system = a.system;
oa.x = x/a.x;
return oa;
}
abscissa operator /(explicit abscissa a, real x)
{
abscissa oa;
oa.system = a.system;
oa.x = a.x/x;
return oa;
}

abscissa operator /(int x, explicit abscissa a)
{
return ((real)x)/a;
}

abscissa relabscissa(real x)
{/*Return a relative abscissa.View index ordered by [name] - [type]*/
return (abscissa)(Relative(x));
}
abscissa relabscissa(int x)
{
return (abscissa)(Relative(x));
}

abscissa curabscissa(real x)
{/*Return a curvilinear abscissa.View index ordered by [name] - [type]*/
return (abscissa)((position)x);
}
abscissa curabscissa(int x)
{
return (abscissa)((position)x);
}

abscissa angabscissa(real x, polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return a angular abscissa.View index ordered by [name] - [type]*/
abscissa oarcc;
oarcc.x = x;
oarcc.polarconicroutine = polarconicroutine;
oarcc.system = angularsystem;
return oarcc;
}
abscissa angabscissa(int x, polarconicroutine polarconicroutine = currentpolarconicroutine)
{
return angabscissa((real)x, polarconicroutine);
}

abscissa nodabscissa(real x)
{/*Return an abscissa as time on the path.View index ordered by [name] - [type]*/
abscissa oarcc;
oarcc.x = x;
oarcc.system = nodesystem;
return oarcc;
}
abscissa nodabscissa(int x)
{
return nodabscissa((real)x);
}

abscissa operator cast(real x)
{/*Cast real to abscissa, precisely 'nodabscissa'.View index ordered by [name] - [type]*/
return nodabscissa(x);
}
abscissa operator cast(int x)
{
return nodabscissa((real)x);
}

point point(circle c, abscissa l)
{/*Return the point of 'c' which has the abscissa 'l.x'
according to the abscissa system 'l.system'.View index ordered by [name] - [type]*/
coordsys R = c.C.coordsys;
if (l.system == nodesystem)
return point(R, point((path)c, l.x)/R);
if (l.system == relativesystem)
return c.C + point(R, R.polar(c.r, 2 * pi * l.x));
if (l.system == curvilinearsystem)
return c.C + point(R, R.polar(c.r, l.x/c.r));
if (l.system == angularsystem)
return c.C + point(R, R.polar(c.r, radians(l.x)));
return (0, 0);
}

point point(ellipse el, abscissa l)
{/*Return the point of 'el' which has the abscissa 'l.x'
according to the abscissa system 'l.system'.View index ordered by [name] - [type]*/
if(el.e == 0) return point((circle)el, l);
coordsys R = coordsys(el);
if (l.system == nodesystem)
return point(R, point((path)el, l.x)/R);
if (l.system == relativesystem) {
return point(el, curabscissa((l.x%1) * arclength(el)));
}
if (l.system == curvilinearsystem) {
real a1 = 0, a2 = 360, cx = 0;
real aout = a1;
real x = abs(l.x)%arclength(el);
while (abs(cx - x) > epsgeo) {
aout = (a1 + a2)/2;
cx = arclength(el, 0, aout, CCW, fromCenter); //fromCenter is speeder
if(cx > x) a2 = (a1 + a2)/2; else a1 = (a1 + a2)/2;
}
path pel = fromCenter(el, sgn(l.x) * aout, sgn(l.x) * aout, 1, CCW);
return point(R, point(pel, 0)/R);
}
if (l.system == angularsystem) {
return point(R, point(l.polarconicroutine(el, l.x, l.x, 1, CCW), 0)/R);
}
return (0, 0);
}

point point(parabola p, abscissa l)
{/*Return the point of 'p' which has the abscissa 'l.x'
according to the abscissa system 'l.system'.View index ordered by [name] - [type]*/
coordsys R = coordsys(p);
if (l.system == nodesystem)
return point(R, point((path)p, l.x)/R);
if (l.system == relativesystem) {
real[] b = bangles(p);
real al = sgn(l.x) > 0 ? arclength(p, 180, b[1]) : arclength(p, 180, b[0]);
return point(p, curabscissa(abs(l.x) * al));
}
if (l.system == curvilinearsystem) {
real a1 = 1e-3, a2 = 360 - 1e-3, cx = infinity;
while (abs(cx - l.x) > epsgeo) {
cx = arclength(p, 180, (a1 + a2)/2);
if(cx > l.x) a2 = (a1 + a2)/2; else a1 = (a1 + a2)/2;
}
path pp = fromFocus(p, a1, a1, 1, CCW);
return point(R, point(pp, 0)/R);
}
if (l.system == angularsystem) {
return point(R, point(fromFocus(p, l.x, l.x, 1, CCW), 0)/R);
}
return (0, 0);
}

point point(hyperbola h, abscissa l)
{/*Return the point of 'h' which has the abscissa 'l.x'
according to the abscissa system 'l.system'.View index ordered by [name] - [type]*/
coordsys R = coordsys(h);
if (l.system == nodesystem)
return point(R, point((path)h, l.x)/R);
if (l.system == relativesystem) {
abort("point(hyperbola, relativeSystem) is not implemented...
Try relpoint((path)your_hyperbola, x);");
}
if (l.system == curvilinearsystem) {
abort("point(hyperbola, curvilinearSystem) is not implemented...");
}
if (l.system == angularsystem) {
return point(R, point(l.polarconicroutine(h, l.x, l.x, 1, CCW), 0)/R);
}
return (0, 0);
}

point point(explicit conic co, abscissa l)
{/*Return the curvilinear abscissa of 'M' on the conic 'co'.View index ordered by [name] - [type]*/
if(co.e == 0) return point((circle)co, l);
if(co.e < 1) return point((ellipse)co, l);
if(co.e == 1) return point((parabola)co, l);
return point((hyperbola)co, l);
}

point point(line l, abscissa x)
{/*Return the point of 'l' which has the abscissa 'l.x' according to the abscissa system 'l.system'.
Note that the origin is l.A, and point(l, relabscissa(x)) returns l.A + x.x * vector(l.B - l.A).View index ordered by [name] - [type]*/
coordsys R = l.A.coordsys;
if (x.system == nodesystem)
return l.A + (x.x < 0 ? 0 : x.x > 1 ? 1 : x.x) * vector(l.B - l.A);
if (x.system == relativesystem)
return l.A + x.x * vector(l.B - l.A);
if (x.system == curvilinearsystem)
return l.A + x.x * l.u;
if (x.system == angularsystem)
abort("point: what the meaning of angular abscissa on line ?.");
return (0, 0);
}

point point(line l, explicit real x)
{/*Return the point between node l.A and l.B (x <= 0 means l.A, x >=1 means l.B).View index ordered by [name] - [type]*/
return point(l, nodabscissa(x));
}
point point(line l, explicit int x)
{
return point(l, nodabscissa(x));
}

point point(explicit circle c, explicit real x)
{/*Return the point between node floor(x) and floor(x) + 1.View index ordered by [name] - [type]*/
return point(c, nodabscissa(x));
}
point point(explicit circle c, explicit int x)
{
return point(c, nodabscissa(x));
}

point point(explicit ellipse el, explicit real x)
{/*Return the point between node floor(x) and floor(x) + 1.View index ordered by [name] - [type]*/
return point(el, nodabscissa(x));
}
point point(explicit ellipse el, explicit int x)
{
return point(el, nodabscissa(x));
}

point point(explicit parabola p, explicit real x)
{/*Return the point between node floor(x) and floor(x) + 1.View index ordered by [name] - [type]*/
return point(p, nodabscissa(x));
}
point point(explicit parabola p, explicit int x)
{
return point(p, nodabscissa(x));
}

point point(explicit hyperbola h, explicit real x)
{/*Return the point between node floor(x) and floor(x) + 1.View index ordered by [name] - [type]*/
return point(h, nodabscissa(x));
}
point point(explicit hyperbola h, explicit int x)
{
return point(h, nodabscissa(x));
}

point point(explicit conic co, explicit real x)
{/*Return the point between node floor(x) and floor(x) + 1.View index ordered by [name] - [type]*/
point op;
if(co.e == 0) op = point((circle)co, nodabscissa(x));
else if(co.e < 1) op = point((ellipse)co, nodabscissa(x));
else if(co.e == 1) op = point((parabola)co, nodabscissa(x));
else op = point((hyperbola)co, nodabscissa(x));
return op;
}
point point(explicit conic co, explicit int x)
{
return point(co, (real)x);
}

point relpoint(line l, real x)
{/*Return the relative point of 'l' (0 means l.A,
1 means l.B, x means l.A + x * vector(l.B - l.A) ).View index ordered by [name] - [type]*/
return point(l, Relative(x));
}

point relpoint(explicit circle c, real x)
{/*Return the relative point of 'c' (0 means origin, 1 means end).
Origin is c.center + c.r * (1, 0).View index ordered by [name] - [type]*/
return point(c, Relative(x));
}

point relpoint(explicit ellipse el, real x)
{/*Return the relative point of 'el' (0 means origin, 1 means end).View index ordered by [name] - [type]*/
return point(el, Relative(x));
}

point relpoint(explicit parabola p, real x)
{/*Return the relative point of the path of the parabola
bounded by the bounding box of the current picture.
0 means origin, 1 means end, where the origin is the vertex of 'p'.View index ordered by [name] - [type]*/
return point(p, Relative(x));
}

point relpoint(explicit hyperbola h, real x)
{/*Not yet implemented... View index ordered by [name] - [type]*/
return point(h, Relative(x));
}

point relpoint(explicit conic co, explicit real x)
{/*Return the relative point of 'co' (0 means origin, 1 means end).View index ordered by [name] - [type]*/
point op;
if(co.e == 0) op = point((circle)co, Relative(x));
else if(co.e < 1) op = point((ellipse)co, Relative(x));
else if(co.e == 1) op = point((parabola)co, Relative(x));
else op = point((hyperbola)co, Relative(x));
return op;
}
point relpoint(explicit conic co, explicit int x)
{
return relpoint(co, (real)x);
}

point angpoint(explicit circle c, real x)
{/*Return the point of 'c' in the direction 'x' measured in degrees.View index ordered by [name] - [type]*/
return point(c, angabscissa(x));
}

point angpoint(explicit ellipse el, real x,
polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return the point of 'el' in the direction 'x'
measured in degrees according to 'polarconicroutine'.View index ordered by [name] - [type]*/
return el.e == 0 ? angpoint((circle) el, x) : point(el, angabscissa(x, polarconicroutine));
}

point angpoint(explicit parabola p, real x)
{/*Return the point of 'p' in the direction 'x' measured in degrees.View index ordered by [name] - [type]*/
return point(p, angabscissa(x));
}

point angpoint(explicit hyperbola h, real x,
polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return the point of 'h' in the direction 'x'
measured in degrees according to 'polarconicroutine'.View index ordered by [name] - [type]*/
return point(h, angabscissa(x, polarconicroutine));
}

point curpoint(line l, real x)
{/*Return the point of 'l' which has the curvilinear abscissa 'x'.
Origin is l.A.View index ordered by [name] - [type]*/
return point(l, curabscissa(x));
}

point curpoint(explicit circle c, real x)
{/*Return the point of 'c' which has the curvilinear abscissa 'x'.
Origin is c.center + c.r * (1, 0).View index ordered by [name] - [type]*/
return point(c, curabscissa(x));
}

point curpoint(explicit ellipse el, real x)
{/*Return the point of 'el' which has the curvilinear abscissa 'el'.View index ordered by [name] - [type]*/
return point(el, curabscissa(x));
}

point curpoint(explicit parabola p, real x)
{/*Return the point of 'p' which has the curvilinear abscissa 'x'.
Origin is the vertex of 'p'.View index ordered by [name] - [type]*/
return point(p, curabscissa(x));
}

point curpoint(conic co, real x)
{/*Return the point of 'co' which has the curvilinear abscissa 'x'.View index ordered by [name] - [type]*/
point op;
if(co.e == 0) op = point((circle)co, curabscissa(x));
else if(co.e < 1) op = point((ellipse)co, curabscissa(x));
else if(co.e == 1) op = point((parabola)co, curabscissa(x));
else op = point((hyperbola)co, curabscissa(x));
return op;
}

abscissa angabscissa(circle c, point M)
{/*Return the angular abscissa of 'M' on the circle 'c'.View index ordered by [name] - [type]*/
if(!(M @ c)) abort("angabscissa: the point is not on the circle.");
abscissa oa;
oa.system = angularsystem;
oa.x = degrees(M - c.C);
if(oa.x < 0) oa.x+=360;
return oa;
}

abscissa angabscissa(ellipse el, point M,
polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return the angular abscissa of 'M' on the ellipse 'el' according to 'polarconicroutine'.View index ordered by [name] - [type]*/
if(!(M @ el)) abort("angabscissa: the point is not on the ellipse.");
abscissa oa;
oa.system = angularsystem;
oa.polarconicroutine = polarconicroutine;
oa.x = polarconicroutine == fromCenter ? degrees(M - el.C) : degrees(M - el.F1);
oa.x -= el.angle;
if(oa.x < 0) oa.x += 360;
return oa;
}

abscissa angabscissa(hyperbola h, point M,
polarconicroutine polarconicroutine = currentpolarconicroutine)
{/*Return the angular abscissa of 'M' on the hyperbola 'h' according to 'polarconicroutine'.View index ordered by [name] - [type]*/
if(!(M @ h)) abort("angabscissa: the point is not on the hyperbola.");
abscissa oa;
oa.system = angularsystem;
oa.polarconicroutine = polarconicroutine;
oa.x = polarconicroutine == fromCenter ? degrees(M - h.C) : degrees(M - h.F1) + 180;
oa.x -= h.angle;
if(oa.x < 0) oa.x += 360;
return oa;
}

abscissa angabscissa(parabola p, point M)
{/*Return the angular abscissa of 'M' on the parabola 'p'.View index ordered by [name] - [type]*/
if(!(M @ p)) abort("angabscissa: the point is not on the parabola.");
abscissa oa;
oa.system = angularsystem;
oa.polarconicroutine = fromFocus;// Not used
oa.x = degrees(M - p.F);
oa.x -= p.angle;
if(oa.x < 0) oa.x += 360;
return oa;
}

abscissa angabscissa(explicit conic co, point M)
{/*Return the angular abscissa of 'M' on the conic 'co'.View index ordered by [name] - [type]*/
if(co.e == 0) return angabscissa((circle)co, M);
if(co.e < 1) return angabscissa((ellipse)co, M);
if(co.e == 1) return angabscissa((parabola)co, M);
return angabscissa((hyperbola)co, M);
}

abscissa curabscissa(line l, point M)
{/*Return the curvilinear abscissa of 'M' on the line 'l'.View index ordered by [name] - [type]*/
if(!(M @ extend(l))) abort("curabscissa: the point is not on the line.");
abscissa oa;
oa.system = curvilinearsystem;
oa.x = sgn(dot(M - l.A, l.B - l.A)) * abs(M - l.A);
return oa;
}

abscissa curabscissa(circle c, point M)
{/*Return the curvilinear abscissa of 'M' on the circle 'c'.View index ordered by [name] - [type]*/
if(!(M @ c)) abort("curabscissa: the point is not on the circle.");
abscissa oa;
oa.system = curvilinearsystem;
oa.x = pi * angabscissa(c, M).x * c.r/180;
return oa;
}

abscissa curabscissa(ellipse el, point M)
{/*Return the curvilinear abscissa of 'M' on the ellipse 'el'.View index ordered by [name] - [type]*/
if(!(M @ el)) abort("curabscissa: the point is not on the ellipse.");
abscissa oa;
oa.system = curvilinearsystem;
real a = angabscissa(el, M, fromCenter).x;
oa.x = arclength(el, 0, a, fromCenter);
oa.polarconicroutine = fromCenter;
return oa;
}

abscissa curabscissa(parabola p, point M)
{/*Return the curvilinear abscissa of 'M' on the parabola 'p'.View index ordered by [name] - [type]*/
if(!(M @ p)) abort("curabscissa: the point is not on the parabola.");
abscissa oa;
oa.system = curvilinearsystem;
real a = angabscissa(p, M).x;
oa.x = arclength(p, 180, a);
oa.polarconicroutine = fromFocus; // Not used.
return oa;
}

abscissa curabscissa(conic co, point M)
{/*Return the curvilinear abscissa of 'M' on the conic 'co'.View index ordered by [name] - [type]*/
if(co.e > 1) abort("curabscissa: not implemented for this hyperbola.");
if(co.e == 0) return curabscissa((circle)co, M);
if(co.e < 1) return curabscissa((ellipse)co, M);
return curabscissa((parabola)co, M);
}

abscissa nodabscissa(line l, point M)
{/*Return the node abscissa of 'M' on the line 'l'.View index ordered by [name] - [type]*/
if(!(M @ (segment)l)) abort("nodabscissa: the point is not on the segment.");
abscissa oa;
oa.system = nodesystem;
oa.x = abs(M - l.A)/abs(l.A - l.B);
return oa;
}

abscissa nodabscissa(circle c, point M)
{/*Return the node abscissa of 'M' on the circle 'c'.View index ordered by [name] - [type]*/
if(!(M @ c)) abort("nodabscissa: the point is not on the circle.");
abscissa oa;
oa.system = nodesystem;
oa.x = intersect((path)c, locate(M))[0];
return oa;
}

abscissa nodabscissa(ellipse el, point M)
{/*Return the node abscissa of 'M' on the ellipse 'el'.View index ordered by [name] - [type]*/
if(!(M @ el)) abort("nodabscissa: the point is not on the ellipse.");
abscissa oa;
oa.system = nodesystem;
oa.x = intersect((path)el, M)[0];
return oa;
}

abscissa nodabscissa(parabola p, point M)
{/*Return the node abscissa OF 'M' on the parabola 'p'.View index ordered by [name] - [type]*/
if(!(M @ p)) abort("nodabscissa: the point is not on the parabola.");
abscissa oa;
oa.system = nodesystem;
path pg = p;
real[] t = intersect(pg, M, 1e-5);
if(t.length == 0) abort("nodabscissa: the point is not on the path of the parabola.");
oa.x = t[0];
return oa;
}

abscissa nodabscissa(conic co, point M)
{/*Return the node abscissa of 'M' on the conic 'co'.View index ordered by [name] - [type]*/
if(co.e > 1) abort("nodabscissa: not implemented for hyperbola.");
if(co.e == 0) return nodabscissa((circle)co, M);
if(co.e < 1) return nodabscissa((ellipse)co, M);
return nodabscissa((parabola)co, M);
}

abscissa relabscissa(line l, point M)
{/*Return the relative abscissa of 'M' on the line 'l'.View index ordered by [name] - [type]*/
if(!(M @ extend(l))) abort("relabscissa: the point is not on the line.");
abscissa oa;
oa.system = relativesystem;
oa.x = sgn(dot(M - l.A, l.B - l.A)) * abs(M - l.A)/abs(l.A - l.B);
return oa;
}

abscissa relabscissa(circle c, point M)
{/*Return the relative abscissa of 'M' on the circle 'c'.View index ordered by [name] - [type]*/
if(!(M @ c)) abort("relabscissa: the point is not on the circle.");
abscissa oa;
oa.system = relativesystem;
oa.x = angabscissa(c, M).x/360;
return oa;
}

abscissa relabscissa(ellipse el, point M)
{/*Return the relative abscissa of 'M' on the ellipse 'el'.View index ordered by [name] - [type]*/
if(!(M @ el)) abort("relabscissa: the point is not on the ellipse.");
abscissa oa;
oa.system = relativesystem;
oa.x = curabscissa(el, M).x/arclength(el);
oa.polarconicroutine = fromFocus;
return oa;
}

abscissa relabscissa(conic co, point M)
{/*Return the relative abscissa of 'M'
on the conic 'co'.View index ordered by [name] - [type]*/
if(co.e > 1) abort("relabscissa: not implemented for hyperbola and parabola.");
if(co.e == 1) return relabscissa((parabola)co, M);
if(co.e == 0) return relabscissa((circle)co, M);
return relabscissa((ellipse)co, M);
}
// *.......................ABSCISSA........................*
// *=======================================================*

// *=======================================================*
// *.........................ARCS..........................*

struct arc {
/*Implement oriented ellipse (included circle) arcs.
All the calculus with this structure will be as exact as Asymptote can do.
For a full precision, you must not cast 'arc' to 'path' excepted for drawing routines.
*/
ellipse el;/*The support of the arc.*/
restricted real angle0 = 0;/*Internal use: rotating a circle does not modify the origin point,this variable stocks the eventual angle rotation. This value is not used for ellipses which are not circles.*/
restricted  real angle1, angle2;/*Values (in degrees) in ]-360, 360[.*/
bool direction = CCW;/*The arc will be drawn from 'angle1' to 'angle2' rotating in the direction 'direction'.*/
polarconicroutine polarconicroutine = currentpolarconicroutine;/*The routine to which the angles refer.
If 'el' is a circle 'fromCenter' is always used.*/

void setangles(real a0, real a1, real a2)
{/*Set the angles.*/
if (a1 < 0 && a2 < 0) {
a1 += 360;
a2 += 360;
}
this.angle0 = a0%(sgnd(a0) * 360);
this.angle1 = a1%(sgnd(a1) * 360);
this.angle2 = a2%(sgnd(2) * 360);
}

void init(ellipse el, real angle0 = 0, real angle1, real angle2,
polarconicroutine polarconicroutine,
bool direction = CCW)
{/*Constructor.*/
if(abs(angle1 - angle2) > 360) abort("arc: |angle1 - angle2| > 360.");
this.el = el;
this.setangles(angle0, angle1, angle2);
this.polarconicroutine = polarconicroutine;
this.direction = direction;
}

arc copy()
{/*Copy the arc.*/
arc oa = new arc;
oa.el = this.el;
oa.direction = this.direction;
oa.polarconicroutine = this.polarconicroutine;
oa.angle1 = this.angle1;
oa.angle2 = this.angle2;
oa.angle0 = this.angle0;
return oa;
}
}/*View index ordered by [name] - [type]*/

polarconicroutine polarconicroutine(conic co)
{/*Return the default routine used to draw a conic.View index ordered by [name] - [type]*/
if(co.e == 0) return fromCenter;
if(co.e == 1) return fromFocus;
return currentpolarconicroutine;
}

arc arc(ellipse el, real angle1, real angle2,
polarconicroutine polarconicroutine = polarconicroutine(el),
bool direction = CCW)
{/*Return the ellipse arc from 'angle1' to 'angle2' with respect to 'polarconicroutine' and rotating in the direction 'direction'.View index ordered by [name] - [type]*/
arc oa;
oa.init(el, 0, angle1, angle2, polarconicroutine, direction);
return oa;
}

arc complementary(arc a)
{/*Return the complementary of 'a'.View index ordered by [name] - [type]*/
arc oa;
oa.init(a.el, a.angle0, a.angle2, a.angle1, a.polarconicroutine, a.direction);
return oa;
}

arc reverse(arc a)
{/*Return arc 'a' oriented in reverse direction.View index ordered by [name] - [type]*/
arc oa;
oa.init(a.el, a.angle0, a.angle2, a.angle1, a.polarconicroutine, !a.direction);
return oa;
}

real degrees(arc a)
{/*Return the measure in degrees of the oriented arc 'a'.View index ordered by [name] - [type]*/
real or;
real da = a.angle2 - a.angle1;
if(a.direction) {
or = a.angle1 < a.angle2 ? da : 360 + da;
} else {
or = a.angle1 < a.angle2 ? -360 + da : da;
}
return or;
}

real angle(arc a)
{/*Return the measure in radians of the oriented arc 'a'.View index ordered by [name] - [type]*/
}

int arcnodesnumber(explicit arc a)
{/*Return the number of nodes to draw the arc 'a'.View index ordered by [name] - [type]*/
return ellipsenodesnumber(a.el.a, a.el.b, a.angle1, a.angle2, a.direction);
}

private path arctopath(arc a, int n)
{
if(a.el.e == 0) return arcfromcenter(a.el, a.angle0 + a.angle1, a.angle0 + a.angle2, a.direction, n);
if(a.el.e != 1) return a.polarconicroutine(a.el, a.angle1, a.angle2, n, a.direction);
return arcfromfocus(a.el, a.angle1, a.angle2, n, a.direction);
}

point angpoint(arc a, real angle)
{/*Return the point given by its angular position (in degrees) relative to the arc 'a'.
If 'angle > degrees(a)' or 'angle < 0' the returned point is on the extended arc.View index ordered by [name] - [type]*/
pair p;
if(a.el.e == 0) {
real gle = a.angle0 + a.angle1 + (a.direction ? angle : -angle);
p = point(arcfromcenter(a.el, gle, gle, CCW, 1), 0);
}
else {
real gle = a.angle1 + (a.direction ? angle : -angle);
p = point(a.polarconicroutine(a.el, gle, gle, 1, CCW), 0);
}
return point(coordsys(a.el), p/coordsys(a.el));
}

path operator cast(explicit arc a)
{/*Cast arc to path.View index ordered by [name] - [type]*/
return arctopath(a, arcnodesnumber(a));
}

guide operator cast(explicit arc a)
{/*Cast arc to guide.View index ordered by [name] - [type]*/
return arctopath(a, arcnodesnumber(a));
}

arc operator *(transform t, explicit arc a)
{/*Provide transform * arc.View index ordered by [name] - [type]*/
pair[] P, PP;
path g = arctopath(a, 3);
real a0, a1 = a.angle1, a2 = a.angle2, ap1, ap2;
bool dir = a.direction;
P[0] = t * point(g, 0);
P[1] = t * point(g, 2);
ellipse el = t * a.el;
arc oa;
a0 = (a.angle0 + angle(shiftless(t)))%360;
pair C;
if(a.polarconicroutine == fromCenter) C = el.C; else C = el.F1;
real d = abs(locate(el.F2 - el.F1)) > epsgeo ?
degrees(locate(el.F2 - el.F1)) : a0 + degrees(el.C.coordsys.i);
ap1 = (degrees(P[0]-C, false) - d)%360;
ap2 = (degrees(P[1]-C, false) - d)%360;
oa.init(el, a0, ap1, ap2, a.polarconicroutine, dir);
g = arctopath(oa, 3);
PP[0] = point(g, 0);
PP[1] = point(g, 2);
if((a1 - a2) * (ap1 - ap2) < 0) {// Handle reflection.
dir=!a.direction;
oa.init(el, a0, ap1, ap2, a.polarconicroutine, dir);
}
return oa;
}

arc operator *(real x, explicit arc a)
{/*Provide real * arc.
Return the arc subtracting and adding '(x - 1) * degrees(a)/2' to 'a.angle1' and 'a.angle2' respectively.View index ordered by [name] - [type]*/
real a1, a2, gle;
gle = (x - 1) * degrees(a)/2;
a1 = a.angle1 - gle;
a2 = a.angle2 + gle;
arc oa;
oa.init(a.el, a.angle0, a1, a2, a.polarconicroutine, a.direction);
return oa;
}
arc operator *(int x, explicit arc a){return (real)x * a;}

arc operator /(explicit arc a, real x)
{/*Provide arc/real.
Return the arc subtracting and adding '(1/x - 1) * degrees(a)/2' to 'a.angle1' and 'a.angle2' respectively.View index ordered by [name] - [type]*/
return (1/x) * a;
}

arc operator +(explicit arc a, point M)
{/*Provide arc + point.
Return shifted arc.
'operator +(explicit arc, point)', 'operator +(explicit arc, vector)' and 'operator -(explicit arc, vector)' are also defined.View index ordered by [name] - [type]*/
return shift(M) * a;
}
arc operator -(explicit arc a, point M){return a + (-M);}
arc operator +(explicit arc a, vector v){return shift(locate(v)) * a;}
arc operator -(explicit arc a, vector v){return a + (-v);}

bool operator @(point M, arc a)
{/*Return true iff 'M' is on the arc 'a'.View index ordered by [name] - [type]*/
if (!(M @ a.el)) return false;
coordsys R = defaultcoordsys;
path ap = arctopath(a, 3);
line l = line(point(R, point(ap, 0)), point(R, point(ap, 2)));
return sameside(M, point(R, point(ap, 1)), l);
}

void draw(picture pic = currentpicture, Label L = "", arc a,
align align = NoAlign, pen p = currentpen,
arrowbar arrow = None, arrowbar bar = None, margin margin = NoMargin,
Label legend = "", marker marker = nomarker)
{/*Draw 'arc' adding the pen returned by 'addpenarc(p)' to the pen 'p'.
Look at addpenarcView index ordered by [name] - [type]*/
draw(pic, L, (path)a, align, addpenarc(p), arrow, bar, margin, legend, marker);
}

real arclength(arc a)
{/*The arc length of 'a'.View index ordered by [name] - [type]*/
return arclength(a.el, a.angle1, a.angle2, a.direction, a.polarconicroutine);
}

private point ppoint(arc a, real x)
{// Return the point of the arc proportionally to its length.
point oP;
if(a.el.e == 0) { // Case of circle.
oP = angpoint(a, x * abs(degrees(a)));
} else { // Ellipse and not circle.
if(!a.direction) {
transform t = reflect(line(a.el.F1, a.el.F2));
return t * ppoint(t * a, x);
}

real angle1 = a.angle1, angle2 = a.angle2;
if(a.polarconicroutine == fromFocus) {
//       dot(point(fromFocus(a.el, angle1, angle1, 1, CCW), 0), 2mm + blue);
//       dot(point(fromFocus(a.el, angle2, angle2, 1, CCW), 0), 2mm + blue);
//       write("fromfocus1 = ", angle1);
//       write("fromfocus2 = ", angle2);
real gle1 = focusToCenter(a.el, angle1);
real gle2 = focusToCenter(a.el, angle2);
if((gle1 - gle2) * (angle1 - angle2) > 0) {
angle1 = gle1; angle2 = gle2;
} else {
angle1 = gle2; angle2 = gle1;
}
//       write("fromcenter1 = ", angle1);
//       write("fromcenter2 = ", angle2);
//       dot(point(fromCenter(a.el, angle1, angle1, 1, CCW), 0), 1mm + red);
//       dot(point(fromCenter(a.el, angle2, angle2, 1, CCW), 0), 1mm + red);
}

if(angle1 > angle2) {
arc ta = a.copy();
ta.polarconicroutine = fromCenter;
ta.setangles(a0 = a.angle0, a1 = angle1 - 360, a2 = angle2);
return ppoint(ta, x);
}
ellipse co = a.el;
real gle, a1, a2, cx = 0;
bool direction;
if(x >= 0) {
a1 = angle1;
a2 = a1 + 360;
direction = CCW;
} else {
a1 = angle1 - 360;
a2 = a1 - 360;
direction = CW;
}
gle = a1;
real L = arclength(co, angle1, angle2, a.direction, fromCenter);
real tx = L * abs(x)%arclength(co);
real aout = a1;
while(abs(cx - tx) > epsgeo) {
aout = (a1 + a2)/2;
cx = abs(arclength(co, gle, aout, direction, fromCenter));
if(cx > tx) a2 = (a1 + a2)/2 ; else a1 = (a1 + a2)/2;
}
pair p = point(arcfromcenter(co, aout, aout, CCW, 1), 0);
oP = point(coordsys(co), p/coordsys(co));
}
return oP;
}

point point(arc a, abscissa l)
{/*Return the point of 'a' which has the abscissa 'l.x'
according to the abscissa system 'l.system'.
Note that 'a.polarconicroutine' is used instead of 'l.polarconicroutine'.
Look at struct abscissaView index ordered by [name] - [type]*/
real posx;
arc ta = a.copy();
ellipse co = a.el;
if (l.system == relativesystem) {
posx = l.x;
} else
if (l.system == curvilinearsystem) {
real tl;
if(co.e == 0) {
tl = curabscissa(a.el, angpoint(a.el, a.angle0 + a.angle1)).x;
return curpoint(a.el, tl + (a.direction ? l.x : -l.x));
} else {
tl = curabscissa(a.el, angpoint(a.el, a.angle1, a.polarconicroutine)).x;
return curpoint(a.el, tl + (a.direction ? l.x : -l.x));
}
} else
if (l.system == nodesystem) {
coordsys R = coordsys(co);
return point(R, point((path)a, l.x)/R);
} else
if (l.system == angularsystem) {
return angpoint(a, l.x);
} else abort("point: bad abscissa system.");
return ppoint(ta, posx);
}

point point(arc a, real x)
{/*Return the point between node floor(t) and floor(t) + 1.View index ordered by [name] - [type]*/
return point(a, nodabscissa(x));
}
pair point(explicit arc a, int x)
{
return point(a, nodabscissa(x));
}

point relpoint(arc a, real x)
{/*Return the relative point of 'a'.
If x > 1 or x < 0, the returned point is on the extended arc.View index ordered by [name] - [type]*/
return point(a, relabscissa(x));
}

point curpoint(arc a, real x)
{/*Return the point of 'a' which has the curvilinear abscissa 'x'.
If x < 0 or x > arclength(a), the returned point is on the extended arc.View index ordered by [name] - [type]*/
return point(a, curabscissa(x));
}

abscissa angabscissa(arc a, point M)
{/*Return the angular abscissa of 'M' according to the arc 'a'.View index ordered by [name] - [type]*/
if(!(M @ a.el))
abort("angabscissa: the point is not on the extended arc.");
abscissa oa;
oa.system = angularsystem;
oa.polarconicroutine = a.polarconicroutine;
real am = angabscissa(a.el, M, a.polarconicroutine).x;
oa.x = (am - a.angle1 - (a.el.e == 0 ? a.angle0 : 0))%360;
oa.x = a.direction ? oa.x : 360 - oa.x;
return oa;
}

abscissa curabscissa(arc a, point M)
{/*Return the curvilinear abscissa according to the arc 'a'.View index ordered by [name] - [type]*/
ellipse el = a.el;
if(!(M @ el))
abort("angabscissa: the point is not on the extended arc.");
abscissa oa;
oa.system = curvilinearsystem;
real xm = curabscissa(el, M).x;
real a0 = el.e == 0 ? a.angle0 : 0;
real am = curabscissa(el, angpoint(el, a.angle1 + a0, a.polarconicroutine)).x;
real l = arclength(el);
oa.x = (xm - am)%l;
oa.x = a.direction ? oa.x : l - oa.x;
return oa;
}

abscissa nodabscissa(arc a, point M)
{/*Return the node abscissa according to the arc 'a'.View index ordered by [name] - [type]*/
if(!(M @ a))
abort("nodabscissa: the point is not on the arc.");
abscissa oa;
oa.system = nodesystem;
oa.x = intersect((path)a, M)[0];
return oa;
}

abscissa relabscissa(arc a, point M)
{/*Return the relative abscissa according to the arc 'a'.View index ordered by [name] - [type]*/
ellipse el = a.el;
if(!( M @ el))
abort("relabscissa: the point is not on the prolonged arc.");
abscissa oa;
oa.system = relativesystem;
oa.x = curabscissa(a, M).x/arclength(a);
return oa;
}

void markarc(picture pic = currentpicture,
Label L = "",
int n = 1, real radius = 0, real space = 0,
arc a,
pen sectorpen = currentpen,
pen markpen = sectorpen,
margin margin = NoMargin,
arrowbar arrow = None,
marker marker = nomarker)
{/*View index ordered by [name] - [type]*/
real Da = degrees(a);
pair p1 = point(a, 0);
pair p2 = relpoint(a, 1);
pair c = a.polarconicroutine == fromCenter ? locate(a.el.C) : locate(a.el.F1);
draw(c--p1^^c--p2, sectorpen);
markangle(pic = pic, L = L, n = n, radius = radius, space = space,
A = p1, O = c, B = p2,
arrow = arrow, p = markpen, margin = margin,
marker = marker);
}
// *.........................ARCS..........................*
// *=======================================================*

// *=======================================================*
// *........................MASSES.........................*

struct mass {/**/
point M;/**/
real m;/**/
}/*View index ordered by [name] - [type]*/

mass mass(point M, real m)
{/*Constructor of mass point.View index ordered by [name] - [type]*/
mass om;
om.M = M;
om.m = m;
return om;
}

point operator cast(mass m)
{/*Cast mass point to point.View index ordered by [name] - [type]*/
point op;
op = m.M;
op.m = m.m;
return op;
}

point point(explicit mass m){return m;}/*Cast
'm' to pointView index ordered by [name] - [type]*/

mass operator cast(point M)
{/*Cast point to mass point.View index ordered by [name] - [type]*/
mass om;
om.M = M;
om.m = M.m;
return om;
}

mass mass(explicit point P)
{/*Cast 'P' to mass.View index ordered by [name] - [type]*/
return mass(P, P.m);
}

point[] operator cast(mass[] m)
{/*Cast mass[] to point[].View index ordered by [name] - [type]*/
point[] op;
for(mass am : m) op.push(point(am));
return op;
}

mass[] operator cast(point[] P)
{/*Cast point[] to mass[].View index ordered by [name] - [type]*/
mass[] om;
for(point op : P) om.push(mass(op));
return om;
}

mass mass(coordsys R, explicit pair p, real m)
{/*Return the mass which has coordinates
'p' with respect to 'R' and weight 'm'.View index ordered by [name] - [type]*/
return point(R, p, m);// Using casting.
}

mass operator cast(pair m){return mass((point)m, 1);}/*Cast pair to mass point.View index ordered by [name] - [type]*/

path operator cast(mass M){return M.M;}/*Cast mass point to path.View index ordered by [name] - [type]*/

guide operator cast(mass M){return M.M;}/*Cast mass to guide.View index ordered by [name] - [type]*/

mass operator +(mass M1, mass M2)
{/*Provide mass + mass.
mass - mass is also defined.View index ordered by [name] - [type]*/
return mass(M1.M + M2.M, M1.m + M2.m);
}
mass operator -(mass M1, mass M2)
{
return mass(M1.M - M2.M, M1.m - M2.m);
}

mass operator *(real x, explicit mass M)
{/*Provide real * mass.
The resulted mass is the mass of 'M' multiplied by 'x' .
mass/real, mass + real and mass - real are also defined.View index ordered by [name] - [type]*/
return mass(M.M, x * M.m);
}
mass operator *(int x, explicit mass M){return mass(M.M, x * M.m);}
mass operator /(explicit mass M, real x){return mass(M.M, M.m/x);}
mass operator /(explicit mass M, int x){return mass(M.M, M.m/x);}
mass operator +(explicit mass M, real x){return mass(M.M, M.m + x);}
mass operator +(explicit mass M, int x){return mass(M.M, M.m + x);}
mass operator -(explicit mass M, real x){return mass(M.M, M.m - x);}
mass operator -(explicit mass M, int x){return mass(M.M, M.m - x);}

mass operator *(transform t, mass M)
{/*Provide transform * mass.View index ordered by [name] - [type]*/
return mass(t * M.M, M.m);
}

mass masscenter(... mass[] M)
{/*Return the center of the masses 'M'.View index ordered by [name] - [type]*/
point[] P;
for (int i = 0; i < M.length; ++i)
P.push(M[i].M);
P = standardizecoordsys(currentcoordsys, true ... P);
real m = M[0].m;
point oM = M[0].m * P[0];
for (int i = 1; i < M.length; ++i) {
oM += M[i].m * P[i];
m += M[i].m;
}
if (m == 0) abort("masscenter: the sum of masses is null.");
return mass(oM/m, m);
}

string massformat(string format = defaultmassformat,
string s, mass M)
{/*Return the string formated by 'format' with the mass value.
In the parameter 'format', %L will be replaced by 's'.
Look at defaultmassformat.View index ordered by [name] - [type]*/
return format == "" ? s :
format(replace(format, "%L", replace(s, "", "")), M.m); } void label(picture pic = currentpicture, Label L, explicit mass M, align align = NoAlign, string format = defaultmassformat, pen p = nullpen, filltype filltype = NoFill) {/*Draw label returned by massformat(format, L, M) at coordinates of M. Look at massformat(string, string, mass).View index ordered by [name] - [type]*/ Label lL = L.copy(); lL.s = massformat(format, lL.s, M); Label L = Label(lL, M.M, align, p, filltype); add(pic, L); } void dot(picture pic = currentpicture, Label L, explicit mass M, align align = NoAlign, string format = defaultmassformat, pen p = currentpen) {/*Draw a dot with label 'L' as label(picture, Label, explicit mass, align, string, pen, filltype) does. Look at label(picture, Label, mass, align, string, pen, filltype).View index ordered by [name] - [type]*/ Label lL = L.copy(); lL.s = massformat(format, lL.s, M); lL.position(locate(M.M)); lL.align(align, E); lL.p(p); dot(pic, M.M, p); add(pic, lL); } // *........................MASSES.........................* // *=======================================================* // *=======================================================* // *.......................TRIANGLES.......................* point orthocentercenter(point A, point B, point C) {/*Return the orthocenter of the triangle ABC.View index ordered by [name] - [type]*/ point[] P = standardizecoordsys(A, B, C); coordsys R = P[0].coordsys; pair pp = extension(A, projection(P[1], P[2]) * P[0], B, projection(P[0], P[2]) * P[1]); return point(R, pp/R); } point centroid(point A, point B, point C) {/*Return the centroid of the triangle ABC.View index ordered by [name] - [type]*/ return (A + B + C)/3; } point incenter(point A, point B, point C) {/*Return the center of the incircle of the triangle ABC.View index ordered by [name] - [type]*/ point[] P = standardizecoordsys(A, B, C); coordsys R = P[0].coordsys; pair a = A, b = B, c = C; pair pp = extension(a, a + dir(a--b, a--c), b, b + dir(b--a, b--c)); return point(R, pp/R); } real inradius(point A, point B, point C) {/*Return the radius of the incircle of the triangle ABC.View index ordered by [name] - [type]*/ point IC = incenter(A, B, C); return abs(IC - projection(A, B) * IC); } circle incircle(point A, point B, point C) {/*Return the incircle of the triangle ABC.View index ordered by [name] - [type]*/ point IC = incenter(A, B, C); return circle(IC, abs(IC - projection(A, B) * IC)); } point excenter(point A, point B, point C) {/*Return the center of the excircle of the triangle tangent with (AB).View index ordered by [name] - [type]*/ point[] P = standardizecoordsys(A, B, C); coordsys R = P[0].coordsys; pair a = A, b = B, c = C; pair pp = extension(a, a + rotate(90) * dir(a--b, a--c), b, b + rotate(90) * dir(b--a, b--c)); return point(R, pp/R); } real exradius(point A, point B, point C) {/*Return the radius of the excircle of the triangle ABC with (AB).View index ordered by [name] - [type]*/ point EC = excenter(A, B, C); return abs(EC - projection(A, B) * EC); } circle excircle(point A, point B, point C) {/*Return the excircle of the triangle ABC tangent with (AB).View index ordered by [name] - [type]*/ point center = excenter(A, B, C); real radius = abs(center - projection(B, C) * center); return circle(center, radius); } private int[] numarray = {1, 2, 3}; numarray.cyclic = true; struct triangle {/**/ struct vertex {/*Structure used to communicate the vertex of a triangle.*/ int n;/*1 means VA,2 means VB,3 means VC,4 means VA etc...*/ triangle t;/*The triangle to which the vertex refers.*/ }/*View index ordered by [name] - [type]*/ restricted point A, B, C;/*The vertices of the triangle (as point).*/ restricted vertex VA, VB, VC;/*The vertices of the triangle (as vertex). Note that the vertex structure contains the triangle to wish it refers.*/ VA.n = 1;VB.n = 2;VC.n = 3; vertex vertex(int n) {/*Return numbered vertex. 'n' is 1 means VA, 2 means VB, 3 means VC, 4 means VA etc...*/ n = numarray[n - 1]; if(n == 1) return VA; else if(n == 2) return VB; return VC; } point point(int n) {/*Return numbered point. n is 1 means A, 2 means B, 3 means C, 4 means A etc...*/ n = numarray[n - 1]; if(n == 1) return A; else if(n == 2) return B; return C; } void init(point A, point B, point C) {/*Constructor.*/ point[] P = standardizecoordsys(A, B, C); this.A = P[0]; this.B = P[1]; this.C = P[2]; VA.t = this; VB.t = this; VC.t = this; } void operator init(point A, point B, point C) {/*For backward compatibility. Provide the routine 'triangle(point A, point B, point C)'.*/ this.init(A, B, C); } void operator init(real b, real alpha, real c, real angle = 0, point A = (0, 0)) {/*For backward compatibility. Provide the routine 'triangle(real b, real alpha, real c, real angle = 0, point A = (0, 0)) which returns the triangle ABC rotated by 'angle' (in degrees) and where b = AC, degrees(A) = alpha, AB = c.*/ coordsys R = A.coordsys; this.init(A, A + R.polar(c, radians(angle)), A + R.polar(b, radians(angle + alpha))); } real a() {/*Return the length BC. b() and c() are also defined and return the length AC and AB respectively.*/ return length(C - B); } real b() {return length(A - C);} real c() {return length(B - A);} private real det(pair a, pair b) {return a.x * b.y - a.y * b.x;} real area() {/**/ pair a = locate(A), b = locate(B), c = locate(C); return 0.5 * abs(det(a, b) + det(b, c) + det(c, a)); } real alpha() {/*Return the measure (in degrees) of the angle A. beta() and gamma() are also defined and return the measure of the angles B and C respectively.*/ return degrees(acos((b()^2 + c()^2 - a()^2)/(2b() * c()))); } real beta() {return degrees(acos((c()^2 + a()^2 - b()^2)/(2c() * a())));} real gamma() {return degrees(acos((a()^2 + b()^2 - c()^2)/(2a() * b())));} path Path() {/*The path of the triangle.*/ return A--C--B--cycle; } struct side {/*Structure used to communicate the side of a triangle.*/ int n;/*1 or 0 means [AB],-1 means [BA],2 means [BC],-2 means [CB] etc.*/ triangle t;/*The triangle to which the side refers.*/ }/*View index ordered by [name] - [type]*/ side AB;/*For the routines using the structure 'side', triangle.AB means 'side AB'. BA, AC, CA etc are also defined.*/ AB.n = 1; AB.t = this; side BA; BA.n = -1; BA.t = this; side BC; BC.n = 2; BC.t = this; side CB; CB.n = -2; CB.t = this; side CA; CA.n = 3; CA.t = this; side AC; AC.n = -3; AC.t = this; side side(int n) {/*Return numbered side. n is 1 means AB, -1 means BA, 2 means BC, -2 means CB, etc.*/ if(n == 0) abort('Invalid side number.'); int an = numarray[abs(n)-1]; if(an == 1) return n > 0 ? AB : BA; else if(an == 2) return n > 0 ? BC : CB; return n > 0 ? CA : AC; } line line(int n) {/*Return the numbered line.*/ if(n == 0) abort('Invalid line number.'); int an = numarray[abs(n)-1]; if(an == 1) return n > 0 ? line(A, B) : line(B, A); else if(an == 2) return n > 0 ? line(B, C) : line(C, B); return n > 0 ? line(C, A) : line(A, C); } }/*View index ordered by [name] - [type]*/ from triangle unravel side; // The structure 'side' is now available outside the triangle structure. from triangle unravel vertex; // The structure 'vertex' is now available outside the triangle structure. triangle[] operator ^^(triangle[] t1, triangle t2) { triangle[] T; for (int i = 0; i < t1.length; ++i) T.push(t1[i]); T.push(t2); return T; } triangle[] operator ^^(... triangle[] t) { triangle[] T; for (int i = 0; i < t.length; ++i) { T.push(t[i]); } return T; } line operator cast(side side) {/*Cast side to (infinite) line. Most routine with line parameters works with side parameters. One can use the code 'segment(a_side)' to obtain a line segment.View index ordered by [name] - [type]*/ triangle t = side.t; return t.line(side.n); } line line(explicit side side) {/*Return 'side' as line.View index ordered by [name] - [type]*/ return (line)side; } segment segment(explicit side side) {/*Return 'side' as segment.View index ordered by [name] - [type]*/ return (segment)(line)side; } point operator cast(vertex V) {/*Cast vertex to point. Most routine with point parameters works with vertex parameters.View index ordered by [name] - [type]*/ return V.t.point(V.n); } point point(explicit vertex V) {/*Return the point corresponding to the vertex 'V'.View index ordered by [name] - [type]*/ return (point)V; } side opposite(vertex V) {/*Return the opposite side of vertex 'V'.View index ordered by [name] - [type]*/ return V.t.side(numarray[abs(V.n)]); } vertex opposite(side side) {/*Return the opposite vertex of side 'side'.View index ordered by [name] - [type]*/ return side.t.vertex(numarray[abs(side.n) + 1]); } point midpoint(side side) {/*View index ordered by [name] - [type]*/ return midpoint(segment(side)); } triangle operator *(transform T, triangle t) {/*Provide transform * triangle.View index ordered by [name] - [type]*/ return triangle(T * t.A, T * t.B, T * t.C); } triangle triangleAbc(real alpha, real b, real c, real angle = 0, point A = (0, 0)) {/*Return the triangle ABC rotated by 'angle' with BAC = alpha, AC = b and AB = c.View index ordered by [name] - [type]*/ triangle T; coordsys R = A.coordsys; T.init(A, A + R.polar(c, radians(angle)), A + R.polar(b, radians(angle + alpha))); return T; } triangle triangleabc(real a, real b, real c, real angle = 0, point A = (0, 0)) {/*Return the triangle ABC rotated by 'angle' with BC = a, AC = b and AB = c.View index ordered by [name] - [type]*/ triangle T; coordsys R = A.coordsys; T.init(A, A + R.polar(c, radians(angle)), A + R.polar(b, radians(angle) + acos((b^2 + c^2 - a^2)/(2 * b * c)))); return T; } triangle triangle(line l1, line l2, line l3) {/*Return the triangle defined by three line.View index ordered by [name] - [type]*/ point P1, P2, P3; P1 = intersectionpoint(l1, l2); P2 = intersectionpoint(l1, l3); P3 = intersectionpoint(l2, l3); if(!(defined(P1) && defined(P2) && defined(P3))) abort("triangle: two lines are parallel."); return triangle(P1, P2, P3); } point foot(vertex V) {/*Return the endpoint of the altitude from V.View index ordered by [name] - [type]*/ return projection((line)opposite(V)) * ((point)V); } point foot(side side) {/*Return the endpoint of the altitude on 'side'.View index ordered by [name] - [type]*/ return projection((line)side) * point(opposite(side)); } line altitude(vertex V) {/*Return the altitude passing through 'V'.View index ordered by [name] - [type]*/ return line(point(V), foot(V)); } line altitude(side side) {/*Return the altitude cutting 'side'.View index ordered by [name] - [type]*/ return altitude(opposite(side)); } point orthocentercenter(triangle t) {/*Return the orthocenter of the triangle t.View index ordered by [name] - [type]*/ return orthocentercenter(t.A, t.B, t.C); } point centroid(triangle t) {/*Return the centroid of the triangle 't'.View index ordered by [name] - [type]*/ return (t.A + t.B + t.C)/3; } point circumcenter(triangle t) {/*Return the circumcenter of the triangle 't'.View index ordered by [name] - [type]*/ return circumcenter(t.A, t.B, t.C); } circle circle(triangle t) {/*Return the circumcircle of the triangle 't'.View index ordered by [name] - [type]*/ return circle(t.A, t.B, t.C); } circle circumcircle(triangle t) {/*Return the circumcircle of the triangle 't'.View index ordered by [name] - [type]*/ return circle(t.A, t.B, t.C); } point incenter(triangle t) {/*Return the center of the incircle of the triangle 't'.View index ordered by [name] - [type]*/ return incenter(t.A, t.B, t.C); } real inradius(triangle t) {/*Return the radius of the incircle of the triangle 't'.View index ordered by [name] - [type]*/ return inradius(t.A, t.B, t.C); } circle incircle(triangle t) {/*Return the the incircle of the triangle 't'.View index ordered by [name] - [type]*/ return incircle(t.A, t.B, t.C); } point excenter(side side) {/*Return the center of the excircle tangent with the side 'side' of its triangle. side = 0 means AB, 1 means AC, other means BC. One must use the predefined sides t.AB, t.AC where 't' is a triangle....View index ordered by [name] - [type]*/ point op; triangle t = side.t; int n = numarray[abs(side.n) - 1]; if(n == 1) op = excenter(t.A, t.B, t.C); else if(n == 2) op = excenter(t.B, t.C, t.A); else op = excenter(t.C, t.A, t.B); return op; } real exradius(side side) {/*Return radius of the excircle tangent with the side 'side' of its triangle. side = 0 means AB, 1 means BC, other means CA. One must use the predefined sides t.AB, t.AC where 't' is a triangle....View index ordered by [name] - [type]*/ real or; triangle t = side.t; int n = numarray[abs(side.n) - 1]; if(n == 1) or = exradius(t.A, t.B, t.C); else if(n == 2) or = exradius(t.B, t.C, t.A); else or = exradius(t.A, t.C, t.B); return or; } circle excircle(side side) {/*Return the excircle tangent with the side 'side' of its triangle. side = 0 means AB, 1 means AC, other means BC. One must use the predefined sides t.AB, t.AC where 't' is a triangle....View index ordered by [name] - [type]*/ circle oc; int n = numarray[abs(side.n) - 1]; triangle t = side.t; if(n == 1) oc = excircle(t.A, t.B, t.C); else if(n == 2) oc = excircle(t.B, t.C, t.A); else oc = excircle(t.A, t.C, t.B); return oc; } struct trilinear {/*Trilinear coordinates 'a:b:c' relative to triangle 't'. http://mathworld.wolfram.com/TrilinearCoordinates.html*/ real a,b,c;/*The trilinear coordinates.*/ triangle t;/*The reference triangle.*/ }/*View index ordered by [name] - [type]*/ trilinear trilinear(triangle t, real a, real b, real c) {/*Return the trilinear coordinates relative to 't'. http://mathworld.wolfram.com/TrilinearCoordinates.htmlView index ordered by [name] - [type]*/ trilinear ot; ot.a = a; ot.b = b; ot.c = c; ot.t = t; return ot; } trilinear trilinear(triangle t, point M) {/*Return the trilinear coordinates of 'M' relative to 't'. http://mathworld.wolfram.com/TrilinearCoordinates.htmlView index ordered by [name] - [type]*/ trilinear ot; pair m = locate(M); int sameside(pair A, pair B, pair m, pair p) {// Return 1 if 'm' and 'p' are same side of line (AB) else return -1. pair mil = (A + B)/2; pair mA = rotate(90, mil) * A; pair mB = rotate(-90, mil) * A; return (abs(m - mA) <= abs(m - mB)) == (abs(p - mA) <= abs(p - mB)) ? 1 : -1; } real det(pair a, pair b) {return a.x * b.y - a.y * b.x;} real area(pair a, pair b, pair c){return 0.5 * abs(det(a, b) + det(b, c) + det(c, a));} pair A = t.A, B = t.B, C = t.C; real t1 = area(B, C, m), t2 = area(C, A, m), t3 = area(A, B, m); ot.a = sameside(B, C, A, m) * t1/t.a(); ot.b = sameside(A, C, B, m) * t2/t.b(); ot.c = sameside(A, B, C, m) * t3/t.c(); ot.t = t; return ot; } void write(trilinear tri) {/*View index ordered by [name] - [type]*/ write(format("%f : ", tri.a) + format("%f : ", tri.b) + format("%f", tri.c)); } point point(trilinear tri) {/*Return the trilinear coordinates relative to 't'. http://mathworld.wolfram.com/TrilinearCoordinates.htmlView index ordered by [name] - [type]*/ triangle t = tri.t; return masscenter(0.5 * t.a() * mass(t.A, tri.a), 0.5 * t.b() * mass(t.B, tri.b), 0.5 * t.c() * mass(t.C, tri.c)); } int[] tricoef(side side) {/*Return an array of integer (values are 0 or 1) which represents 'side'. For example, side = t.BC will be represented by {0, 1, 1}.View index ordered by [name] - [type]*/ int[] oi; int n = numarray[abs(side.n) - 1]; oi.push((n == 1 || n == 3) ? 1 : 0); oi.push((n == 1 || n == 2) ? 1 : 0); oi.push((n == 2 || n == 3) ? 1 : 0); return oi; } point operator cast(trilinear tri) {/*Cast trilinear to point. One may use the routine 'point(trilinear)' to force the casting.View index ordered by [name] - [type]*/ return point(tri); } typedef real centerfunction(real, real, real);/*http://mathworld.wolfram.com/TriangleCenterFunction.html*/ trilinear trilinear(triangle t, centerfunction f, real a = t.a(), real b = t.b(), real c = t.c()) {/*http://mathworld.wolfram.com/TriangleCenterFunction.htmlView index ordered by [name] - [type]*/ return trilinear(t, f(a, b, c), f(b, c, a), f(c, a, b)); } point symmedian(triangle t) {/*Return the symmedian point of 't'.View index ordered by [name] - [type]*/ point A, B, C; real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, 0, b, c); B = trilinear(t, a, 0, c); return intersectionpoint(line(t.A, A), line(t.B, B)); } point symmedian(side side) {/*The symmedian point on the side 'side'.View index ordered by [name] - [type]*/ triangle t = side.t; int n = numarray[abs(side.n) - 1]; if(n == 1) return trilinear(t, t.a(), t.b(), 0); if(n == 2) return trilinear(t, 0, t.b(), t.c()); return trilinear(t, t.a(), 0, t.c()); } line symmedian(vertex V) {/*Return the symmedian passing through 'V'.View index ordered by [name] - [type]*/ return line(point(V), symmedian(V.t)); } triangle cevian(triangle t, point P) {/*Return the Cevian triangle with respect of 'P' http://mathworld.wolfram.com/CevianTriangle.html.View index ordered by [name] - [type]*/ trilinear tri = trilinear(t, locate(P)); point A = point(trilinear(t, 0, tri.b, tri.c)); point B = point(trilinear(t, tri.a, 0, tri.c)); point C = point(trilinear(t, tri.a, tri.b, 0)); return triangle(A, B, C); } point cevian(side side, point P) {/*Return the Cevian point on 'side' with respect of 'P'.View index ordered by [name] - [type]*/ triangle t = side.t; trilinear tri = trilinear(t, locate(P)); int[] s = tricoef(side); return point(trilinear(t, s[0] * tri.a, s[1] * tri.b, s[2] * tri.c)); } line cevian(vertex V, point P) {/*Return line passing through 'V' and its Cevian image with respect of 'P'.View index ordered by [name] - [type]*/ return line(point(V), cevian(opposite(V), P)); } point gergonne(triangle t) {/*Return the Gergonne point of 't'.View index ordered by [name] - [type]*/ real f(real a, real b, real c){return 1/(a * (b + c - a));} return point(trilinear(t, f)); } point[] fermat(triangle t) {/*Return the Fermat points of 't'.View index ordered by [name] - [type]*/ point[] P; real A = t.alpha(), B = t.beta(), C = t.gamma(); P.push(point(trilinear(t, 1/Sin(A + 60), 1/Sin(B + 60), 1/Sin(C + 60)))); P.push(point(trilinear(t, 1/Sin(A - 60), 1/Sin(B - 60), 1/Sin(C - 60)))); return P; } point isotomicconjugate(triangle t, point M) {/*http://mathworld.wolfram.com/IsotomicConjugate.htmlView index ordered by [name] - [type]*/ if(!inside(t.Path(), locate(M))) abort("isotomic: the point must be inside the triangle."); trilinear tr = trilinear(t, M); return point(trilinear(t, 1/(t.a()^2 * tr.a), 1/(t.b()^2 * tr.b), 1/(t.c()^2 * tr.c))); } line isotomic(vertex V, point M) {/*http://mathworld.wolfram.com/IsotomicConjugate.html.View index ordered by [name] - [type]*/ side op = opposite(V); return line(V, rotate(180, midpoint(op)) * cevian(op, M)); } point isotomic(side side, point M) {/*http://mathworld.wolfram.com/IsotomicConjugate.htmlView index ordered by [name] - [type]*/ return intersectionpoint(isotomic(opposite(side), M), side); } triangle isotomic(triangle t, point M) {/*http://mathworld.wolfram.com/IsotomicConjugate.htmlView index ordered by [name] - [type]*/ return triangle(isotomic(t.BC, M), isotomic(t.CA, M), isotomic(t.AB, M)); } point isogonalconjugate(triangle t, point M) {/*http://mathworld.wolfram.com/IsogonalConjugate.htmlView index ordered by [name] - [type]*/ trilinear tr = trilinear(t, M); return point(trilinear(t, 1/tr.a, 1/tr.b, 1/tr.c)); } point isogonal(side side, point M) {/*http://mathworld.wolfram.com/IsogonalConjugate.htmlView index ordered by [name] - [type]*/ return cevian(side, isogonalconjugate(side.t, M)); } line isogonal(vertex V, point M) {/*http://mathworld.wolfram.com/IsogonalConjugate.htmlView index ordered by [name] - [type]*/ return line(V, isogonal(opposite(V), M)); } triangle isogonal(triangle t, point M) {/*http://mathworld.wolfram.com/IsogonalConjugate.htmlView index ordered by [name] - [type]*/ return triangle(isogonal(t.BC, M), isogonal(t.CA, M), isogonal(t.AB, M)); } triangle pedal(triangle t, point M) {/*Return the pedal triangle of 'M' in 't'. http://mathworld.wolfram.com/PedalTriangle.htmlView index ordered by [name] - [type]*/ return triangle(projection(t.BC) * M, projection(t.AC) * M, projection(t.AB) * M); } line pedal(side side, point M) {/*Return the pedal line of 'M' cutting 'side'. http://mathworld.wolfram.com/PedalTriangle.htmlView index ordered by [name] - [type]*/ return line(M, projection(side) * M); } triangle antipedal(triangle t, point M) {/*http://mathworld.wolfram.com/AntipedalTriangle.htmlView index ordered by [name] - [type]*/ trilinear Tm = trilinear(t, M); real a = Tm.a, b = Tm.b, c = Tm.c; real CA = Cos(t.alpha()), CB = Cos(t.beta()), CC = Cos(t.gamma()); point A = trilinear(t, -(b + a * CC) * (c + a * CB), (c + a * CB) * (a + b * CC), (b + a * CC) * (a + c * CB)); point B = trilinear(t, (c + b * CA) * (b + a * CC), -(c + b * CA) * (a + b * CC), (a + b * CC) * (b + c * CA)); point C = trilinear(t, (b + c * CA) * (c + a * CB), (a + c * CB) * (c + b * CA), -(a + c * CB) * (b + c * CA)); return triangle(A, B, C); } triangle extouch(triangle t) {/*Return the extouch triangle of the triangle 't'. The extouch triangle of 't' is the triangle formed by the points of tangency of a triangle 't' with its excircles.View index ordered by [name] - [type]*/ point A, B, C; real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, 0, (a - b + c)/b, (a + b - c)/c); B = trilinear(t, (-a + b + c)/a, 0, (a + b - c)/c); C = trilinear(t, (-a + b + c)/a, (a - b + c)/b, 0); return triangle(A, B, C); } triangle incentral(triangle t) {/*Return the incentral triangle of the triangle 't'. It is the triangle whose vertices are determined by the intersections of the reference triangle's angle bisectors with the respective opposite sides.View index ordered by [name] - [type]*/ point A, B, C; // real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, 0, 1, 1); B = trilinear(t, 1, 0, 1); C = trilinear(t, 1, 1, 0); return triangle(A, B, C); } triangle extouch(side side) {/*Return the triangle formed by the points of tangency of the triangle referenced by 'side' with its excircles. One vertex of the returned triangle is on the segment 'side'.View index ordered by [name] - [type]*/ triangle t = side.t; transform p1 = projection((line)t.AB); transform p2 = projection((line)t.AC); transform p3 = projection((line)t.BC); point EP = excenter(side); return triangle(p3 * EP, p2 * EP, p1 * EP); } point bisectorpoint(side side) {/*The intersection point of the angle bisector from the opposite point of 'side' with the side 'side'.View index ordered by [name] - [type]*/ triangle t = side.t; int n = numarray[abs(side.n) - 1]; if(n == 1) return trilinear(t, 1, 1, 0); if(n == 2) return trilinear(t, 0, 1, 1); return trilinear(t, 1, 0, 1); } line bisector(vertex V, real angle = 0) {/*Return the interior bisector passing through 'V' rotated by angle (in degrees) around 'V'.View index ordered by [name] - [type]*/ return rotate(angle, point(V)) * line(point(V), incenter(V.t)); } line bisector(side side) {/*Return the bisector of the line segment 'side'.View index ordered by [name] - [type]*/ return bisector(segment(side)); } point intouch(side side) {/*The point of tangency on the side 'side' of its incircle.View index ordered by [name] - [type]*/ triangle t = side.t; real a = t.a(), b = t.b(), c = t.c(); int n = numarray[abs(side.n) - 1]; if(n == 1) return trilinear(t, b * c/(-a + b + c), a * c/(a - b + c), 0); if(n == 2) return trilinear(t, 0, a * c/(a - b + c), a * b/(a + b - c)); return trilinear(t, b * c/(-a + b + c), 0, a * b/(a + b - c)); } triangle intouch(triangle t) {/*Return the intouch triangle of the triangle 't'. The intouch triangle of 't' is the triangle formed by the points of tangency of a triangle 't' with its incircles.View index ordered by [name] - [type]*/ point A, B, C; real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, 0, a * c/(a - b + c), a * b/(a + b - c)); B = trilinear(t, b * c/(-a + b + c), 0, a * b/(a + b - c)); C = trilinear(t, b * c/(-a + b + c), a * c/(a - b + c), 0); return triangle(A, B, C); } triangle tangential(triangle t) {/*Return the tangential triangle of the triangle 't'. The tangential triangle of 't' is the triangle formed by the lines tangent to the circumcircle of the given triangle 't' at its vertices.View index ordered by [name] - [type]*/ point A, B, C; real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, -a, b, c); B = trilinear(t, a, -b, c); C = trilinear(t, a, b, -c); return triangle(A, B, C); } triangle medial(triangle t) {/*Return the triangle whose vertices are midpoints of the sides of 't'.View index ordered by [name] - [type]*/ return triangle(midpoint(t.BC), midpoint(t.AC), midpoint(t.AB)); } line median(vertex V) {/*Return median from 'V'.View index ordered by [name] - [type]*/ return line(point(V), midpoint(segment(opposite(V)))); } line median(side side) {/*Return median from the opposite vertex of 'side'.View index ordered by [name] - [type]*/ return median(opposite(side)); } triangle orthic(triangle t) {/*Return the triangle whose vertices are endpoints of the altitudes from each of the vertices of 't'.View index ordered by [name] - [type]*/ return triangle(foot(t.BC), foot(t.AC), foot(t.AB)); } triangle symmedial(triangle t) {/*Return the symmedial triangle of 't'.View index ordered by [name] - [type]*/ point A, B, C; real a = t.a(), b = t.b(), c = t.c(); A = trilinear(t, 0, b, c); B = trilinear(t, a, 0, c); C = trilinear(t, a, b, 0); return triangle(A, B, C); } triangle anticomplementary(triangle t) {/*Return the triangle which has the given triangle 't' as its medial triangle.View index ordered by [name] - [type]*/ real a = t.a(), b = t.b(), c = t.c(); real ab = a * b, bc = b * c, ca = c * a; point A = trilinear(t, -bc, ca, ab); point B = trilinear(t, bc, -ca, ab); point C = trilinear(t, bc, ca, -ab); return triangle(A, B, C); } point[] intersectionpoints(triangle t, line l, bool extended = false) {/*Return the intersection points. If 'extended' is true, the sides are lines else the sides are segments. intersectionpoints(line, triangle, bool) is also defined.View index ordered by [name] - [type]*/ point[] OP; void addpoint(point P) { if(defined(P)) { bool exist = false; for (int i = 0; i < OP.length; ++i) { if(P == OP[i]) {exist = true; break;} } if(!exist) OP.push(P); } } if(extended) { for (int i = 1; i <= 3; ++i) { addpoint(intersectionpoint(t.line(i), l)); } } else { for (int i = 1; i <= 3; ++i) { addpoint(intersectionpoint((segment)t.line(i), l)); } } return OP; } point[] intersectionpoints(line l, triangle t, bool extended = false) { return intersectionpoints(t, l, extended); } vector dir(vertex V) {/*The direction (towards the outside of the triangle) of the interior angle bisector of 'V'.View index ordered by [name] - [type]*/ triangle t = V.t; if(V.n == 1) return vector(defaultcoordsys, (-dir(t.A--t.B, t.A--t.C))); if(V.n == 2) return vector(defaultcoordsys, (-dir(t.B--t.A, t.B--t.C))); return vector(defaultcoordsys, (-dir(t.C--t.A, t.C--t.B))); } void label(picture pic = currentpicture, Label L, vertex V, pair align = dir(V), real alignFactor = 1, pen p = nullpen, filltype filltype = NoFill) {/*Draw 'L' on picture 'pic' at vertex 'V' aligned by 'alignFactor * align'.View index ordered by [name] - [type]*/ label(pic, L, locate(point(V)), alignFactor * align, p, filltype); } void label(picture pic = currentpicture, Label LA = "A$", Label LB = "$B$", Label LC = "$C", triangle t, real alignAngle = 0, real alignFactor = 1, pen p = nullpen, filltype filltype = NoFill) {/*Draw labels LA, LB and LC aligned in the rotated (by 'alignAngle' in degrees) direction (towards the outside of the triangle) of the interior angle bisector of vertices. One can individually modify the alignment by setting the Label parameter 'align'.View index ordered by [name] - [type]*/ Label lla = LA.copy(); lla.align(lla.align, rotate(alignAngle) * locate(dir(t.VA))); label(pic, LA, t.VA, align = lla.align.dir, alignFactor = alignFactor, p, filltype); Label llb = LB.copy(); llb.align(llb.align, rotate(alignAngle) * locate(dir(t.VB))); label(pic, llb, t.VB, align = llb.align.dir, alignFactor = alignFactor, p, filltype); Label llc = LC.copy(); llc.align(llc.align, rotate(alignAngle) * locate(dir(t.VC))); label(pic, llc, t.VC, align = llc.align.dir, alignFactor = alignFactor, p, filltype); } void show(picture pic = currentpicture, Label LA = "A$", Label LB = "$B$", Label LC = "$C$", Label La = "$a$", Label Lb = "$b$", Label Lc = "$c\$",
triangle t, pen p = currentpen, filltype filltype = NoFill)
{/*Draw triangle and labels of sides and vertices.View index ordered by [name] - [type]*/
pair a = locate(t.A), b = locate(t.B), c = locate(t.C);
draw(pic, a--b--c--cycle, p);
label(pic, LA, a, -dir(a--b, a--c), p, filltype);
label(pic, LB, b, -dir(b--a, b--c), p, filltype);
label(pic, LC, c, -dir(c--a, c--b), p, filltype);
pair aligna = I * unit(c - b), alignb = I * unit(c - a), alignc = I * unit(b - a);
pair mAB = locate(midpoint(t.AB)), mAC = locate(midpoint(t.AC)), mBC = locate(midpoint(t.BC));
draw(pic, La, b--c, align = rotate(dot(a - mBC, aligna) > 0 ? 180 :0) * aligna, p);
draw(pic, Lb, a--c, align = rotate(dot(b - mAC, alignb) > 0 ? 180 :0) * alignb, p);
draw(pic, Lc, a--b, align = rotate(dot(c - mAB, alignc) > 0 ? 180 :0) * alignc, p);
}

void draw(picture pic = currentpicture, triangle t, pen p = currentpen, marker marker = nomarker)
{/*Draw sides of the triangle 't' on picture 'pic' using pen 'p'.View index ordered by [name] - [type]*/
draw(pic, t.Path(), p, marker);
}

void draw(picture pic = currentpicture, triangle[] t, pen p = currentpen, marker marker = nomarker)
{/*Draw sides of the triangles 't' on picture 'pic' using pen 'p'.View index ordered by [name] - [type]*/
for(int i = 0; i < t.length; ++i) draw(pic, t[i], p, marker);
}

void drawline(picture pic = currentpicture, triangle t, pen p = currentpen)
{/*Draw lines of the triangle 't' on picture 'pic' using pen 'p'.View index ordered by [name] - [type]*/
draw(t, p);
draw(pic, line(t.A, t.B), p);
draw(pic, line(t.A, t.C), p);
draw(pic, line(t.B, t.C), p);
}

void dot(picture pic = currentpicture, triangle t, pen p = currentpen)
{/*Draw a dot at each vertex of 't'.View index ordered by [name] - [type]*/
dot(pic, t.A^^t.B^^t.C, p);
}
// *.......................TRIANGLES.......................*
// *=======================================================*

// *=======================================================*
// *.......................INVERSIONS......................*

point inverse(real k, point A, point M)
{/*Return the inverse point of 'M' with respect to point A and inversion radius 'k'.View index ordered by [name] - [type]*/
return A + k/conj(M - A);
}

{/*http://fr.wikipedia.org/wiki/Puissance_d'un_point_par_rapport_%C3%A0_un_cercleView index ordered by [name] - [type]*/
point[] P = standardizecoordsys(c1.C, c2.C);
real k = c1.r^2 - c2.r^2;
pair C1 = locate(c1.C);
pair C2 = locate(c2.C);
pair oop = C2 - C1;
pair K = (abs(oop) == 0) ?
(infinity, infinity) :
midpoint(C1--C2) + 0.5 * k * oop/dot(oop, oop);
return point(P[0].coordsys, K/P[0].coordsys);
}

{/*http://fr.wikipedia.org/wiki/Puissance_d'un_point_par_rapport_%C3%A0_un_cercleView index ordered by [name] - [type]*/
if (c1.C == c2.C) abort("radicalline: the centers must be distinct");
}

point radicalcenter(circle c1, circle c2, circle c3)
{/*http://fr.wikipedia.org/wiki/Puissance_d'un_point_par_rapport_%C3%A0_un_cercleView index ordered by [name] - [type]*/
}

struct inversion
{/*http://mathworld.wolfram.com/Inversion.html*/
point C;
real k;
}/*View index ordered by [name] - [type]*/

inversion inversion(real k, point C)
{/*Return the inversion with respect to 'C' having inversion radius 'k'.View index ordered by [name] - [type]*/
inversion oi;
oi.k = k;
oi.C = C;
return oi;
}

inversion inversion(point C, real k)
{/*Return the inversion with respect to 'C' having inversion radius 'k'.View index ordered by [name] - [type]*/
return inversion(k, C);
}

inversion inversion(circle c1, circle c2, real sgn = 1)
{/*Return the inversion which transforms 'c1' to
. 'c2' and positive inversion radius if 'sgn > 0';
. 'c2' and negative inversion radius if 'sgn < 0';
. 'c1' and 'c2' to 'c2' if 'sgn = 0'.View index ordered by [name] - [type]*/
if(sgn == 0) {
return inversion(O^c1, O);
}
real a = abs(c1.r/c2.r);
if(sgn > 0) {
point O = c1.C + a/abs(1 - a) * (c2.C - c1.C);
return inversion(a * abs(abs(O - c2.C)^2 - c2.r^2), O);
}
point O = c1.C + a/abs(1 + a) * (c2.C - c1.C);
return inversion(-a * abs(abs(O - c2.C)^2 - c2.r^2), O);
}

inversion inversion(circle c1, circle c2, circle c3)
{/*Return the inversion which transform 'c1' to 'c1', 'c2' to 'c2' and 'c3' to 'c3'.View index ordered by [name] - [type]*/
point Rc = radicalcenter(c1, c2, c3);
return inversion(Rc, Rc^c1);
}

circle operator cast(inversion i){return circle(i.C, sgn(i.k) * sqrt(abs(i.k)));}

circle circle(inversion i)
{/*Return the inversion circle of 'i'.View index ordered by [name] - [type]*/
return i;
}

inversion operator cast(circle c)
{
return inversion(sgn(c.r) * c.r^2, c.C);
}

inversion inversion(circle c)
{/*Return the inversion represented by the circle of 'c'.View index ordered by [name] - [type]*/
return c;
}

point operator *(inversion i, point P)
{/*Provide inversion * point.View index ordered by [name] - [type]*/
return inverse(i.k, i.C, P);
}

void lineinversion()
{
warning("lineinversion", "the inversion of the line is not a circle.
The returned circle has an infinite radius, circle.l has been set.");
}

circle inverse(real k, point A, line l)
{/*Return the inverse circle of 'l' with
respect to point 'A' and inversion radius 'k'.View index ordered by [name] - [type]*/
if(A @ l) {
lineinversion();
circle C = circle(A, infinity);
C.l = l;
return C;
}
point Ap = inverse(k, A, l.A), Bp = inverse(k, A, l.B);
return circle(A, Ap, Bp);
}

circle operator *(inversion i, line l)
{/*Provide inversion * line for lines that don't pass through the inversion center.View index ordered by [name] - [type]*/
return inverse(i.k, i.C, l);
}

circle inverse(real k, point A, circle c)
{/*Return the inverse circle of 'c' with
respect to point A and inversion radius 'k'.View index ordered by [name] - [type]*/
if(degenerate(c)) return inverse(k, A, c.l);
if(A @ c) {
lineinversion();
point M = rotate(180, c.C) * A, Mp = rotate(90, c.C) * A;
circle oc = circle(A, infinity);
oc.l = line(inverse(k, A, M), inverse(k, A, Mp));
return oc;
}
point[] P = standardizecoordsys(A, c.C);
real s = k/((P[1].x - P[0].x)^2 + (P[1].y - P[0].y)^2 - c.r^2);
return circle(P[0] + s * (P[1]-P[0]), abs(s) * c.r);
}

circle operator *(inversion i, circle c)
{/*Provide inversion * circle.View index ordered by [name] - [type]*/
return inverse(i.k, i.C, c);
}
// *.......................INVERSIONS......................*
// *=======================================================*

// *=======================================================*
// *........................FOOTER.........................*

point[] intersectionpoints(line l, circle c)
{/*Note that the line 'l' may be a segment by casting.
intersectionpoints(circle, line) is also defined.View index ordered by [name] - [type]*/
if(degenerate(c)) return new point[]{intersectionpoint(l, c.l)};
point[] op;
coordsys R = samecoordsys(l.A, c.C) ?
l.A.coordsys : defaultcoordsys;
coordsys Rp = defaultcoordsys;
circle cc = circle(changecoordsys(Rp, c.C), c.r);
point proj = projection(l) * c.C;
if(proj @ cc) { // The line is a tangente of the circle.
if(proj @ l) op.push(proj);// line may be a segement...
} else {
coordsys Rc = cartesiansystem(c.C, (1, 0), (0, 1));
line ll = changecoordsys(Rc, l);
pair[] P = intersectionpoints(ll.A.coordinates, ll.B.coordinates,
1, 0, 1, 0, 0, -c.r^2);
for (int i = 0; i < P.length; ++i) {
point inter = changecoordsys(R, point(Rc, P[i]));
if(inter @ l) op.push(inter);
}
}
return op;
}

point[] intersectionpoints(circle c, line l)
{
return intersectionpoints(l, c);
}

point[] intersectionpoints(line l, ellipse el)
{/*Note that the line 'l' may be a segment by casting.
intersectionpoints(ellipse, line) is also defined.View index ordered by [name] - [type]*/
if(el.e == 0) return intersectionpoints(l, (circle)el);
if(degenerate(el)) return new point[]{intersectionpoint(l, el.l)};
point[] op;
coordsys R = samecoordsys(l.A, el.C) ? l.A.coordsys : defaultcoordsys;
coordsys Rp = defaultcoordsys;
line ll = changecoordsys(Rp, l);
ellipse ell = changecoordsys(Rp, el);
circle C = circle(ell.C, ell.a);
point[] Ip = intersectionpoints(ll, C);
if (Ip.length > 0 &&
(perpendicular(ll, line(ell.F1, Ip[0])) ||
perpendicular(ll, line(ell.F2, Ip[0])))) {
// http://www.mathcurve.com/courbes2d/ellipse/ellipse.shtml
//  Définition tangentielle par antipodaire de cercle.
// 'l' is a tangent of 'el'
transform t = scale(el.a/el.b, el.F1, el.F2, el.C, rotate(90, el.C) * el.F1);
point inter = inverse(t) * intersectionpoints(C, t * ll)[0];
if(inter @ l) op.push(inter);
} else {
coordsys Rc = canonicalcartesiansystem(el);
line ll = changecoordsys(Rc, l);
pair[] P = intersectionpoints(ll.A.coordinates, ll.B.coordinates,
1/el.a^2, 0, 1/el.b^2, 0, 0, -1);
for (int i = 0; i < P.length; ++i) {
point inter = changecoordsys(R, point(Rc, P[i]));
if(inter @ l) op.push(inter);
}
}
return op;
}

point[] intersectionpoints(ellipse el, line l)
{
return intersectionpoints(l, el);
}

point[] intersectionpoints(line l, parabola p)
{/*Note that the line 'l' may be a segment by casting.
intersectionpoints(parabola, line) is also defined.View index ordered by [name] - [type]*/
point[] op;
coordsys R = coordsys(p);
bool tgt = false;
line ll = changecoordsys(R, l),
lv = parallel(p.V, p.D);
point M = intersectionpoint(lv, ll), tgtp;
if(finite(M)) {// Test if 'l' is tangent to 'p'
line l1 = bisector(line(M, p.F));
line l2 = rotate(90, M) * lv;
point P = intersectionpoint(l1, l2);
tgtp = rotate(180, P) * p.F;
tgt = (tgtp @ l);
}
if(tgt) {
if(tgtp @ l) op.push(tgtp);
} else {
real[] eq = changecoordsys(defaultcoordsys, equation(p)).a;
pair[] tp = intersectionpoints(locate(l.A), locate(l.B), eq);
point inter;
for (int i = 0; i < tp.length; ++i) {
inter = point(R, tp[i]/R);
if(inter @ l) op.push(inter);
}
}
return op;
}

point[] intersectionpoints(parabola p, line l)
{
return intersectionpoints(l, p);
}

point[] intersectionpoints(line l, hyperbola h)
{/*Note that the line 'l' may be a segment by casting.
intersectionpoints(hyperbola, line) is also defined.View index ordered by [name] - [type]*/
point[] op;
coordsys R = coordsys(h);
point A = intersectionpoint(l, h.A1), B = intersectionpoint(l, h.A2);
point M = midpoint(segment(A, B));
bool tgt = M @ h;
if(tgt) {
if(M @ l) op.push(M);
} else {
real[] eq = changecoordsys(defaultcoordsys, equation(h)).a;
pair[] tp = intersectionpoints(locate(l.A), locate(l.B), eq);
point inter;
for (int i = 0; i < tp.length; ++i) {
inter = point(R, tp[i]/R);
if(inter @ l) op.push(inter);
}
}
return op;
}

point[] intersectionpoints(hyperbola h, line l)
{
return intersectionpoints(l, h);
}

point[] intersectionpoints(line l, conic co)
{/*Note that the line 'l' may be a segment by casting.
intersectionpoints(conic, line) is also defined.View index ordered by [name] - [type]*/
point[] op;
if(co.e < 1) op = intersectionpoints((ellipse)co, l);
else
if(co.e == 1) op = intersectionpoints((parabola)co, l);
else op = intersectionpoints((hyperbola)co, l);
return op;
}

point[] intersectionpoints(conic co, line l)
{
return intersectionpoints(l, co);
}

point[] intersectionpoints(conic co1, conic co2)
{/*Return the intersection points of the two conics.View index ordered by [name] - [type]*/
if(degenerate(co1)) return intersectionpoints(co1.l[0], co2);
if(degenerate(co2)) return intersectionpoints(co1, co2.l[0]);
return intersectionpoints(equation(co1), equation(co2));
}

point[] intersectionpoints(triangle t, conic co, bool extended = false)
{/*Return the intersection points.
If 'extended' is true, the sides are lines else the sides are segments.
intersectionpoints(conic, triangle, bool) is also defined.View index ordered by [name] - [type]*/
if(degenerate(co)) return intersectionpoints(t, co.l[0], extended);
point[] OP;
{
for (int i = 0; i < P.length; ++i) {
if(defined(P[i])) {
bool exist = false;
for (int j = 0; j < OP.length; ++j) {
if(P[i] == OP[j]) {exist = true; break;}
}
if(!exist) OP.push(P[i]);
}}}
if(extended) {
for (int i = 1; i <= 3; ++i) {
}
} else {
for (int i = 1; i <= 3; ++i) {
}
}
return OP;
}

point[] intersectionpoints(conic co, triangle t, bool extended = false)
{
return intersectionpoints(t, co, extended);
}

point[] intersectionpoints(ellipse a, ellipse b)
{/*View index ordered by [name] - [type]*/
// if(degenerate(a)) return intersectionpoints(a.l, b);
// if(degenerate(b)) return intersectionpoints(a, b.l);;
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(ellipse a, circle b)
{/*View index ordered by [name] - [type]*/
// if(degenerate(a)) return intersectionpoints(a.l, b);
// if(degenerate(b)) return intersectionpoints(a, b.l);;
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(circle a, ellipse b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints(b, a);
}

point[] intersectionpoints(ellipse a, parabola b)
{/*View index ordered by [name] - [type]*/
// if(degenerate(a)) return intersectionpoints(a.l, b);
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(parabola a, ellipse b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints(b, a);
}

point[] intersectionpoints(ellipse a, hyperbola b)
{/*View index ordered by [name] - [type]*/
// if(degenerate(a)) return intersectionpoints(a.l, b);
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(hyperbola a, ellipse b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints(b, a);
}

point[] intersectionpoints(circle a, parabola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(parabola a, circle b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(circle a, hyperbola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(hyperbola a, circle b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(parabola a, parabola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(parabola a, hyperbola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(hyperbola a, parabola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(hyperbola a, hyperbola b)
{/*View index ordered by [name] - [type]*/
return intersectionpoints((conic)a, (conic)b);
}

point[] intersectionpoints(circle c1, circle c2)
{/*View index ordered by [name] - [type]*/
if(degenerate(c1))
return degenerate(c2) ?
new point[]{intersectionpoint(c1.l, c2.l)} : intersectionpoints(c1.l, c2);
if(degenerate(c2)) return intersectionpoints(c1, c2.l);
return (c1.C == c2.C) ?
new point[] :
}

line tangent(circle c, abscissa x)
{/*Return the tangent of 'c' at 'point(c, x)'.View index ordered by [name] - [type]*/
if(c.r == 0) abort("tangent: a circle with a radius equals zero has no tangent.");
point M = point(c, x);
return line(rotate(90, M) * c.C, M);
}

line[] tangents(circle c, point M)
{/*Return the tangents of 'c' passing through 'M'.View index ordered by [name] - [type]*/
line[] ol;
if(inside(c, M)) return ol;
if(M @ c) {
ol.push(tangent(c, relabscissa(c, M)));
} else {
circle cc = circle(c.C, M);
point[] inter = intersectionpoints(c, cc);
for (int i = 0; i < inter.length; ++i)
ol.push(tangents(c, inter[i])[0]);
}
return ol;
}

point point(circle c, point M)
{/*Return the intersection point of 'c'
with the half-line '[c.C M)'.View index ordered by [name] - [type]*/
return intersectionpoints(c, line(c.C, false, M))[0];
}

line tangent(circle c, point M)
{/*Return the tangent of 'c' at the
intersection point of the half-line'[c.C M)'.View index ordered by [name] - [type]*/
return tangents(c, point(c, M))[0];
}

point point(circle c, explicit vector v)
{/*Return the intersection point of 'c'
with the half-line '[c.C v)'.View index ordered by [name] - [type]*/
return point(c, c.C + v);
}

line tangent(circle c, explicit vector v)
{/*Return the tangent of 'c' at the
point M so that vec(c.C M) is collinear to 'v' with the same sense.View index ordered by [name] - [type]*/
line ol = tangent(c, c.C + v);
return dot(ol.v, v) > 0 ? ol : reverse(ol);
}

line tangent(ellipse el, abscissa x)
{/*Return the tangent of 'el' at 'point(el, x)'.View index ordered by [name] - [type]*/
point M = point(el, x);
line l1 = line(el.F1, M);
line l2 = line(el.F2, M);
line ol = (l1 == l2) ? perpendicular(M, l1) : bisector(l1, l2, 90, false);
return ol;
}

line[] tangents(ellipse el, point M)
{/*Return the tangents of 'el' passing through 'M'.View index ordered by [name] - [type]*/
line[] ol;
if(inside(el, M)) return ol;
if(M @ el) {
ol.push(tangent(el, relabscissa(el, M)));
} else {
point Mp = samecoordsys(M, el.F2) ?
M : changecoordsys(el.F2.coordsys, M);
circle c = circle(Mp, abs(el.F1 - Mp));
circle cc = circle(el.F2, 2 * el.a);
point[] inter = intersectionpoints(c, cc);
for (int i = 0; i < inter.length; ++i) {
line tl = line(inter[i], el.F2, false);
point[] P = intersectionpoints(tl, el);
ol.push(line(Mp, P[0]));
}
}
return ol;
}

line tangent(parabola p, abscissa x)
{/*Return the tangent of 'p' at 'point(p, x)' (use the Wells method).View index ordered by [name] - [type]*/
line lt = rotate(90, p.V) * line(p.V, p.F);
point P = point(p, x);
if(P == p.V) return lt;
point M = midpoint(segment(P, p.F));
line l = rotate(90, M) * line(P, p.F);
return line(P, projection(lt) * M);
}

line[] tangents(parabola p, point M)
{/*Return the tangent of 'p' at 'M' (use the Wells method).View index ordered by [name] - [type]*/
line[] ol;
if(inside(p, M)) return ol;
if(M @ p) {
ol.push(tangent(p, angabscissa(p, M)));
}
else {
point Mt = changecoordsys(coordsys(p), M);
circle c = circle(Mt, p.F);
line l = rotate(90, p.V) * line(p.V, p.F);
point[] R = intersectionpoints(l, c);
for (int i = 0; i < R.length; ++i) {
ol.push(line(Mt, R[i]));
}
// An other method: http://www.du.edu/~jcalvert/math/parabola.htm
//   point[] R = intersectionpoints(p.directrix, c);
//   for (int i = 0; i < R.length; ++i) {
//     ol.push(bisector(segment(p.F, R[i])));
//   }
}
return ol;
}

line tangent(hyperbola h, abscissa x)
{/*Return the tangent of 'h' at 'point(p, x)'.View index ordered by [name] - [type]*/
point M = point(h, x);
line ol = bisector(line(M, h.F1), line(M, h.F2));
if(sameside(h.F1, h.F2, ol) || ol == line(h.F1, h.F2)) ol = rotate(90, M) * ol;
return ol;
}

line[] tangents(hyperbola h, point M)
{/*Return the tangent of 'h' at 'M'.View index ordered by [name] - [type]*/
line[] ol;
if(M @ h) {
ol.push(tangent(h, angabscissa(h, M, fromCenter)));
} else {
coordsys cano = canonicalcartesiansystem(h);
bqe bqe = changecoordsys(cano, equation(h));
real a = abs(1/(bqe.a[5] * bqe.a[0])), b = abs(1/(bqe.a[5] * bqe.a[2]));
point Mp = changecoordsys(cano, M);
real x0 = Mp.x, y0 = Mp.y;
if(abs(x0) > epsgeo) {
real c0 = a * y0^2/(b * x0)^2 - 1/b,
c1 = 2 * a * y0/(b * x0^2), c2 = a/x0^2 - 1;
real[] sol = quadraticroots(c0, c1, c2);
for (real y:sol) {
point tmp = changecoordsys(coordsys(h), point(cano, (a * (1 + y * y0/b)/x0, y)));
ol.push(line(M, tmp));
}
} else if(abs(y0) > epsgeo) {
real y = -b/y0, x = sqrt(a * (1 + b/y0^2));
ol.push(line(M, changecoordsys(coordsys(h), point(cano, (x, y)))));
ol.push(line(M, changecoordsys(coordsys(h), point(cano, (-x, y)))));
}}
return ol;
}

point[] intersectionpoints(conic co, arc a)
{/*intersectionpoints(arc, circle) is also defined.View index ordered by [name] - [type]*/
point[] op;
point[] tp = intersectionpoints(co, (conic)a.el);
for (int i = 0; i < tp.length; ++i)
if(tp[i] @ a) op.push(tp[i]);
return op;
}

point[] intersectionpoints(arc a, conic co)
{
return intersectionpoints(co, a);
}

point[] intersectionpoints(arc a1, arc a2)
{/*View index ordered by [name] - [type]*/
point[] op;
point[] tp = intersectionpoints(a1.el, a2.el);
for (int i = 0; i < tp.length; ++i)
if(tp[i] @ a1 && tp[i] @ a2) op.push(tp[i]);
return op;
}

point[] intersectionpoints(line l, arc a)
{/*intersectionpoints(arc, line) is also defined.View index ordered by [name] - [type]*/
point[] op;
point[] tp = intersectionpoints(a.el, l);
for (int i = 0; i < tp.length; ++i)
if(tp[i] @ a && tp[i] @ l) op.push(tp[i]);
return op;
}

point[] intersectionpoints(arc a, line l)
{
return intersectionpoints(l, a);
}

point arcsubtendedcenter(point A, point B, real angle)
{/*Return the center of the arc retuned
by the 'arcsubtended' routine.View index ordered by [name] - [type]*/
point OM;
point[] P = standardizecoordsys(A, B);
angle = angle%(sgnd(angle) * 180);
line bis = bisector(P[0], P[1]);
line AB = line(P[0], P[1]);
return intersectionpoint(bis, rotate(90 - angle, A) * AB);
}

arc arcsubtended(point A, point B, real angle)
{/*Return the arc circle from which the segment AB is saw with
the angle 'angle'.
If the point 'M' is on this arc, the oriented angle (MA, MB) is
equal to 'angle'.View index ordered by [name] - [type]*/
point[] P = standardizecoordsys(A, B);
line AB = line(P[0], P[1]);
angle = angle%(sgnd(angle) * 180);
point C = arcsubtendedcenter(P[0], P[1], angle);
real BC = degrees(B - C)%360;
real AC = degrees(A - C)%360;
return arc(circle(C, abs(B - C)), BC, AC, angle > 0 ? CCW : CW);
}

arc arccircle(point A, point M, point B)
{/*Return the CCW arc circle 'AB' passing through 'M'.View index ordered by [name] - [type]*/
circle tc = circle(A, M, B);
real a = degrees(A - tc.C);
real b = degrees(B - tc.C);
real m = degrees(M - tc.C);

arc oa = arc(tc, a, b);
// TODO : use cross product to determine CWW or CW
if (!(M @ oa)) {
oa.direction = !oa.direction;
}

return oa;
}

arc arc(ellipse el, explicit abscissa x1, explicit abscissa x2, bool direction = CCW)
{/*Return the arc from 'point(c, x1)' to 'point(c, x2)' in the direction 'direction'.View index ordered by [name] - [type]*/
real a = degrees(point(el, x1) - el.C);
real b = degrees(point(el, x2) - el.C);
arc oa = arc(el, a - el.angle, b - el.angle, fromCenter, direction);
return oa;
}

arc arc(ellipse el, point M, point N, bool direction = CCW)
{/*Return the arc from 'M' to 'N' in the direction 'direction'.
The points 'M' and 'N' must belong to the ellipse 'el'.View index ordered by [name] - [type]*/
return arc(el, relabscissa(el, M), relabscissa(el, N), direction);
}

arc arccircle(point A, point B, real angle, bool direction = CCW)
{/*Return the arc circle centered on A
from B to rotate(angle, A) * B in the direction 'direction'.View index ordered by [name] - [type]*/
point M = rotate(angle, A) * B;
return arc(circle(A, abs(A - B)), B, M, direction);
}

arc arc(explicit arc a, abscissa x1, abscissa x2)
{/*Return the arc from 'point(a, x1)' to 'point(a, x2)' traversed in the direction of the arc direction.View index ordered by [name] - [type]*/
real a1 = angabscissa(a.el, point(a, x1), a.polarconicroutine).x;
real a2 = angabscissa(a.el, point(a, x2), a.polarconicroutine).x;
return arc(a.el, a1, a2, a.polarconicroutine, a.direction);
}

arc arc(explicit arc a, point M, point N)
{/*Return the arc from 'M' to 'N'.
The points 'M' and 'N' must belong to the arc 'a'.View index ordered by [name] - [type]*/
return arc(a, relabscissa(a, M), relabscissa(a, N));
}

arc inverse(real k, point A, segment s)
{/*Return the inverse arc circle of 's'
with respect to point A and inversion radius 'k'.View index ordered by [name] - [type]*/
point Ap = inverse(k, A, s.A), Bp = inverse(k, A, s.B),
M = inverse(k, A, midpoint(s));
return arccircle(Ap, M, Bp);
}

arc operator *(inversion i, segment s)
{/*Provide
inversion * segment.View index ordered by [name] - [type]*/
return inverse(i.k, i.C, s);
}

path operator *(inversion i, triangle t)
{/*Provide inversion * triangle.View index ordered by [name] - [type]*/
return (path)(i * segment(t.AB))--
(path)(i * segment(t.BC))--
(path)(i * segment(t.CA))&cycle;
}

path compassmark(pair O, pair A, real position, real angle = 10)
{/*Return an arc centered on O with the angle 'angle' so that the position
of 'A' on this arc makes an angle 'position * angle'.View index ordered by [name] - [type]*/
real a = degrees(A - O);
real pa = (a - position * angle)%360,
pb = (a - (position - 1) * angle)%360;
real t1 = intersect(unitcircle, (0, 0)--2 * dir(pa))[0];
real t2 = intersect(unitcircle, (0, 0)--2 * dir(pb))[0];
int n = length(unitcircle);
if(t1 >= t2) t1 -= n;
return shift(O) * scale(abs(O - A)) * subpath(unitcircle, t1, t2);
}

line tangent(explicit arc a, abscissa x)
{/*Return the tangent of 'a' at 'point(a, x)'.View index ordered by [name] - [type]*/
abscissa ag = angabscissa(a, point(a, x));
return tangent(a.el, ag + a.angle1 + (a.el.e == 0 ? a.angle0 : 0));
}

line tangent(explicit arc a, point M)
{/*Return the tangent of 'a' at 'M'.
The points 'M' must belong to the arc 'a'.View index ordered by [name] - [type]*/
return tangent(a, angabscissa(a, M));
}

// *=======================================================*
// *.......Routines for compatibility with original geometry module........*

path square(pair z1, pair z2)
{
pair v = z2 - z1;
pair z3 = z2 + I * v;
pair z4 = z3 - v;
return z1--z2--z3--z4--cycle;
}

// Draw a perpendicular symbol at z aligned in the direction align
// relative to the path z--z + dir.
void perpendicular(picture pic = currentpicture, pair z, pair align,
pair dir = E, real size = 0, pen p = currentpen,
margin margin = NoMargin, filltype filltype = NoFill)
{
perpendicularmark(pic, (point) z, align, dir, size, p, margin, filltype);
}

// Draw a perpendicular symbol at z aligned in the direction align
// relative to the path z--z + dir(g, 0)
void perpendicular(picture pic = currentpicture, pair z, pair align, path g,
real size = 0, pen p = currentpen, margin margin = NoMargin,
filltype filltype = NoFill)
{
perpendicularmark(pic, (point) z, align, dir(g, 0), size, p, margin, filltype);
}

// Return an interior arc BAC of triangle ABC, given a radius r > 0.
// If r < 0, return the corresponding exterior arc of radius |r|.
path arc(explicit pair B, explicit pair A, explicit pair C, real r)
{
return arc(A, r, degrees(B - A), degrees(C - A));
}

// *.......End of compatibility routines........*
// *=======================================================*

// *........................FOOTER.........................*
// *=======================================================*