1 /*
  2     Copyright 2008-2025
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Bianca Valentin,
  7         Andreas Walter,
  8         Alfred Wassermann,
  9         Peter Wilfahrt
 10 
 11     This file is part of JSXGraph.
 12 
 13     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
 14 
 15     You can redistribute it and/or modify it under the terms of the
 16 
 17       * GNU Lesser General Public License as published by
 18         the Free Software Foundation, either version 3 of the License, or
 19         (at your option) any later version
 20       OR
 21       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
 22 
 23     JSXGraph is distributed in the hope that it will be useful,
 24     but WITHOUT ANY WARRANTY; without even the implied warranty of
 25     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 26     GNU Lesser General Public License for more details.
 27 
 28     You should have received a copy of the GNU Lesser General Public License and
 29     the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/>
 30     and <https://opensource.org/licenses/MIT/>.
 31  */
 32 
 33 /*global JXG: true, define: true*/
 34 /*jslint nomen: true, plusplus: true*/
 35 
 36 /**
 37  * @fileoverview This file contains the Math.Geometry namespace for calculating algebraic/geometric
 38  * stuff like intersection points, angles, midpoint, and so on.
 39  */
 40 
 41 import JXG from "../jxg.js";
 42 import Const from "../base/constants.js";
 43 import Coords from "../base/coords.js";
 44 import Mat from "./math.js";
 45 import Stat from "../math/statistics.js";
 46 import Numerics from "./numerics.js";
 47 import Type from "../utils/type.js";
 48 import Expect from "../utils/expect.js";
 49 
 50 /**
 51  * Math.Geometry namespace definition. This namespace holds geometrical algorithms,
 52  * especially intersection algorithms.
 53  * @name JXG.Math.Geometry
 54  * @exports Mat.Geometry as JXG.Math.Geometry
 55  * @namespace
 56  */
 57 Mat.Geometry = {};
 58 
 59 // the splitting is necessary due to the shortcut for the circumcircleMidpoint method to circumcenter.
 60 
 61 JXG.extend(
 62     Mat.Geometry,
 63     /** @lends JXG.Math.Geometry */ {
 64         /* ***************************************/
 65         /* *** GENERAL GEOMETRIC CALCULATIONS ****/
 66         /* ***************************************/
 67 
 68         /**
 69          * Calculates the angle defined by the points A, B, C.
 70          * @param {JXG.Point|Array} A A point  or [x,y] array.
 71          * @param {JXG.Point|Array} B Another point or [x,y] array.
 72          * @param {JXG.Point|Array} C A circle - no, of course the third point or [x,y] array.
 73          * @deprecated Use {@link JXG.Math.Geometry.rad} instead.
 74          * @see JXG.Math.Geometry.rad
 75          * @see JXG.Math.Geometry.trueAngle
 76          * @returns {Number} The angle in radian measure.
 77          */
 78         angle: function (A, B, C) {
 79             var u,
 80                 v,
 81                 s,
 82                 t,
 83                 a = [],
 84                 b = [],
 85                 c = [];
 86 
 87             JXG.deprecated("Geometry.angle()", "Geometry.rad()");
 88             if (A.coords) {
 89                 a[0] = A.coords.usrCoords[1];
 90                 a[1] = A.coords.usrCoords[2];
 91             } else {
 92                 a[0] = A[0];
 93                 a[1] = A[1];
 94             }
 95 
 96             if (B.coords) {
 97                 b[0] = B.coords.usrCoords[1];
 98                 b[1] = B.coords.usrCoords[2];
 99             } else {
100                 b[0] = B[0];
101                 b[1] = B[1];
102             }
103 
104             if (C.coords) {
105                 c[0] = C.coords.usrCoords[1];
106                 c[1] = C.coords.usrCoords[2];
107             } else {
108                 c[0] = C[0];
109                 c[1] = C[1];
110             }
111 
112             u = a[0] - b[0];
113             v = a[1] - b[1];
114             s = c[0] - b[0];
115             t = c[1] - b[1];
116 
117             return Math.atan2(u * t - v * s, u * s + v * t);
118         },
119 
120         /**
121          * Calculates the angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
122          * @param {JXG.Point|Array} A Point or [x,y] array
123          * @param {JXG.Point|Array} B Point or [x,y] array
124          * @param {JXG.Point|Array} C Point or [x,y] array
125          * @see JXG.Math.Geometry.rad
126          * @returns {Number} The angle in degrees.
127          */
128         trueAngle: function (A, B, C) {
129             return this.rad(A, B, C) * 57.295779513082323; // *180.0/Math.PI;
130         },
131 
132         /**
133          * Calculates the internal angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
134          * @param {JXG.Point|Array} A Point or [x,y] array
135          * @param {JXG.Point|Array} B Point or [x,y] array
136          * @param {JXG.Point|Array} C Point or [x,y] array
137          * @see JXG.Math.Geometry.trueAngle
138          * @returns {Number} Angle in radians.
139          */
140         rad: function (A, B, C) {
141             var ax, ay, bx, by, cx, cy, phi;
142 
143             if (A.coords) {
144                 ax = A.coords.usrCoords[1];
145                 ay = A.coords.usrCoords[2];
146             } else {
147                 ax = A[0];
148                 ay = A[1];
149             }
150 
151             if (B.coords) {
152                 bx = B.coords.usrCoords[1];
153                 by = B.coords.usrCoords[2];
154             } else {
155                 bx = B[0];
156                 by = B[1];
157             }
158 
159             if (C.coords) {
160                 cx = C.coords.usrCoords[1];
161                 cy = C.coords.usrCoords[2];
162             } else {
163                 cx = C[0];
164                 cy = C[1];
165             }
166 
167             phi = Math.atan2(cy - by, cx - bx) - Math.atan2(ay - by, ax - bx);
168 
169             if (phi < 0) {
170                 phi += 6.2831853071795862;
171             }
172 
173             return phi;
174         },
175 
176         /**
177          * Calculates a point on the bisection line between the three points A, B, C.
178          * As a result, the bisection line is defined by two points:
179          * Parameter B and the point with the coordinates calculated in this function.
180          * Does not work for ideal points.
181          * @param {JXG.Point} A Point
182          * @param {JXG.Point} B Point
183          * @param {JXG.Point} C Point
184          * @param [board=A.board] Reference to the board
185          * @returns {JXG.Coords} Coordinates of the second point defining the bisection.
186          */
187         angleBisector: function (A, B, C, board) {
188             var phiA,
189                 phiC,
190                 phi,
191                 Ac = A.coords.usrCoords,
192                 Bc = B.coords.usrCoords,
193                 Cc = C.coords.usrCoords,
194                 x,
195                 y;
196 
197             if (!Type.exists(board)) {
198                 board = A.board;
199             }
200 
201             // Parallel lines
202             if (Bc[0] === 0) {
203                 return new Coords(
204                     Const.COORDS_BY_USER,
205                     [1, (Ac[1] + Cc[1]) * 0.5, (Ac[2] + Cc[2]) * 0.5],
206                     board
207                 );
208             }
209 
210             // Non-parallel lines
211             x = Ac[1] - Bc[1];
212             y = Ac[2] - Bc[2];
213             phiA = Math.atan2(y, x);
214 
215             x = Cc[1] - Bc[1];
216             y = Cc[2] - Bc[2];
217             phiC = Math.atan2(y, x);
218 
219             phi = (phiA + phiC) * 0.5;
220 
221             if (phiA > phiC) {
222                 phi += Math.PI;
223             }
224 
225             x = Math.cos(phi) + Bc[1];
226             y = Math.sin(phi) + Bc[2];
227 
228             return new Coords(Const.COORDS_BY_USER, [1, x, y], board);
229         },
230 
231         // /**
232         //  * Calculates a point on the m-section line between the three points A, B, C.
233         //  * As a result, the m-section line is defined by two points:
234         //  * Parameter B and the point with the coordinates calculated in this function.
235         //  * The m-section generalizes the bisector to any real number.
236         //  * For example, the trisectors of an angle are simply the 1/3-sector and the 2/3-sector.
237         //  * Does not work for ideal points.
238         //  * @param {JXG.Point} A Point
239         //  * @param {JXG.Point} B Point
240         //  * @param {JXG.Point} C Point
241         //  * @param {Number} m Number
242         //  * @param [board=A.board] Reference to the board
243         //  * @returns {JXG.Coords} Coordinates of the second point defining the bisection.
244         //  */
245         // angleMsector: function (A, B, C, m, board) {
246         //     var phiA, phiC, phi,
247         //         Ac = A.coords.usrCoords,
248         //         Bc = B.coords.usrCoords,
249         //         Cc = C.coords.usrCoords,
250         //         x, y;
251 
252         //     if (!Type.exists(board)) {
253         //         board = A.board;
254         //     }
255 
256         //     // Parallel lines
257         //     if (Bc[0] === 0) {
258         //         return new Coords(Const.COORDS_BY_USER,
259         //             [1, (Ac[1] + Cc[1]) * m, (Ac[2] + Cc[2]) * m], board);
260         //     }
261 
262         //     // Non-parallel lines
263         //     x = Ac[1] - Bc[1];
264         //     y = Ac[2] - Bc[2];
265         //     phiA =  Math.atan2(y, x);
266 
267         //     x = Cc[1] - Bc[1];
268         //     y = Cc[2] - Bc[2];
269         //     phiC =  Math.atan2(y, x);
270 
271         //     phi = phiA + ((phiC - phiA) * m);
272 
273         //     if (phiA - phiC > Math.PI) {
274         //         phi += 2*m*Math.PI;
275         //     }
276 
277         //     x = Math.cos(phi) + Bc[1];
278         //     y = Math.sin(phi) + Bc[2];
279 
280         //     return new Coords(Const.COORDS_BY_USER, [1, x, y], board);
281         // },
282 
283         /**
284          * Reflects the point along the line.
285          * @param {JXG.Line} line Axis of reflection.
286          * @param {JXG.Point} point Point to reflect.
287          * @param [board=point.board] Reference to the board
288          * @returns {JXG.Coords} Coordinates of the reflected point.
289          */
290         reflection: function (line, point, board) {
291             // (v,w) defines the slope of the line
292             var x0,
293                 y0,
294                 x1,
295                 y1,
296                 v,
297                 w,
298                 mu,
299                 pc = point.coords.usrCoords,
300                 p1c = line.point1.coords.usrCoords,
301                 p2c = line.point2.coords.usrCoords;
302 
303             if (!Type.exists(board)) {
304                 board = point.board;
305             }
306 
307             v = p2c[1] - p1c[1];
308             w = p2c[2] - p1c[2];
309 
310             x0 = pc[1] - p1c[1];
311             y0 = pc[2] - p1c[2];
312 
313             mu = (v * y0 - w * x0) / (v * v + w * w);
314 
315             // point + mu*(-y,x) is the perpendicular foot
316             x1 = pc[1] + 2 * mu * w;
317             y1 = pc[2] - 2 * mu * v;
318 
319             return new Coords(Const.COORDS_BY_USER, [x1, y1], board);
320         },
321 
322         /**
323          * Computes the new position of a point which is rotated
324          * around a second point (called rotpoint) by the angle phi.
325          * @param {JXG.Point} rotpoint Center of the rotation
326          * @param {JXG.Point} point point to be rotated
327          * @param {Number} phi rotation angle in arc length
328          * @param {JXG.Board} [board=point.board] Reference to the board
329          * @returns {JXG.Coords} Coordinates of the new position.
330          */
331         rotation: function (rotpoint, point, phi, board) {
332             var x0,
333                 y0,
334                 c,
335                 s,
336                 x1,
337                 y1,
338                 pc = point.coords.usrCoords,
339                 rotpc = rotpoint.coords.usrCoords;
340 
341             if (!Type.exists(board)) {
342                 board = point.board;
343             }
344 
345             x0 = pc[1] - rotpc[1];
346             y0 = pc[2] - rotpc[2];
347 
348             c = Math.cos(phi);
349             s = Math.sin(phi);
350 
351             x1 = x0 * c - y0 * s + rotpc[1];
352             y1 = x0 * s + y0 * c + rotpc[2];
353 
354             return new Coords(Const.COORDS_BY_USER, [x1, y1], board);
355         },
356 
357         /**
358          * Calculates the coordinates of a point on the perpendicular to the given line through
359          * the given point.
360          * @param {JXG.Line} line A line.
361          * @param {JXG.Point} point Point which is projected to the line.
362          * @param {JXG.Board} [board=point.board] Reference to the board
363          * @returns {Array} Array of length two containing coordinates of a point on the perpendicular to the given line
364          *                  through the given point and boolean flag "change".
365          */
366         perpendicular: function (line, point, board) {
367             var x,
368                 y,
369                 change,
370                 c,
371                 z,
372                 A = line.point1.coords.usrCoords,
373                 B = line.point2.coords.usrCoords,
374                 C = point.coords.usrCoords;
375 
376             if (!Type.exists(board)) {
377                 board = point.board;
378             }
379 
380             // special case: point is the first point of the line
381             if (point === line.point1) {
382                 x = A[1] + B[2] - A[2];
383                 y = A[2] - B[1] + A[1];
384                 z = A[0] * B[0];
385 
386                 if (Math.abs(z) < Mat.eps) {
387                     x = B[2];
388                     y = -B[1];
389                 }
390                 c = [z, x, y];
391                 change = true;
392 
393                 // special case: point is the second point of the line
394             } else if (point === line.point2) {
395                 x = B[1] + A[2] - B[2];
396                 y = B[2] - A[1] + B[1];
397                 z = A[0] * B[0];
398 
399                 if (Math.abs(z) < Mat.eps) {
400                     x = A[2];
401                     y = -A[1];
402                 }
403                 c = [z, x, y];
404                 change = false;
405 
406                 // special case: point lies somewhere else on the line
407             } else if (Math.abs(Mat.innerProduct(C, line.stdform, 3)) < Mat.eps) {
408                 x = C[1] + B[2] - C[2];
409                 y = C[2] - B[1] + C[1];
410                 z = B[0];
411 
412                 if (Math.abs(z) < Mat.eps) {
413                     x = B[2];
414                     y = -B[1];
415                 }
416 
417                 change = true;
418                 if (
419                     Math.abs(z) > Mat.eps &&
420                     Math.abs(x - C[1]) < Mat.eps &&
421                     Math.abs(y - C[2]) < Mat.eps
422                 ) {
423                     x = C[1] + A[2] - C[2];
424                     y = C[2] - A[1] + C[1];
425                     change = false;
426                 }
427                 c = [z, x, y];
428 
429                 // general case: point does not lie on the line
430                 // -> calculate the foot of the dropped perpendicular
431             } else {
432                 c = [0, line.stdform[1], line.stdform[2]];
433                 c = Mat.crossProduct(c, C); // perpendicuar to line
434                 c = Mat.crossProduct(c, line.stdform); // intersection of line and perpendicular
435                 change = true;
436             }
437 
438             return [new Coords(Const.COORDS_BY_USER, c, board), change];
439         },
440 
441         /**
442          * @deprecated Please use {@link JXG.Math.Geometry.circumcenter} instead.
443          */
444         circumcenterMidpoint: function () {
445             JXG.deprecated("Geometry.circumcenterMidpoint()", "Geometry.circumcenter()");
446             this.circumcenter.apply(this, arguments);
447         },
448 
449         /**
450          * Calculates the center of the circumcircle of the three given points.
451          * @param {JXG.Point} point1 Point
452          * @param {JXG.Point} point2 Point
453          * @param {JXG.Point} point3 Point
454          * @param {JXG.Board} [board=point1.board] Reference to the board
455          * @returns {JXG.Coords} Coordinates of the center of the circumcircle of the given points.
456          */
457         circumcenter: function (point1, point2, point3, board) {
458             var u,
459                 v,
460                 m1,
461                 m2,
462                 A = point1.coords.usrCoords,
463                 B = point2.coords.usrCoords,
464                 C = point3.coords.usrCoords;
465 
466             if (!Type.exists(board)) {
467                 board = point1.board;
468             }
469 
470             u = [B[0] - A[0], -B[2] + A[2], B[1] - A[1]];
471             v = [(A[0] + B[0]) * 0.5, (A[1] + B[1]) * 0.5, (A[2] + B[2]) * 0.5];
472             m1 = Mat.crossProduct(u, v);
473 
474             u = [C[0] - B[0], -C[2] + B[2], C[1] - B[1]];
475             v = [(B[0] + C[0]) * 0.5, (B[1] + C[1]) * 0.5, (B[2] + C[2]) * 0.5];
476             m2 = Mat.crossProduct(u, v);
477 
478             return new Coords(Const.COORDS_BY_USER, Mat.crossProduct(m1, m2), board);
479         },
480 
481         /**
482          * Calculates the Euclidean distance for two given arrays of the same length.
483          * @param {Array} array1 Array of Number
484          * @param {Array} array2 Array of Number
485          * @param {Number} [n] Length of the arrays. Default is the minimum length of the given arrays.
486          * @returns {Number} Euclidean distance of the given vectors.
487          */
488         distance: function (array1, array2, n) {
489             var i,
490                 sum = 0;
491 
492             if (!n) {
493                 n = Math.min(array1.length, array2.length);
494             }
495 
496             for (i = 0; i < n; i++) {
497                 sum += (array1[i] - array2[i]) * (array1[i] - array2[i]);
498             }
499 
500             return Math.sqrt(sum);
501         },
502 
503         /**
504          * Calculates Euclidean distance for two given arrays of the same length.
505          * If one of the arrays contains a zero in the first coordinate, and the Euclidean distance
506          * is different from zero it is a point at infinity and we return Infinity.
507          * @param {Array} array1 Array containing elements of type number.
508          * @param {Array} array2 Array containing elements of type number.
509          * @param {Number} [n] Length of the arrays. Default is the minimum length of the given arrays.
510          * @returns {Number} Euclidean (affine) distance of the given vectors.
511          */
512         affineDistance: function (array1, array2, n) {
513             var d;
514 
515             d = this.distance(array1, array2, n);
516 
517             if (
518                 d > Mat.eps &&
519                 (Math.abs(array1[0]) < Mat.eps || Math.abs(array2[0]) < Mat.eps)
520             ) {
521                 return Infinity;
522             }
523 
524             return d;
525         },
526 
527         /**
528          * Affine ratio of three collinear points a, b, c: (c - a) / (b - a).
529          * If r > 1 or r < 0 then c is outside of the segment ab.
530          *
531          * @param {Array|JXG.Coords} a
532          * @param {Array|JXG.Coords} b
533          * @param {Array|JXG.Coords} c
534          * @returns {Number} affine ratio (c - a) / (b - a)
535          */
536         affineRatio: function (a, b, c) {
537             var r = 0.0,
538                 dx;
539 
540             if (Type.exists(a.usrCoords)) {
541                 a = a.usrCoords;
542             }
543             if (Type.exists(b.usrCoords)) {
544                 b = b.usrCoords;
545             }
546             if (Type.exists(c.usrCoords)) {
547                 c = c.usrCoords;
548             }
549 
550             dx = b[1] - a[1];
551 
552             if (Math.abs(dx) > Mat.eps) {
553                 r = (c[1] - a[1]) / dx;
554             } else {
555                 r = (c[2] - a[2]) / (b[2] - a[2]);
556             }
557             return r;
558         },
559 
560         /**
561          * Sort vertices counter clockwise starting with the first point.
562          * Used in Polygon.sutherlandHodgman, Geometry.signedPolygon.
563          *
564          * @param {Array} p An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
565          *
566          * @returns {Array}
567          */
568         sortVertices: function (p) {
569             var ll,
570                 ps = Expect.each(p, Expect.coordsArray),
571                 N = ps.length,
572                 lastPoint = null;
573 
574             // If the last point equals the first point, we take the last point out of the array.
575             // It may be that the several points at the end of the array are equal to the first point.
576             // The polygonal chain is been closed by JSXGraph, but this may also have been done by the user.
577             // Therefore, we use a while loop to pop the last points.
578             while (
579                 ps[0][0] === ps[N - 1][0] &&
580                 ps[0][1] === ps[N - 1][1] &&
581                 ps[0][2] === ps[N - 1][2]
582             ) {
583                 lastPoint = ps.pop();
584                 N--;
585             }
586 
587             ll = ps[0];
588             // Sort ps in increasing order of the angle between a point and the first point ll.
589             // If a point is equal to the first point ll, the angle is defined to be -Infinity.
590             // Otherwise, atan2 would return zero, which is a value which also attained by points
591             // on the same horizontal line.
592             ps.sort(function (a, b) {
593                 var rad1 =
594                         (a[2] === ll[2] && a[1] === ll[1])
595                             ? -Infinity
596                             : Math.atan2(a[2] - ll[2], a[1] - ll[1]),
597                     rad2 =
598                         (b[2] === ll[2] && b[1] === ll[1])
599                             ? -Infinity
600                             : Math.atan2(b[2] - ll[2], b[1] - ll[1]);
601                 return rad1 - rad2;
602             });
603 
604             // If the last point has been taken out of the array, we put it in again.
605             if (lastPoint !== null) {
606                 ps.push(lastPoint);
607             }
608 
609             return ps;
610         },
611 
612         /**
613          * Signed triangle area of the three points given. It can also be used
614          * to test the orientation of the triangle.
615          * <ul>
616          * <li> If the return value is < 0, then the point p2 is left of the line [p1, p3] (i.e p3 is right from [p1, p2]).
617          * <li> If the return value is > 0, then the point p2 is right of the line [p1, p3] (i.e p3 is left from [p1, p2]).
618          * <li> If the return value is = 0, then the points p1, p2, p3 are collinear.
619          * </ul>
620          *
621          * @param {JXG.Point|JXG.Coords|Array} p1
622          * @param {JXG.Point|JXG.Coords|Array} p2
623          * @param {JXG.Point|JXG.Coords|Array} p3
624          *
625          * @returns {Number}
626          */
627         signedTriangle: function (p1, p2, p3) {
628             var A = Expect.coordsArray(p1),
629                 B = Expect.coordsArray(p2),
630                 C = Expect.coordsArray(p3);
631             return 0.5 * ((B[1] - A[1]) * (C[2] - A[2]) - (B[2] - A[2]) * (C[1] - A[1]));
632         },
633 
634         /**
635          * Determine the signed area of a non-self-intersecting polygon.
636          * Surveyor's Formula
637          *
638          * @param {Array} p An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
639          * @param {Boolean} [sort=true]
640          *
641          * @returns {Number}
642          */
643         signedPolygon: function (p, sort) {
644             var i,
645                 N,
646                 A = 0,
647                 ps = Expect.each(p, Expect.coordsArray);
648 
649             if (sort === undefined) {
650                 sort = true;
651             }
652 
653             if (!sort) {
654                 ps = this.sortVertices(ps);
655             } else {
656                 // Make sure the polygon is closed. If it is already closed this won't change the sum because the last
657                 // summand will be 0.
658                 ps.unshift(ps[ps.length - 1]);
659             }
660 
661             N = ps.length;
662 
663             for (i = 1; i < N; i++) {
664                 A += ps[i - 1][1] * ps[i][2] - ps[i][1] * ps[i - 1][2];
665             }
666 
667             return 0.5 * A;
668         },
669 
670         /**
671          * Calculate the complex hull of a point cloud by the Graham scan algorithm.
672          *
673          * @param {Array} points An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
674          *
675          * @returns {Array} List of objects <pre>{i: index, c: coords}</pre> containing the convex hull points
676          *  in form of the index in the original input array and a coords array.
677          *
678          * @example
679          *     // Static example
680          *
681          *     var i, hull,
682          *       p = [],
683          *       q = [];
684          *
685          *     p.push( board.create('point', [4, 0], {withLabel:false }) );
686          *     p.push( board.create('point', [0, 4], {withLabel:false }) );
687          *     p.push( board.create('point', [0, 0], {withLabel:false }) );
688          *     p.push([-1, 0]);
689          *     p.push([-3, -3]);
690          *
691          *     hull = JXG.Math.Geometry.GrahamScan(p);
692          *     for (i = 0; i < hull.length; i++) {
693          *       console.log(hull[i]);
694          *       q.push(hull[i].c);
695          *     }
696          *     board.create('polygon', q);
697          *     // Output:
698          *     // { i: 4, c: [1, -3, 3]}
699          *     // { i: 0, c: [1, 4, 0]}
700          *     // { i: 1, c: [1, 0, 4]}
701          *
702          * </pre><div id="JXGb310b874-595e-4020-b0c2-566482797836" class="jxgbox" style="width: 300px; height: 300px;"></div>
703          * <script type="text/javascript">
704          *     (function() {
705          *         var board = JXG.JSXGraph.initBoard('JXGb310b874-595e-4020-b0c2-566482797836',
706          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
707          *         var i, hull,
708          *           p = [],
709          *           q = [];
710          *
711          *         p.push( board.create('point', [4, 0], {withLabel:false }) );
712          *         p.push( board.create('point', [0, 4], {withLabel:false }) );
713          *         p.push( board.create('point', [0, 0], {withLabel:false }) );
714          *         p.push([-1, 0]);
715          *         p.push([-3, -3]);
716          *
717          *         hull = JXG.Math.Geometry.GrahamScan(p);
718          *         for (i = 0; i < hull.length; i++) {
719          *           console.log(hull[i]);
720          *           q.push(hull[i].c);
721          *         }
722          *         board.create('polygon', q);
723          *
724          *     })();
725          *
726          * </script><pre>
727          *
728          */
729         GrahamScan: function (points) {
730             var i, M, o,
731                 mi_idx,
732                 mi_x, mi_y, ma_x, ma_y,
733                 mi_xpy, mi_xmy, ma_xpy, ma_xmy,
734                 mi_x_i, ma_x_i, mi_y_i, ma_y_i,
735                 mi_xpy_i, mi_xmy_i, ma_xpy_i, ma_xmy_i,
736                 v, c,
737                 eps = Mat.eps * Mat.eps,
738                 that = this,
739                 ps_idx = [],
740                 stack = [],
741                 ps = Expect.each(points, Expect.coordsArray), // New array object, i.e. a copy of the input array.
742                 N,
743                 AklToussaint = 1024;  // This is a rough threshold where the heuristic pays off.
744 
745             N = ps.length;
746             if (N === 0) {
747                 return [];
748             }
749 
750             if (N > AklToussaint) {
751                 //
752                 // Akl-Toussaint heuristic
753                 // Determine an irregular convex octagon whose inside can be discarded.
754                 //
755                 mi_x = ps[0][1];
756                 ma_x = mi_x;
757                 mi_y = ps[0][2];
758                 ma_y = mi_y;
759 
760                 mi_xmy = ps[0][1] - ps[0][2];
761                 ma_xmy = mi_xmy;
762                 mi_xpy = ps[0][1] + ps[0][2];
763                 ma_xpy = mi_xpy;
764 
765                 mi_x_i = 0;
766                 ma_x_i = 0;
767                 mi_y_i = 0;
768                 ma_y_i = 0;
769 
770                 mi_xmy_i = 0;
771                 ma_xmy_i = 0;
772                 mi_xpy_i = 0;
773                 ma_xpy_i = 0;
774                 for (i = 1; i < N; i++) {
775                     v = ps[i][1];
776                     if (v < mi_x) {
777                         mi_x = v;
778                         mi_x_i = i;
779                     } else if (v > ma_x) {
780                         ma_x = v;
781                         ma_x_i = i;
782                     }
783 
784                     v = ps[i][2];
785                     if (v < mi_y) {
786                         mi_y = v;
787                         mi_y_i = i;
788                     } else if (v > ma_y) {
789                         ma_y = v;
790                         ma_y_i = i;
791                     }
792 
793                     v = ps[i][1] - ps[i][2];
794                     if (v < mi_xmy) {
795                         mi_xmy = v;
796                         mi_xmy_i = i;
797                     } else if (v > ma_xmy) {
798                         ma_xmy = v;
799                         ma_xmy_i = i;
800                     }
801 
802                     v = ps[i][1] + ps[i][2];
803                     if (v < mi_xpy) {
804                         mi_xpy = v;
805                         mi_xpy_i = i;
806                     } else if (v > ma_xpy) {
807                         ma_xpy = v;
808                         ma_xpy_i = i;
809                     }
810                 }
811             }
812 
813             // Keep track of the indices of the input points.
814             for (i = 0; i < N; i++) {
815                 c = ps[i];
816                 if (N <= AklToussaint ||
817                     // Discard inside of the octagon according to the Akl-Toussaint heuristic
818                     i in [mi_x_i, ma_x_i, mi_y_i, ma_y_i, mi_xpy_i, mi_xmy_i, ma_xpy_i, ma_xmy_i] ||
819                     (mi_x_i !== mi_xmy_i && this.signedTriangle(ps[mi_x_i], ps[mi_xmy_i], c) >= -eps) ||
820                     (mi_xmy_i !== ma_y_i && this.signedTriangle(ps[mi_xmy_i], ps[ma_y_i], c) >= -eps) ||
821                     (ma_y_i !== ma_xpy_i && this.signedTriangle(ps[ma_y_i], ps[ma_xpy_i], c) >= -eps) ||
822                     (ma_xpy_i !== ma_x_i && this.signedTriangle(ps[ma_xpy_i], ps[ma_x_i], c) >= -eps) ||
823                     (ma_x_i !== ma_xmy_i && this.signedTriangle(ps[ma_x_i], ps[ma_xmy_i], c) >= -eps) ||
824                     (ma_xmy_i !== mi_y_i && this.signedTriangle(ps[ma_xmy_i], ps[mi_y_i], c) >= -eps) ||
825                     (mi_y_i !== mi_xpy_i && this.signedTriangle(ps[mi_y_i], ps[mi_xpy_i], c) >= -eps) ||
826                     (mi_xpy_i !== mi_x_i && this.signedTriangle(ps[mi_xpy_i], ps[mi_x_i], c) >= -eps)
827                 ) {
828                     ps_idx.push({
829                         i: i,
830                         c: c
831                     });
832                 }
833             }
834             N = ps_idx.length;
835 
836             // Find the point with the lowest y value
837             mi_idx = 0;
838             mi_x = ps_idx[0].c[1];
839             mi_y = ps_idx[0].c[2];
840             for (i = 1; i < N; i++) {
841                 if ((ps_idx[i].c[2] < mi_y) || (ps_idx[i].c[2] === mi_y && ps_idx[i].c[1] < mi_x)) {
842                     mi_x = ps_idx[i].c[1];
843                     mi_y = ps_idx[i].c[2];
844                     mi_idx = i;
845                 }
846             }
847             ps_idx = Type.swap(ps_idx, mi_idx, 0);
848 
849             // Our origin o, i.e. the first point.
850             o = ps_idx[0].c;
851 
852             // Sort according to the angle around o.
853             ps_idx.sort(function(a_obj, b_obj) {
854                 var a = a_obj.c,
855                     b = b_obj.c,
856                     v = that.signedTriangle(o, a, b);
857 
858                 if (v === 0) {
859                     // if o, a, b are collinear, the point which is further away
860                     // from o is considered greater.
861                     return Mat.hypot(a[1] - o[1], a[2] - o[2]) - Mat.hypot(b[1] - o[1], b[2] - o[2]);
862                 }
863 
864                 // if v < 0, a is to the left of [o, b], i.e. angle(a) > angle(b)
865                 return -v;
866             });
867 
868             // Do the Graham scan.
869             M = 0;
870             for (i = 0; i < N; i++) {
871                 while (M > 1 && this.signedTriangle(stack[M - 2].c, stack[M - 1].c, ps_idx[i].c) <= 0) {
872                     // stack[M - 1] is to the left of stack[M - 1], ps[i]: discard it
873                     stack.pop();
874                     M--;
875                 }
876                 stack.push(ps_idx[i]);
877                 M++;
878             }
879 
880             return stack;
881         },
882 
883         // Original method
884         // GrahamScan: function (points, indices) {
885         //     var i,
886         //         M = 1,
887         //         ps = Expect.each(points, Expect.coordsArray),
888         //         N = ps.length;
889         //     ps = this.sortVertices(ps);
890         //     N = ps.length;
891         //     for (i = 2; i < N; i++) {
892         //         while (this.signedTriangle(ps[M - 1], ps[M], ps[i]) <= 0) {
893         //             if (M > 1) {
894         //                 M -= 1;
895         //             } else if (i === N - 1) {
896         //                 break;
897         //             }
898         //             i += 1;
899         //         }
900         //         M += 1;
901         //         ps = Type.swap(ps, M, i);
902         //         indices = Type.swap(indices, M, i);
903         //     }
904         //     return ps.slice(0, M);
905         // },
906 
907         /**
908          * Calculate the complex hull of a point cloud by the Graham scan algorithm.
909          *
910          * @param {Array} points An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
911          * @param {Boolean} [returnCoords=false] If true, return an array of coords. Otherwise return a list of pointers
912          * to the input list elements. That is, if the input is a list of {@link JXG.Point} elements, the returned list
913          * will contain the points that form the convex hull.
914          * @returns {Array} List containing the convex hull. Format depends on returnCoords.
915          * @see JXG.Math.Geometry.GrahamScan
916          *
917          * @example
918          *     // Static example
919          *     var i, hull,
920          *         p = [];
921          *
922          *     p.push( board.create('point', [4, 0], {withLabel:false }) );
923          *     p.push( board.create('point', [0, 4], {withLabel:false }) );
924          *     p.push( board.create('point', [0, 0], {withLabel:false }) );
925          *     p.push( board.create('point', [1, 1], {withLabel:false }) );
926          *     hull = JXG.Math.Geometry.convexHull(p);
927          *     for (i = 0; i < hull.length; i++) {
928          *       hull[i].setAttribute({color: 'blue'});
929          *     }
930          *
931          * </pre><div id="JXGdfc76123-81b8-4250-96f9-419253bd95dd" class="jxgbox" style="width: 300px; height: 300px;"></div>
932          * <script type="text/javascript">
933          *     (function() {
934          *         var board = JXG.JSXGraph.initBoard('JXGdfc76123-81b8-4250-96f9-419253bd95dd',
935          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
936          *         var i, hull,
937          *             p = [];
938          *
939          *         p.push( board.create('point', [4, 0], {withLabel:false }) );
940          *         p.push( board.create('point', [0, 4], {withLabel:false }) );
941          *         p.push( board.create('point', [0, 0], {withLabel:false }) );
942          *         p.push( board.create('point', [1, 1], {withLabel:false }) );
943          *         hull = JXG.Math.Geometry.convexHull(p);
944          *         for (i = 0; i < hull.length; i++) {
945          *           hull[i].setAttribute({color: 'blue'});
946          *         }
947          *
948          *     })();
949          *
950          * </script><pre>
951          *
952          * @example
953          *     // Dynamic version using returnCoords==true: drag the points
954          *     var p = [];
955          *
956          *     p.push( board.create('point', [4, 0], {withLabel:false }) );
957          *     p.push( board.create('point', [0, 4], {withLabel:false }) );
958          *     p.push( board.create('point', [0, 0], {withLabel:false }) );
959          *     p.push( board.create('point', [1, 1], {withLabel:false }) );
960          *
961          *     var c = board.create('curve', [[], []], {fillColor: 'yellow', fillOpacity: 0.3});
962          *     c.updateDataArray = function() {
963          *       var i,
964          *         hull = JXG.Math.Geometry.convexHull(p, true);
965          *
966          *       this.dataX = [];
967          *       this.dataY = [];
968          *
969          *       for (i = 0; i < hull.length; i ++) {
970          *         this.dataX.push(hull[i][1]);
971          *         this.dataY.push(hull[i][2]);
972          *       }
973          *       this.dataX.push(hull[0][1]);
974          *       this.dataY.push(hull[0][2]);
975          *     };
976          *     board.update();
977          *
978          * </pre><div id="JXG61e51909-da0b-432f-9aa7-9fb0c8bb01c9" class="jxgbox" style="width: 300px; height: 300px;"></div>
979          * <script type="text/javascript">
980          *     (function() {
981          *         var board = JXG.JSXGraph.initBoard('JXG61e51909-da0b-432f-9aa7-9fb0c8bb01c9',
982          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
983          *         var p = [];
984          *
985          *         p.push( board.create('point', [4, 0], {withLabel:false }) );
986          *         p.push( board.create('point', [0, 4], {withLabel:false }) );
987          *         p.push( board.create('point', [0, 0], {withLabel:false }) );
988          *         p.push( board.create('point', [1, 1], {withLabel:false }) );
989          *
990          *         var c = board.create('curve', [[], []], {fillColor: 'yellow', fillOpacity: 0.3});
991          *         c.updateDataArray = function() {
992          *           var i,
993          *             hull = JXG.Math.Geometry.convexHull(p, true);
994          *
995          *           this.dataX = [];
996          *           this.dataY = [];
997          *
998          *           for (i = 0; i < hull.length; i ++) {
999          *             this.dataX.push(hull[i][1]);
1000          *             this.dataY.push(hull[i][2]);
1001          *           }
1002          *           this.dataX.push(hull[0][1]);
1003          *           this.dataY.push(hull[0][2]);
1004          *         };
1005          *         board.update();
1006          *
1007          *
1008          *     })();
1009          *
1010          * </script><pre>
1011          *
1012          */
1013         convexHull: function(points, returnCoords) {
1014             var i, hull,
1015                 res = [];
1016 
1017             hull = this.GrahamScan(points);
1018             for (i = 0; i < hull.length; i++) {
1019                 if (returnCoords) {
1020                     res.push(hull[i].c);
1021                 } else {
1022                     res.push(points[hull[i].i]);
1023                 }
1024             }
1025             return res;
1026         },
1027 
1028         // /**
1029         //  * Determine if a polygon or a path element is convex, non-convex or complex which are defined like this:
1030         //  * <ul>
1031         //  * <li> A polygon is convex if for every pair of points, the line segment connecting them does not intersect
1032         //  * an edge of the polygon in one point.
1033         //  * A single line segment or a a single point is considered as convex. A necessary condition for a polygon
1034         //  * to be convex that the angle sum of its interior angles equals ± 2 π.
1035         //  * <li> A polygon is non-convex, if it does not self-intersect, but is not convex.
1036         //  * <li> A polygon is complex if its the angle sum is not equal to ± 2 π.
1037         //  * That is, there must be self-intersection (contiguous coincident points in the path are not treated as self-intersection).
1038         //  * </ul>
1039         //  * A path  element might be specified as an array of coordinate arrays or {@link JXG.Coords}.
1040         //  *
1041         //  * @param {Array|Polygon|PolygonalChain} points Polygon or list of coordinates
1042         //  * @returns {Number} -1: if complex, 0: if non-convex, 1: if convex
1043         //  */
1044         /**
1045          * Determine if a polygon or a path element is convex:
1046          * <p>
1047          * A polygon is convex if for every pair of points, the line segment connecting them does not intersect
1048          * an edge of the polygon in one point.
1049          * A single line segment, a single point, or the empty set is considered as convex. A necessary condition for a polygon
1050          * to be convex that the angle sum of its interior angles equals ± 2 π.
1051          * <p>
1052          * A path  element might be specified as an array of coordinate arrays or {@link JXG.Coords}.
1053          * See the discussion at <a href="https://stackoverflow.com/questions/471962/how-do-i-efficiently-determine-if-a-polygon-is-convex-non-convex-or-complex">stackoverflow</a>.
1054          *
1055          * @param {Array|Polygon|PolygonalChain} points Polygon or list of coordinates
1056          * @returns {Boolean} true if convex
1057          *
1058          * @example
1059          * var pol = board.create('polygon', [
1060          *     [-1, -1],
1061          *     [3, -1],
1062          *     [4, 2],
1063          *     [3, 3],
1064          *     [0, 4],
1065          *     [-3, 1]
1066          * ], {
1067          *     vertices: {
1068          *         color: 'blue',
1069          *         snapToGrid: true
1070          *     }
1071          * });
1072          *
1073          * console.log(JXG.Math.Geometry.isConvex(pol));
1074          * // > true
1075          *
1076          *
1077          *
1078          * </pre><div id="JXG9b43cc53-15b4-49be-92cc-2a1dfc06665b" class="jxgbox" style="width: 300px; height: 300px;"></div>
1079          * <script type="text/javascript">
1080          *     (function() {
1081          *         var board = JXG.JSXGraph.initBoard('JXG9b43cc53-15b4-49be-92cc-2a1dfc06665b',
1082          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1083          *     var pol = board.create('polygon', [
1084          *         [-1, -1],
1085          *         [3, -1],
1086          *         [4, 2],
1087          *         [3, 3],
1088          *         [0, 4],
1089          *         [-3, 1]
1090          *     ], {
1091          *         vertices: {
1092          *             color: 'blue',
1093          *             snapToGrid: true
1094          *         }
1095          *     });
1096          *
1097          *     console.log(JXG.Math.Geometry.isConvex(pol));
1098          *
1099          *
1100          *
1101          *     })();
1102          *
1103          * </script><pre>
1104          *
1105          */
1106         isConvex: function(points) {
1107             var ps, le, i,
1108                 eps = Mat.eps * Mat.eps,
1109                 old_x, old_y, old_dir,
1110                 new_x, new_y, new_dir,
1111                 angle,
1112                 orient,
1113                 angle_sum = 0.0;
1114 
1115             if (Type.isArray(points)) {
1116                 ps = Expect.each(points, Expect.coordsArray);
1117             } else if (Type.exists(points.type) && points.type === Const.OBJECT_TYPE_POLYGON) {
1118                 ps = Expect.each(points.vertices, Expect.coordsArray);
1119             }
1120             le = ps.length;
1121             if (le === 0) {
1122                 // Empty set is convex
1123                 return true;
1124             }
1125             if (le < 3) {
1126                 // Segments and points are convex
1127                 return true;
1128             }
1129 
1130             orient = null;
1131             old_x = ps[le - 2][1];
1132             old_y = ps[le - 2][2];
1133             new_x = ps[le - 1][1];
1134             new_y = ps[le - 1][2];
1135             new_dir = Math.atan2(new_y - old_y, new_x - old_x);
1136             for (i = 0; i < le; i++) {
1137                 old_x = new_x;
1138                 old_y = new_y;
1139                 old_dir = new_dir;
1140                 new_x = ps[i][1];
1141                 new_y = ps[i][2];
1142                 if (old_x === new_x && old_y === new_y) {
1143                     // Repeated consecutive points are ignored
1144                     continue;
1145                 }
1146                 new_dir = Math.atan2(new_y - old_y, new_x - old_x);
1147                 angle = new_dir - old_dir;
1148                 if (angle <= -Math.PI) {
1149                     angle += 2 * Math.PI;
1150                 } else if (angle > Math.PI) {
1151                     angle -= 2 * Math.PI;
1152                 }
1153                 if (orient === null) {
1154                     if (angle === 0.0) {
1155                         continue;
1156                     }
1157                     orient = (angle > 0) ? 1 : -1;
1158                 } else {
1159                     if (orient * angle < -eps) {
1160                         return false;
1161                     }
1162                 }
1163                 angle_sum += angle;
1164             }
1165 
1166             if ((Math.abs(angle_sum / (2 * Math.PI)) - 1) < eps) {
1167                 return true;
1168             }
1169             return false;
1170         },
1171 
1172         /**
1173          * A line can be a segment, a straight, or a ray. So it is not always delimited by point1 and point2
1174          * calcStraight determines the visual start point and end point of the line. A segment is only drawn
1175          * from start to end point, a straight line is drawn until it meets the boards boundaries.
1176          * @param {JXG.Line} el Reference to a line object, that needs calculation of start and end point.
1177          * @param {JXG.Coords} point1 Coordinates of the point where line drawing begins. This value is calculated and
1178          * set by this method.
1179          * @param {JXG.Coords} point2 Coordinates of the point where line drawing ends. This value is calculated and set
1180          * by this method.
1181          * @param {Number} margin Optional margin, to avoid the display of the small sides of lines.
1182          * @returns null
1183          * @see Line
1184          * @see JXG.Line
1185          */
1186         calcStraight: function (el, point1, point2, margin) {
1187             var takePoint1,
1188                 takePoint2,
1189                 intersection,
1190                 intersect1,
1191                 intersect2,
1192                 straightFirst,
1193                 straightLast,
1194                 c, p1, p2;
1195 
1196             if (!Type.exists(margin)) {
1197                 // Enlarge the drawable region slightly. This hides the small sides
1198                 // of thick lines in most cases.
1199                 margin = 10;
1200             }
1201 
1202             straightFirst = el.evalVisProp('straightfirst');
1203             straightLast = el.evalVisProp('straightlast');
1204 
1205             // If one of the point is an ideal point in homogeneous coordinates
1206             // drawing of line segments or rays are not possible.
1207             if (Math.abs(point1.scrCoords[0]) < Mat.eps) {
1208                 straightFirst = true;
1209             }
1210             if (Math.abs(point2.scrCoords[0]) < Mat.eps) {
1211                 straightLast = true;
1212             }
1213 
1214             // Do nothing in case of line segments (inside or outside of the board)
1215             if (!straightFirst && !straightLast) {
1216                 return;
1217             }
1218 
1219             // Compute the stdform of the line in screen coordinates.
1220             c = [];
1221             c[0] =
1222                 el.stdform[0] -
1223                 (el.stdform[1] * el.board.origin.scrCoords[1]) / el.board.unitX +
1224                 (el.stdform[2] * el.board.origin.scrCoords[2]) / el.board.unitY;
1225             c[1] = el.stdform[1] / el.board.unitX;
1226             c[2] = -el.stdform[2] / el.board.unitY;
1227 
1228             // If p1=p2
1229             if (isNaN(c[0] + c[1] + c[2])) {
1230                 return;
1231             }
1232 
1233             takePoint1 = false;
1234             takePoint2 = false;
1235 
1236             // Line starts at point1 and point1 is inside the board
1237             takePoint1 =
1238                 !straightFirst &&
1239                 Math.abs(point1.usrCoords[0]) >= Mat.eps &&
1240                 point1.scrCoords[1] >= 0.0 &&
1241                 point1.scrCoords[1] <= el.board.canvasWidth &&
1242                 point1.scrCoords[2] >= 0.0 &&
1243                 point1.scrCoords[2] <= el.board.canvasHeight;
1244 
1245             // Line ends at point2 and point2 is inside the board
1246             takePoint2 =
1247                 !straightLast &&
1248                 Math.abs(point2.usrCoords[0]) >= Mat.eps &&
1249                 point2.scrCoords[1] >= 0.0 &&
1250                 point2.scrCoords[1] <= el.board.canvasWidth &&
1251                 point2.scrCoords[2] >= 0.0 &&
1252                 point2.scrCoords[2] <= el.board.canvasHeight;
1253 
1254             // Intersect the line with the four borders of the board.
1255             intersection = this.meetLineBoard(c, el.board, margin);
1256             intersect1 = intersection[0];
1257             intersect2 = intersection[1];
1258 
1259             /**
1260              * At this point we have four points:
1261              * point1 and point2 are the first and the second defining point on the line,
1262              * intersect1, intersect2 are the intersections of the line with border around the board.
1263              */
1264 
1265             /*
1266              * Here we handle rays where both defining points are outside of the board.
1267              */
1268             // If both points are outside and the complete ray is outside we do nothing
1269             if (!takePoint1 && !takePoint2) {
1270                 // Ray starting at point 1
1271                 if (
1272                     !straightFirst &&
1273                     straightLast &&
1274                     !this.isSameDirection(point1, point2, intersect1) &&
1275                     !this.isSameDirection(point1, point2, intersect2)
1276                 ) {
1277                     return;
1278                 }
1279 
1280                 // Ray starting at point 2
1281                 if (
1282                     straightFirst &&
1283                     !straightLast &&
1284                     !this.isSameDirection(point2, point1, intersect1) &&
1285                     !this.isSameDirection(point2, point1, intersect2)
1286                 ) {
1287                     return;
1288                 }
1289             }
1290 
1291             /*
1292              * If at least one of the defining points is outside of the board
1293              * we take intersect1 or intersect2 as one of the end points
1294              * The order is also important for arrows of axes
1295              */
1296             if (!takePoint1) {
1297                 if (!takePoint2) {
1298                     // Two border intersection points are used
1299                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1300                         p1 = intersect1;
1301                         p2 = intersect2;
1302                     } else {
1303                         p2 = intersect1;
1304                         p1 = intersect2;
1305                     }
1306                 } else {
1307                     // One border intersection points is used
1308                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1309                         p1 = intersect1;
1310                     } else {
1311                         p1 = intersect2;
1312                     }
1313                 }
1314             } else {
1315                 if (!takePoint2) {
1316                     // One border intersection points is used
1317                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1318                         p2 = intersect2;
1319                     } else {
1320                         p2 = intersect1;
1321                     }
1322                 }
1323             }
1324 
1325             if (p1) {
1326                 //point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords.slice(1));
1327                 point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords);
1328             }
1329 
1330             if (p2) {
1331                 //point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords.slice(1));
1332                 point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords);
1333             }
1334         },
1335 
1336         /**
1337          * A line can be a segment, a straight, or a ray. so it is not always delimited by point1 and point2.
1338          *
1339          * This method adjusts the line's delimiting points taking into account its nature, the viewport defined
1340          * by the board.
1341          *
1342          * A segment is delimited by start and end point, a straight line or ray is delimited until it meets the
1343          * boards boundaries. However, if the line has infinite ticks, it will be delimited by the projection of
1344          * the boards vertices onto itself.
1345          *
1346          * @param {JXG.Line} el Reference to a line object, that needs calculation of start and end point.
1347          * @param {JXG.Coords} point1 Coordinates of the point where line drawing begins. This value is calculated and
1348          * set by this method.
1349          * @param {JXG.Coords} point2 Coordinates of the point where line drawing ends. This value is calculated and set
1350          * by this method.
1351          * @see Line
1352          * @see JXG.Line
1353          */
1354         calcLineDelimitingPoints: function (el, point1, point2) {
1355             var distP1P2,
1356                 boundingBox,
1357                 lineSlope,
1358                 intersect1,
1359                 intersect2,
1360                 straightFirst,
1361                 straightLast,
1362                 c,
1363                 p1,
1364                 p2,
1365                 takePoint1 = false,
1366                 takePoint2 = false;
1367 
1368             straightFirst = el.evalVisProp('straightfirst');
1369             straightLast = el.evalVisProp('straightlast');
1370 
1371             // If one of the point is an ideal point in homogeneous coordinates
1372             // drawing of line segments or rays are not possible.
1373             if (Math.abs(point1.scrCoords[0]) < Mat.eps) {
1374                 straightFirst = true;
1375             }
1376             if (Math.abs(point2.scrCoords[0]) < Mat.eps) {
1377                 straightLast = true;
1378             }
1379 
1380             // Compute the stdform of the line in screen coordinates.
1381             c = [];
1382             c[0] =
1383                 el.stdform[0] -
1384                 (el.stdform[1] * el.board.origin.scrCoords[1]) / el.board.unitX +
1385                 (el.stdform[2] * el.board.origin.scrCoords[2]) / el.board.unitY;
1386             c[1] = el.stdform[1] / el.board.unitX;
1387             c[2] = -el.stdform[2] / el.board.unitY;
1388 
1389             // p1=p2
1390             if (isNaN(c[0] + c[1] + c[2])) {
1391                 return;
1392             }
1393 
1394             takePoint1 = !straightFirst;
1395             takePoint2 = !straightLast;
1396             // Intersect the board vertices on the line to establish the available visual space for the infinite ticks
1397             // Based on the slope of the line we can optimise and only project the two outer vertices
1398 
1399             // boundingBox = [x1, y1, x2, y2] upper left, lower right vertices
1400             boundingBox = el.board.getBoundingBox();
1401             lineSlope = el.getSlope();
1402             if (lineSlope >= 0) {
1403                 // project vertices (x2,y1) (x1, y2)
1404                 intersect1 = this.projectPointToLine(
1405                     { coords: { usrCoords: [1, boundingBox[2], boundingBox[1]] } },
1406                     el,
1407                     el.board
1408                 );
1409                 intersect2 = this.projectPointToLine(
1410                     { coords: { usrCoords: [1, boundingBox[0], boundingBox[3]] } },
1411                     el,
1412                     el.board
1413                 );
1414             } else {
1415                 // project vertices (x1, y1) (x2, y2)
1416                 intersect1 = this.projectPointToLine(
1417                     { coords: { usrCoords: [1, boundingBox[0], boundingBox[1]] } },
1418                     el,
1419                     el.board
1420                 );
1421                 intersect2 = this.projectPointToLine(
1422                     { coords: { usrCoords: [1, boundingBox[2], boundingBox[3]] } },
1423                     el,
1424                     el.board
1425                 );
1426             }
1427 
1428             /**
1429              * we have four points:
1430              * point1 and point2 are the first and the second defining point on the line,
1431              * intersect1, intersect2 are the intersections of the line with border around the board.
1432              */
1433 
1434             /*
1435              * Here we handle rays/segments where both defining points are outside of the board.
1436              */
1437             if (!takePoint1 && !takePoint2) {
1438                 // Segment, if segment does not cross the board, do nothing
1439                 if (!straightFirst && !straightLast) {
1440                     distP1P2 = point1.distance(Const.COORDS_BY_USER, point2);
1441                     // if  intersect1 not between point1 and point2
1442                     if (
1443                         Math.abs(
1444                             point1.distance(Const.COORDS_BY_USER, intersect1) +
1445                             intersect1.distance(Const.COORDS_BY_USER, point2) -
1446                             distP1P2
1447                         ) > Mat.eps
1448                     ) {
1449                         return;
1450                     }
1451                     // if insersect2 not between point1 and point2
1452                     if (
1453                         Math.abs(
1454                             point1.distance(Const.COORDS_BY_USER, intersect2) +
1455                             intersect2.distance(Const.COORDS_BY_USER, point2) -
1456                             distP1P2
1457                         ) > Mat.eps
1458                     ) {
1459                         return;
1460                     }
1461                 }
1462 
1463                 // If both points are outside and the complete ray is outside we do nothing
1464                 // Ray starting at point 1
1465                 if (
1466                     !straightFirst &&
1467                     straightLast &&
1468                     !this.isSameDirection(point1, point2, intersect1) &&
1469                     !this.isSameDirection(point1, point2, intersect2)
1470                 ) {
1471                     return;
1472                 }
1473 
1474                 // Ray starting at point 2
1475                 if (
1476                     straightFirst &&
1477                     !straightLast &&
1478                     !this.isSameDirection(point2, point1, intersect1) &&
1479                     !this.isSameDirection(point2, point1, intersect2)
1480                 ) {
1481                     return;
1482                 }
1483             }
1484 
1485             /*
1486              * If at least one of the defining points is outside of the board
1487              * we take intersect1 or intersect2 as one of the end points
1488              * The order is also important for arrows of axes
1489              */
1490             if (!takePoint1) {
1491                 if (!takePoint2) {
1492                     // Two border intersection points are used
1493                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1494                         p1 = intersect1;
1495                         p2 = intersect2;
1496                     } else {
1497                         p2 = intersect1;
1498                         p1 = intersect2;
1499                     }
1500                 } else {
1501                     // One border intersection points is used
1502                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1503                         p1 = intersect1;
1504                     } else {
1505                         p1 = intersect2;
1506                     }
1507                 }
1508             } else {
1509                 if (!takePoint2) {
1510                     // One border intersection points is used
1511                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
1512                         p2 = intersect2;
1513                     } else {
1514                         p2 = intersect1;
1515                     }
1516                 }
1517             }
1518 
1519             if (p1) {
1520                 //point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords.slice(1));
1521                 point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords);
1522             }
1523 
1524             if (p2) {
1525                 //point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords.slice(1));
1526                 point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords);
1527             }
1528         },
1529 
1530         /**
1531          * Calculates the visProp.position corresponding to a given angle.
1532          * @param {number} angle angle in radians. Must be in range (-2pi,2pi).
1533          */
1534         calcLabelQuadrant: function (angle) {
1535             var q;
1536             if (angle < 0) {
1537                 angle += 2 * Math.PI;
1538             }
1539             q = Math.floor((angle + Math.PI / 8) / (Math.PI / 4)) % 8;
1540             return ["rt", "urt", "top", "ulft", "lft", "llft", "lrt"][q];
1541         },
1542 
1543         /**
1544          * The vectors <tt>p2-p1</tt> and <tt>i2-i1</tt> are supposed to be collinear. If their cosine is positive
1545          * they point into the same direction otherwise they point in opposite direction.
1546          * @param {JXG.Coords} p1
1547          * @param {JXG.Coords} p2
1548          * @param {JXG.Coords} i1
1549          * @param {JXG.Coords} i2
1550          * @returns {Boolean} True, if <tt>p2-p1</tt> and <tt>i2-i1</tt> point into the same direction
1551          */
1552         isSameDir: function (p1, p2, i1, i2) {
1553             var dpx = p2.usrCoords[1] - p1.usrCoords[1],
1554                 dpy = p2.usrCoords[2] - p1.usrCoords[2],
1555                 dix = i2.usrCoords[1] - i1.usrCoords[1],
1556                 diy = i2.usrCoords[2] - i1.usrCoords[2];
1557 
1558             if (Math.abs(p2.usrCoords[0]) < Mat.eps) {
1559                 dpx = p2.usrCoords[1];
1560                 dpy = p2.usrCoords[2];
1561             }
1562 
1563             if (Math.abs(p1.usrCoords[0]) < Mat.eps) {
1564                 dpx = -p1.usrCoords[1];
1565                 dpy = -p1.usrCoords[2];
1566             }
1567 
1568             return dpx * dix + dpy * diy >= 0;
1569         },
1570 
1571         /**
1572          * If you're looking from point "start" towards point "s" and you can see the point "p", return true.
1573          * Otherwise return false.
1574          * @param {JXG.Coords} start The point you're standing on.
1575          * @param {JXG.Coords} p The point in which direction you're looking.
1576          * @param {JXG.Coords} s The point that should be visible.
1577          * @returns {Boolean} True, if from start the point p is in the same direction as s is, that means s-start = k*(p-start) with k>=0.
1578          */
1579         isSameDirection: function (start, p, s) {
1580             var dx,
1581                 dy,
1582                 sx,
1583                 sy,
1584                 r = false;
1585 
1586             dx = p.usrCoords[1] - start.usrCoords[1];
1587             dy = p.usrCoords[2] - start.usrCoords[2];
1588 
1589             sx = s.usrCoords[1] - start.usrCoords[1];
1590             sy = s.usrCoords[2] - start.usrCoords[2];
1591 
1592             if (Math.abs(dx) < Mat.eps) {
1593                 dx = 0;
1594             }
1595 
1596             if (Math.abs(dy) < Mat.eps) {
1597                 dy = 0;
1598             }
1599 
1600             if (Math.abs(sx) < Mat.eps) {
1601                 sx = 0;
1602             }
1603 
1604             if (Math.abs(sy) < Mat.eps) {
1605                 sy = 0;
1606             }
1607 
1608             if (dx >= 0 && sx >= 0) {
1609                 r = (dy >= 0 && sy >= 0) || (dy <= 0 && sy <= 0);
1610             } else if (dx <= 0 && sx <= 0) {
1611                 r = (dy >= 0 && sy >= 0) || (dy <= 0 && sy <= 0);
1612             }
1613 
1614             return r;
1615         },
1616 
1617         /**
1618          * Determinant of three points in the Euclidean plane.
1619          * Zero, if the points are collinear. Used to determine of a point q is left or
1620          * right to a segment defined by points p1 and p2.
1621          * <p>
1622          * Non-homogeneous version.
1623          *
1624          * @param  {Array|JXG.Point} p1 First point or its coordinates of the segment. Point object or array of length 3. First (homogeneous) coordinate is equal to 1.
1625          * @param  {Array|JXG.Point} p2 Second point or its coordinates of the segment. Point object or array of length 3. First (homogeneous) coordinate is equal to 1.
1626          * @param  {Array|JXG.Point} q Point or its coordinates. Point object or array of length 3. First (homogeneous) coordinate is equal to 1.
1627          * @return {Number} Signed area of the triangle formed by these three points.
1628          *
1629          * @see JXG.Math.Geometry.windingNumber
1630          */
1631         det3p: function (p1, p2, q) {
1632             var pp1, pp2, qq;
1633 
1634             if (Type.isPoint(p1)) {
1635                 pp1 = p1.Coords(true);
1636             } else {
1637                 pp1 = p1;
1638             }
1639             if (Type.isPoint(p2)) {
1640                 pp2 = p2.Coords(true);
1641             } else {
1642                 pp2 = p2;
1643             }
1644             if (Type.isPoint(q)) {
1645                 qq = q.Coords(true);
1646             } else {
1647                 qq = q;
1648             }
1649 
1650             return (pp1[1] - qq[1]) * (pp2[2] - qq[2]) - (pp2[1] - qq[1]) * (pp1[2] - qq[2]);
1651         },
1652 
1653         /**
1654          * Winding number of a point in respect to a polygon path.
1655          *
1656          * The point is regarded outside if the winding number is zero,
1657          * inside otherwise. The algorithm tries to find degenerate cases, i.e.
1658          * if the point is on the path. This is regarded as "outside".
1659          * If the point is a vertex of the path, it is regarded as "inside".
1660          *
1661          * Implementation of algorithm 7 from "The point in polygon problem for
1662          * arbitrary polygons" by Kai Hormann and Alexander Agathos, Computational Geometry,
1663          * Volume 20, Issue 3, November 2001, Pages 131-144.
1664          *
1665          * @param  {Array} usrCoords Homogenous coordinates of the point
1666          * @param  {Array} path      Array of points / coords determining a path, i.e. the vertices of the polygon / path. The array elements
1667          * do not have to be full points, but have to have a subobject "coords" or should be of type JXG.Coords.
1668          * @param  {Boolean} [doNotClosePath=false] If true the last point of the path is not connected to the first point.
1669          * This is necessary if the path consists of two or more closed subpaths, e.g. if the figure has a hole.
1670          *
1671          * @return {Number}          Winding number of the point. The point is
1672          *                           regarded outside if the winding number is zero,
1673          *                           inside otherwise.
1674          */
1675         windingNumber: function (usrCoords, path, doNotClosePath) {
1676             var wn = 0,
1677                 le = path.length,
1678                 x = usrCoords[1],
1679                 y = usrCoords[2],
1680                 p0,
1681                 p1,
1682                 p2,
1683                 d,
1684                 sign,
1685                 i,
1686                 off = 0;
1687 
1688             if (le === 0) {
1689                 return 0;
1690             }
1691 
1692             doNotClosePath = doNotClosePath || false;
1693             if (doNotClosePath) {
1694                 off = 1;
1695             }
1696 
1697             // Infinite points are declared outside
1698             if (isNaN(x) || isNaN(y)) {
1699                 return 1;
1700             }
1701 
1702             if (Type.exists(path[0].coords)) {
1703                 p0 = path[0].coords;
1704                 p1 = path[le - 1].coords;
1705             } else {
1706                 p0 = path[0];
1707                 p1 = path[le - 1];
1708             }
1709             // Handle the case if the point is the first vertex of the path, i.e. inside.
1710             if (p0.usrCoords[1] === x && p0.usrCoords[2] === y) {
1711                 return 1;
1712             }
1713 
1714             for (i = 0; i < le - off; i++) {
1715                 // Consider the edge from p1 = path[i] to p2 = path[i+1]isClosedPath
1716                 if (Type.exists(path[i].coords)) {
1717                     p1 = path[i].coords.usrCoords;
1718                     p2 = path[(i + 1) % le].coords.usrCoords;
1719                 } else {
1720                     p1 = path[i].usrCoords;
1721                     p2 = path[(i + 1) % le].usrCoords;
1722                 }
1723 
1724                 // If one of the two points p1, p2 is undefined or infinite,
1725                 // move on.
1726                 if (
1727                     p1[0] === 0 ||
1728                     p2[0] === 0 ||
1729                     isNaN(p1[1]) ||
1730                     isNaN(p2[1]) ||
1731                     isNaN(p1[2]) ||
1732                     isNaN(p2[2])
1733                 ) {
1734                     continue;
1735                 }
1736 
1737                 if (p2[2] === y) {
1738                     if (p2[1] === x) {
1739                         return 1;
1740                     }
1741                     if (p1[2] === y && p2[1] > x === p1[1] < x) {
1742                         return 0;
1743                     }
1744                 }
1745 
1746                 if (p1[2] < y !== p2[2] < y) {
1747                     // Crossing
1748                     sign = 2 * (p2[2] > p1[2] ? 1 : 0) - 1;
1749                     if (p1[1] >= x) {
1750                         if (p2[1] > x) {
1751                             wn += sign;
1752                         } else {
1753                             d = this.det3p(p1, p2, usrCoords);
1754                             if (d === 0) {
1755                                 // Point is on line, i.e. outside
1756                                 return 0;
1757                             }
1758                             if (d > 0 + Mat.eps === p2[2] > p1[2]) {
1759                                 // Right crossing
1760                                 wn += sign;
1761                             }
1762                         }
1763                     } else {
1764                         if (p2[1] > x) {
1765                             d = this.det3p(p1, p2, usrCoords);
1766                             if (d > 0 + Mat.eps === p2[2] > p1[2]) {
1767                                 // Right crossing
1768                                 wn += sign;
1769                             }
1770                         }
1771                     }
1772                 }
1773             }
1774 
1775             return wn;
1776         },
1777 
1778         /**
1779          * Decides if a point (x,y) is inside of a path / polygon.
1780          * Does not work correct if the path has hole. In this case, windingNumber is the preferred method.
1781          * Implements W. Randolf Franklin's pnpoly method.
1782          *
1783          * See <a href="https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html">https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html</a>.
1784          *
1785          * @param {Number} x_in x-coordinate (screen or user coordinates)
1786          * @param {Number} y_in y-coordinate (screen or user coordinates)
1787          * @param  {Array} path  Array of points / coords determining a path, i.e. the vertices of the polygon / path. The array elements
1788          * do not have to be full points, but have to have a subobject "coords" or should be of type JXG.Coords.
1789          * @param {Number} [coord_type=JXG.COORDS_BY_SCREEN] Type of coordinates used here.
1790          *   Possible values are <b>JXG.COORDS_BY_USER</b> and <b>JXG.COORDS_BY_SCREEN</b>.
1791          *   Default value is JXG.COORDS_BY_SCREEN.
1792          * @param {JXG.Board} board Board object
1793          *
1794          * @returns {Boolean} if (x_in, y_in) is inside of the polygon.
1795          * @see JXG.Polygon#hasPoint
1796          * @see JXG.Polygon#pnpoly
1797          * @see JXG.Math.Geometry.windingNumber
1798          *
1799          * @example
1800          * var pol = board.create('polygon', [[-1,2], [2,2], [-1,4]]);
1801          * var p = board.create('point', [4, 3]);
1802          * var txt = board.create('text', [-1, 0.5, function() {
1803          *   return 'Point A is inside of the polygon = ' +
1804          *     JXG.Math.Geometry.pnpoly(p.X(), p.Y(), pol.vertices, JXG.COORDS_BY_USER, board);
1805          * }]);
1806          *
1807          * </pre><div id="JXG4656ed42-f965-4e35-bb66-c334a4529683" class="jxgbox" style="width: 300px; height: 300px;"></div>
1808          * <script type="text/javascript">
1809          *     (function() {
1810          *         var board = JXG.JSXGraph.initBoard('JXG4656ed42-f965-4e35-bb66-c334a4529683',
1811          *             {boundingbox: [-2, 5, 5,-2], axis: true, showcopyright: false, shownavigation: false});
1812          *     var pol = board.create('polygon', [[-1,2], [2,2], [-1,4]]);
1813          *     var p = board.create('point', [4, 3]);
1814          *     var txt = board.create('text', [-1, 0.5, function() {
1815          *     		return 'Point A is inside of the polygon = ' + JXG.Math.Geometry.pnpoly(p.X(), p.Y(), pol.vertices, JXG.COORDS_BY_USER, board);
1816          *     }]);
1817          *
1818          *     })();
1819          *
1820          * </script><pre>
1821          *
1822          */
1823         pnpoly: function (x_in, y_in, path, coord_type, board) {
1824             var i, j, vi, vj, len,
1825                 x, y, crds,
1826                 v = path,
1827                 isIn = false;
1828 
1829             if (coord_type === Const.COORDS_BY_USER) {
1830                 crds = new Coords(Const.COORDS_BY_USER, [x_in, y_in], board);
1831                 x = crds.scrCoords[1];
1832                 y = crds.scrCoords[2];
1833             } else {
1834                 x = x_in;
1835                 y = y_in;
1836             }
1837 
1838             len = path.length;
1839             for (i = 0, j = len - 2; i < len - 1; j = i++) {
1840                 vi = Type.exists(v[i].coords) ? v[i].coords : v[i];
1841                 vj = Type.exists(v[j].coords) ? v[j].coords : v[j];
1842 
1843                 if (
1844                     vi.scrCoords[2] > y !== vj.scrCoords[2] > y &&
1845                     x <
1846                     ((vj.scrCoords[1] - vi.scrCoords[1]) * (y - vi.scrCoords[2])) /
1847                     (vj.scrCoords[2] - vi.scrCoords[2]) +
1848                     vi.scrCoords[1]
1849                 ) {
1850                     isIn = !isIn;
1851                 }
1852             }
1853 
1854             return isIn;
1855         },
1856 
1857         /****************************************/
1858         /****          INTERSECTIONS         ****/
1859         /****************************************/
1860 
1861         /**
1862          * Generate the function which computes the coordinates of the intersection point.
1863          * Primarily used in {@link JXG.Point.createIntersectionPoint}.
1864          * @param {JXG.Board} board object
1865          * @param {JXG.Line,JXG.Circle_JXG.Line,JXG.Circle_Number|Function} el1,el2,i The result will be a intersection point on el1 and el2.
1866          * i determines the intersection point if two points are available: <ul>
1867          *   <li>i==0: use the positive square root,</li>
1868          *   <li>i==1: use the negative square root.</li></ul>
1869          * @param {Boolean} alwaysintersect. Flag that determines if segments and arc can have an outer intersection point
1870          * on their defining line or circle.
1871          * @returns {Function} Function returning a {@link JXG.Coords} object that determines
1872          * the intersection point.
1873          *
1874          * @see JXG.Point.createIntersectionPoint
1875          */
1876         intersectionFunction: function (board, el1, el2, i, j, alwaysintersect) {
1877             var func,
1878                 that = this,
1879                 el1_isArcType = false,
1880                 el2_isArcType = false;
1881 
1882             el1_isArcType =
1883                 el1.elementClass === Const.OBJECT_CLASS_CURVE &&
1884                     (el1.type === Const.OBJECT_TYPE_ARC || el1.type === Const.OBJECT_TYPE_SECTOR)
1885                     ? true
1886                     : false;
1887             el2_isArcType =
1888                 el2.elementClass === Const.OBJECT_CLASS_CURVE &&
1889                     (el2.type === Const.OBJECT_TYPE_ARC || el2.type === Const.OBJECT_TYPE_SECTOR)
1890                     ? true
1891                     : false;
1892 
1893             if (
1894                 (el1.elementClass === Const.OBJECT_CLASS_CURVE ||
1895                     el2.elementClass === Const.OBJECT_CLASS_CURVE) &&
1896                 (el1.elementClass === Const.OBJECT_CLASS_CURVE ||
1897                     el1.elementClass === Const.OBJECT_CLASS_CIRCLE) &&
1898                 (el2.elementClass === Const.OBJECT_CLASS_CURVE ||
1899                     el2.elementClass === Const.OBJECT_CLASS_CIRCLE) /*&&
1900                 !(el1_isArcType && el2_isArcType)*/
1901             ) {
1902                 // curve - curve
1903                 // with the exception that both elements are arc types
1904                 /** @ignore */
1905                 func = function () {
1906                     return that.meetCurveCurve(el1, el2, i, j, el1.board);
1907                 };
1908             } else if (
1909                 (el1.elementClass === Const.OBJECT_CLASS_CURVE &&
1910                     !el1_isArcType &&
1911                     el2.elementClass === Const.OBJECT_CLASS_LINE) ||
1912                 (el2.elementClass === Const.OBJECT_CLASS_CURVE &&
1913                     !el2_isArcType &&
1914                     el1.elementClass === Const.OBJECT_CLASS_LINE)
1915             ) {
1916                 // curve - line (this includes intersections between conic sections and lines)
1917                 // with the exception that the curve is of arc type
1918                 /** @ignore */
1919                 func = function () {
1920                     return that.meetCurveLine(el1, el2, i, el1.board, Type.evaluate(alwaysintersect));
1921                 };
1922             } else if (
1923                 el1.type === Const.OBJECT_TYPE_POLYGON ||
1924                 el2.type === Const.OBJECT_TYPE_POLYGON
1925             ) {
1926                 // polygon - other
1927                 // Uses the Greiner-Hormann clipping algorithm
1928                 // Not implemented: polygon - point
1929 
1930                 if (el1.elementClass === Const.OBJECT_CLASS_LINE) {
1931                     // line - path
1932                     /** @ignore */
1933                     func = function () {
1934                         var first1 = el1.evalVisProp('straightfirst'),
1935                             last1 = el1.evalVisProp('straightlast'),
1936                             first2 = el2.evalVisProp('straightfirst'),
1937                             last2 = el2.evalVisProp('straightlast'),
1938                             a_not;
1939 
1940                         a_not = (!Type.evaluate(alwaysintersect) && (!first1 || !last1 || !first2 || !last2));
1941                         return that.meetPolygonLine(el2, el1, i, el1.board, a_not);
1942                     };
1943                 } else if (el2.elementClass === Const.OBJECT_CLASS_LINE) {
1944                     // path - line
1945                     /** @ignore */
1946                     func = function () {
1947                         var first1 = el1.evalVisProp('straightfirst'),
1948                             last1 = el1.evalVisProp('straightlast'),
1949                             first2 = el2.evalVisProp('straightfirst'),
1950                             last2 = el2.evalVisProp('straightlast'),
1951                             a_not;
1952 
1953                         a_not = (!Type.evaluate(alwaysintersect) && (!first1 || !last1 || !first2 || !last2));
1954                         return that.meetPolygonLine(el1, el2, i, el1.board, a_not);
1955                     };
1956                 } else {
1957                     // path - path
1958                     /** @ignore */
1959                     func = function () {
1960                         return that.meetPathPath(el1, el2, i, el1.board);
1961                     };
1962                 }
1963             } else if (
1964                 el1.elementClass === Const.OBJECT_CLASS_LINE &&
1965                 el2.elementClass === Const.OBJECT_CLASS_LINE
1966             ) {
1967                 // line - line, lines may also be segments.
1968                 /** @ignore */
1969                 func = function () {
1970                     var res,
1971                         c,
1972                         first1 = el1.evalVisProp('straightfirst'),
1973                         last1 = el1.evalVisProp('straightlast'),
1974                         first2 = el2.evalVisProp('straightfirst'),
1975                         last2 = el2.evalVisProp('straightlast');
1976 
1977                     /**
1978                      * If one of the lines is a segment or ray and
1979                      * the intersection point should disappear if outside
1980                      * of the segment or ray we call
1981                      * meetSegmentSegment
1982                      */
1983                     if (
1984                         !Type.evaluate(alwaysintersect) &&
1985                         (!first1 || !last1 || !first2 || !last2)
1986                     ) {
1987                         res = that.meetSegmentSegment(
1988                             el1.point1.coords.usrCoords,
1989                             el1.point2.coords.usrCoords,
1990                             el2.point1.coords.usrCoords,
1991                             el2.point2.coords.usrCoords
1992                         );
1993 
1994                         if (
1995                             (!first1 && res[1] < 0) ||
1996                             (!last1 && res[1] > 1) ||
1997                             (!first2 && res[2] < 0) ||
1998                             (!last2 && res[2] > 1)
1999                         ) {
2000                             // Non-existent
2001                             c = [0, NaN, NaN];
2002                         } else {
2003                             c = res[0];
2004                         }
2005 
2006                         return new Coords(Const.COORDS_BY_USER, c, el1.board);
2007                     }
2008 
2009                     return that.meet(el1.stdform, el2.stdform, i, el1.board);
2010                 };
2011             } else {
2012                 // All other combinations of circles and lines,
2013                 // Arc types are treated as circles.
2014                 /** @ignore */
2015                 func = function () {
2016                     var res = that.meet(el1.stdform, el2.stdform, i, el1.board),
2017                         has = true,
2018                         first,
2019                         last,
2020                         r;
2021 
2022                     if (Type.evaluate(alwaysintersect)) {
2023                         return res;
2024                     }
2025                     if (el1.elementClass === Const.OBJECT_CLASS_LINE) {
2026                         first = el1.evalVisProp('straightfirst');
2027                         last = el1.evalVisProp('straightlast');
2028                         if (!first || !last) {
2029                             r = that.affineRatio(el1.point1.coords, el1.point2.coords, res);
2030                             if ((!last && r > 1 + Mat.eps) || (!first && r < 0 - Mat.eps)) {
2031                                 return new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board);
2032                             }
2033                         }
2034                     }
2035                     if (el2.elementClass === Const.OBJECT_CLASS_LINE) {
2036                         first = el2.evalVisProp('straightfirst');
2037                         last = el2.evalVisProp('straightlast');
2038                         if (!first || !last) {
2039                             r = that.affineRatio(el2.point1.coords, el2.point2.coords, res);
2040                             if ((!last && r > 1 + Mat.eps) || (!first && r < 0 - Mat.eps)) {
2041                                 return new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board);
2042                             }
2043                         }
2044                     }
2045                     if (el1_isArcType) {
2046                         has = that.coordsOnArc(el1, res);
2047                         if (has && el2_isArcType) {
2048                             has = that.coordsOnArc(el2, res);
2049                         }
2050                         if (!has) {
2051                             return new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board);
2052                         }
2053                     }
2054                     return res;
2055                 };
2056             }
2057 
2058             return func;
2059         },
2060 
2061         otherIntersectionFunction: function (input, others, alwaysintersect, precision) {
2062             var func, board,
2063                 el1, el2,
2064                 that = this;
2065 
2066             el1 = input[0];
2067             el2 = input[1];
2068             board = el1.board;
2069             /** @ignore */
2070             func = function () {
2071                 var i, k, c, d,
2072                     isClose,
2073                     le = others.length,
2074                     eps = Type.evaluate(precision);
2075 
2076                 for (i = le; i >= 0; i--) {
2077                     if (el1.elementClass === Const.OBJECT_CLASS_CIRCLE &&
2078                         [Const.OBJECT_CLASS_CIRCLE, Const.OBJECT_CLASS_LINE].indexOf(el2.elementClass) >= 0) {
2079                         // circle, circle|line
2080                         c = that.meet(el1.stdform, el2.stdform, i, board);
2081                     } else if (el1.elementClass === Const.OBJECT_CLASS_CURVE &&
2082                         [Const.OBJECT_CLASS_CURVE, Const.OBJECT_CLASS_CIRCLE].indexOf(el2.elementClass) >= 0) {
2083                         // curve, circle|curve
2084                         c = that.meetCurveCurve(el1, el2, i, 0, board);
2085                     } else if (el1.elementClass === Const.OBJECT_CLASS_CURVE && el2.elementClass === Const.OBJECT_CLASS_LINE) {
2086                         // curve, line
2087                         if (Type.exists(el1.dataX)) {
2088                             c = JXG.Math.Geometry.meetCurveLine(el1, el2, i, el1.board, Type.evaluate(alwaysintersect));
2089                         } else {
2090                             c = JXG.Math.Geometry.meetCurveLineContinuous(el1, el2, i, el1.board);
2091                         }
2092                     }
2093 
2094                     if (c === undefined) {
2095                         // Intersection point does not exist
2096                         continue;
2097                     }
2098 
2099                     // If the intersection is close to one of the points in other
2100                     // we have to search for another intersection point.
2101                     isClose = false;
2102                     for (k = 0; !isClose && k < le; k++) {
2103                         if (Type.exists(c) && Type.exists(c.distance)) {
2104                             d = c.distance(JXG.COORDS_BY_USER, others[k].coords);
2105                             if (d < eps) {
2106                                 isClose = true;
2107                             }
2108                         }
2109                     }
2110                     if (!isClose) {
2111                         // We are done, the intersection is away from any other
2112                         // intersection point.
2113                         return c;
2114                     }
2115                 }
2116                 // Otherwise we return the last intersection point
2117                 return c;
2118             };
2119             return func;
2120         },
2121 
2122         /**
2123          * Returns true if the coordinates are on the arc element,
2124          * false otherwise. Usually, coords is an intersection
2125          * on the circle line. Now it is decided if coords are on the
2126          * circle restricted to the arc line.
2127          * @param  {Arc} arc arc or sector element
2128          * @param  {JXG.Coords} coords Coords object of an intersection
2129          * @returns {Boolean}
2130          * @private
2131          */
2132         coordsOnArc: function (arc, coords) {
2133             var angle = this.rad(arc.radiuspoint, arc.center, coords.usrCoords.slice(1)),
2134                 alpha = 0.0,
2135                 beta = this.rad(arc.radiuspoint, arc.center, arc.anglepoint),
2136                 ev_s = arc.evalVisProp('selection');
2137 
2138             if ((ev_s === "minor" && beta > Math.PI) || (ev_s === "major" && beta < Math.PI)) {
2139                 alpha = beta;
2140                 beta = 2 * Math.PI;
2141             }
2142             if (angle < alpha || angle > beta) {
2143                 return false;
2144             }
2145             return true;
2146         },
2147 
2148         /**
2149          * Computes the intersection of a pair of lines, circles or both.
2150          * It uses the internal data array stdform of these elements.
2151          * @param {Array} el1 stdform of the first element (line or circle)
2152          * @param {Array} el2 stdform of the second element (line or circle)
2153          * @param {Number|Function} i Index of the intersection point that should be returned.
2154          * @param board Reference to the board.
2155          * @returns {JXG.Coords} Coordinates of one of the possible two or more intersection points.
2156          * Which point will be returned is determined by i.
2157          */
2158         meet: function (el1, el2, i, board) {
2159             var result,
2160                 eps = Mat.eps;
2161 
2162             if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) < eps) {
2163                 // line line
2164                 result = this.meetLineLine(el1, el2, i, board);
2165             } else if (Math.abs(el1[3]) >= eps && Math.abs(el2[3]) < eps) {
2166                 // circle line
2167                 result = this.meetLineCircle(el2, el1, i, board);
2168             } else if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) >= eps) {
2169                 // line circle
2170                 result = this.meetLineCircle(el1, el2, i, board);
2171             } else {
2172                 // circle circle
2173                 result = this.meetCircleCircle(el1, el2, i, board);
2174             }
2175 
2176             return result;
2177         },
2178 
2179         /**
2180          * Intersection of the line with the board
2181          * @param  {Array}     line   stdform of the line in screen coordinates
2182          * @param  {JXG.Board} board  reference to a board.
2183          * @param  {Number}    margin optional margin, to avoid the display of the small sides of lines.
2184          * @returns {Array}            [intersection coords 1, intersection coords 2]
2185          */
2186         meetLineBoard: function (line, board, margin) {
2187             // Intersect the line with the four borders of the board.
2188             var s = [],
2189                 intersect1,
2190                 intersect2,
2191                 i, j;
2192 
2193             if (!Type.exists(margin)) {
2194                 margin = 0;
2195             }
2196 
2197             // top
2198             s[0] = Mat.crossProduct(line, [margin, 0, 1]);
2199             // left
2200             s[1] = Mat.crossProduct(line, [margin, 1, 0]);
2201             // bottom
2202             s[2] = Mat.crossProduct(line, [-margin - board.canvasHeight, 0, 1]);
2203             // right
2204             s[3] = Mat.crossProduct(line, [-margin - board.canvasWidth, 1, 0]);
2205 
2206             // Normalize the intersections
2207             for (i = 0; i < 4; i++) {
2208                 if (Math.abs(s[i][0]) > Mat.eps) {
2209                     for (j = 2; j > 0; j--) {
2210                         s[i][j] /= s[i][0];
2211                     }
2212                     s[i][0] = 1.0;
2213                 }
2214             }
2215 
2216             // line is parallel to "left", take "top" and "bottom"
2217             if (Math.abs(s[1][0]) < Mat.eps) {
2218                 intersect1 = s[0]; // top
2219                 intersect2 = s[2]; // bottom
2220                 // line is parallel to "top", take "left" and "right"
2221             } else if (Math.abs(s[0][0]) < Mat.eps) {
2222                 intersect1 = s[1]; // left
2223                 intersect2 = s[3]; // right
2224                 // left intersection out of board (above)
2225             } else if (s[1][2] < 0) {
2226                 intersect1 = s[0]; // top
2227 
2228                 // right intersection out of board (below)
2229                 if (s[3][2] > board.canvasHeight) {
2230                     intersect2 = s[2]; // bottom
2231                 } else {
2232                     intersect2 = s[3]; // right
2233                 }
2234                 // left intersection out of board (below)
2235             } else if (s[1][2] > board.canvasHeight) {
2236                 intersect1 = s[2]; // bottom
2237 
2238                 // right intersection out of board (above)
2239                 if (s[3][2] < 0) {
2240                     intersect2 = s[0]; // top
2241                 } else {
2242                     intersect2 = s[3]; // right
2243                 }
2244             } else {
2245                 intersect1 = s[1]; // left
2246 
2247                 // right intersection out of board (above)
2248                 if (s[3][2] < 0) {
2249                     intersect2 = s[0]; // top
2250                     // right intersection out of board (below)
2251                 } else if (s[3][2] > board.canvasHeight) {
2252                     intersect2 = s[2]; // bottom
2253                 } else {
2254                     intersect2 = s[3]; // right
2255                 }
2256             }
2257 
2258             return [
2259                 new Coords(Const.COORDS_BY_SCREEN, intersect1.slice(1), board),
2260                 new Coords(Const.COORDS_BY_SCREEN, intersect2.slice(1), board)
2261             ];
2262         },
2263 
2264         /**
2265          * Intersection of two lines.
2266          * @param {Array} l1 stdform of the first line
2267          * @param {Array} l2 stdform of the second line
2268          * @param {number} i unused
2269          * @param {JXG.Board} board Reference to the board.
2270          * @returns {JXG.Coords} Coordinates of the intersection point.
2271          */
2272         meetLineLine: function (l1, l2, i, board) {
2273             var s = isNaN(l1[5] + l2[5]) ? [0, 0, 0] : Mat.crossProduct(l1, l2);
2274 
2275             // Make intersection of parallel lines more robust:
2276             if (Math.abs(s[0]) < 1.0e-14) {
2277                 s[0] = 0;
2278             }
2279             return new Coords(Const.COORDS_BY_USER, s, board);
2280         },
2281 
2282         /**
2283          * Intersection of line and circle.
2284          * @param {Array} lin stdform of the line
2285          * @param {Array} circ stdform of the circle
2286          * @param {number|function} i number of the returned intersection point.
2287          *   i==0: use the positive square root,
2288          *   i==1: use the negative square root.
2289          * @param {JXG.Board} board Reference to a board.
2290          * @returns {JXG.Coords} Coordinates of the intersection point
2291          */
2292         meetLineCircle: function (lin, circ, i, board) {
2293             var a, b, c, d, n, A, B, C, k, t;
2294 
2295             // Radius is zero, return center of circle
2296             if (circ[4] < Mat.eps) {
2297                 if (Math.abs(Mat.innerProduct([1, circ[6], circ[7]], lin, 3)) < Mat.eps) {
2298                     return new Coords(Const.COORDS_BY_USER, circ.slice(6, 8), board);
2299                 }
2300 
2301                 return new Coords(Const.COORDS_BY_USER, [NaN, NaN], board);
2302             }
2303             c = circ[0];
2304             b = circ.slice(1, 3);
2305             a = circ[3];
2306             d = lin[0];
2307             n = lin.slice(1, 3);
2308 
2309             // Line is assumed to be normalized. Therefore, nn==1 and we can skip some operations:
2310             /*
2311              var nn = n[0]*n[0]+n[1]*n[1];
2312              A = a*nn;
2313              B = (b[0]*n[1]-b[1]*n[0])*nn;
2314              C = a*d*d - (b[0]*n[0]+b[1]*n[1])*d + c*nn;
2315              */
2316             A = a;
2317             B = b[0] * n[1] - b[1] * n[0];
2318             C = a * d * d - (b[0] * n[0] + b[1] * n[1]) * d + c;
2319 
2320             k = B * B - 4 * A * C;
2321             if (k > -Mat.eps * Mat.eps) {
2322                 k = Math.sqrt(Math.abs(k));
2323                 t = [(-B + k) / (2 * A), (-B - k) / (2 * A)];
2324 
2325                 return Type.evaluate(i) === 0
2326                     ? new Coords(
2327                         Const.COORDS_BY_USER,
2328                         [-t[0] * -n[1] - d * n[0], -t[0] * n[0] - d * n[1]],
2329                         board
2330                     )
2331                     : new Coords(
2332                         Const.COORDS_BY_USER,
2333                         [-t[1] * -n[1] - d * n[0], -t[1] * n[0] - d * n[1]],
2334                         board
2335                     );
2336             }
2337 
2338             return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
2339         },
2340 
2341         /**
2342          * Intersection of two circles.
2343          * @param {Array} circ1 stdform of the first circle
2344          * @param {Array} circ2 stdform of the second circle
2345          * @param {number|function} i number of the returned intersection point.
2346          *   i==0: use the positive square root,
2347          *   i==1: use the negative square root.
2348          * @param {JXG.Board} board Reference to the board.
2349          * @returns {JXG.Coords} Coordinates of the intersection point
2350          */
2351         meetCircleCircle: function (circ1, circ2, i, board) {
2352             var radicalAxis;
2353 
2354             // Radius is zero, return center of circle, if on other circle
2355             if (circ1[4] < Mat.eps) {
2356                 if (
2357                     Math.abs(this.distance(circ1.slice(6, 2), circ2.slice(6, 8)) - circ2[4]) <
2358                     Mat.eps
2359                 ) {
2360                     return new Coords(Const.COORDS_BY_USER, circ1.slice(6, 8), board);
2361                 }
2362 
2363                 return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
2364             }
2365 
2366             // Radius is zero, return center of circle, if on other circle
2367             if (circ2[4] < Mat.eps) {
2368                 if (
2369                     Math.abs(this.distance(circ2.slice(6, 2), circ1.slice(6, 8)) - circ1[4]) <
2370                     Mat.eps
2371                 ) {
2372                     return new Coords(Const.COORDS_BY_USER, circ2.slice(6, 8), board);
2373                 }
2374 
2375                 return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
2376             }
2377 
2378             radicalAxis = [
2379                 circ2[3] * circ1[0] - circ1[3] * circ2[0],
2380                 circ2[3] * circ1[1] - circ1[3] * circ2[1],
2381                 circ2[3] * circ1[2] - circ1[3] * circ2[2],
2382                 0,
2383                 1,
2384                 Infinity,
2385                 Infinity,
2386                 Infinity
2387             ];
2388             radicalAxis = Mat.normalize(radicalAxis);
2389 
2390             return this.meetLineCircle(radicalAxis, circ1, i, board);
2391         },
2392 
2393         /**
2394          * Segment-wise search for the nr-th intersection of two curves.
2395          * testSegment is always assumed to be true.
2396          *
2397          * @param {JXG.Curve} c1 Curve, Line or Circle
2398          * @param {JXG.Curve} c2 Curve, Line or Circle
2399          * @param {Number} nr the nr-th intersection point will be returned
2400          * @param {JXG.Board} [board=c1.board] Reference to a board object
2401          * @returns {JXG.Coords} intersection as Coords object
2402          *
2403          * @private
2404          * @see JXG.Math.Geometry.meetCurveCurve
2405          */
2406         meetCurveCurveDiscrete: function (c1, c2, nr, board) {
2407             var co,
2408                 i = Type.evaluate(nr);
2409 
2410             if (c1.bezierDegree === 3 || c2.bezierDegree === 3) {
2411                 co = this.meetBezierCurveRedBlueSegments(c1, c2, i);
2412             } else {
2413                 co = this.meetCurveRedBlueSegments(c1, c2, i);
2414             }
2415             return new Coords(Const.COORDS_BY_USER, co, board);
2416         },
2417 
2418         /**
2419          * Apply Newton-Raphson to search for an intersection of two curves
2420          * in a given range of the first curve.
2421          *
2422          * @param {JXG.Curve} c1 Curve, Line or Circle
2423          * @param {JXG.Curve} c2 Curve, Line or Circle
2424          * @param {Array} range Domain for the search of an intersection. The start value
2425          * for the search is chosen to be inside of that range.
2426          * @param {Boolean} testSegment If true require that t1 and t2 are inside of the allowed bounds.
2427          * @returns {Array} [[z, x, y], t1, t2, t, ||c1[t1]-c2[t2]||**2]. The last entry is set to
2428          * 10000 if the intersection is outside of the given domain (range) for the first curve.
2429          * @private
2430          * @see JXG.Math.Geometry._meetCurveCurveRecursive
2431          * @see JXG.Math.Geometry._meetCurveCurveIterative
2432          * @see JXG.Math.Numerics.generalizedDampedNewton
2433          * @see JXG.Math.Geometry.meetCurveCurveCobyla
2434          */
2435         meetCurveCurveNewton: function (c1, c2, range1, range2, testSegment) {
2436             var t1, t2,
2437                 co, r,
2438                 inphi = (Math.sqrt(5) - 1) * 0.5,
2439                 damp = 0.85, // (
2440                 eps3 = Mat.eps * Mat.eps * Mat.eps,
2441                 eps2 = Mat.eps * Mat.eps,
2442 
2443                 ma1 = c1.maxX(),
2444                 mi1 = c1.minX(),
2445                 ma2 = c2.maxX(),
2446                 mi2 = c2.minX(),
2447 
2448                 F = function(t, n) {
2449                     var f1 = c1.Ft(t[0]),
2450                         f2 = c2.Ft(t[1]),
2451                         e = f1[1] - f2[1],
2452                         f = f1[2] - f2[2];
2453 
2454                     return [e, f];
2455                 },
2456                 D = function(t, n) {
2457                     var h = Mat.eps,
2458                         h2 = 2 * h,
2459                         f1_1 = c1.Ft(t[0] - h),
2460                         f1_2 = c1.Ft(t[0] + h),
2461                         f2_1 = c2.Ft(t[1] - h),
2462                         f2_2 = c2.Ft(t[1] + h);
2463                     return [
2464                         [ (f1_2[1] - f1_1[1]) / h2,
2465                          -(f2_2[1] - f2_1[1]) / h2],
2466                         [ (f1_2[2] - f1_1[2]) / h2,
2467                          -(f2_2[2] - f2_1[2]) / h2]
2468                     ];
2469                 };
2470 
2471             t1 = range1[0] + (range1[1] - range1[0]) * (1 - inphi);
2472             t2 = range2[0] + (range2[1] - range2[0]) * (1 - inphi);
2473 
2474             // Use damped Newton
2475             // r = Numerics.generalizedDampedNewtonCurves(c1, c2, t1, t2, damp, eps3);
2476             r = Numerics.generalizedDampedNewton(F, D, 2, [t1, t2], damp, eps3, 40);
2477             // r: [t1, t2, F2]
2478 
2479             t1 = r[0][0];
2480             t2 = r[0][1];
2481             co = c1.Ft(t1);
2482 
2483             if (
2484                 t1 < range1[0] - Mat.eps || t1 > range1[1] + Mat.eps ||
2485                 t2 < range2[0] - Mat.eps || t2 > range2[1] + Mat.eps ||
2486                 (testSegment &&
2487                     (t1 < mi1 - eps2 || t1 > ma1 + eps2 ||
2488                      t2 < mi2 - eps2 || t2 > ma2 + eps2)
2489                 )
2490             ) {
2491                 // Damped-Newton found solution outside of range
2492                 return [co, t1, t2, 10000];
2493             }
2494 // console.log(t1, r[3])
2495 
2496             return [co, t1, t2, r[1]];
2497         },
2498 
2499         /**
2500          * Return a list of the (at most) first i intersection points of two curves.
2501          * Computed iteratively.
2502          *
2503          * @param {JXG.Curve} c1 Curve, Line or Circle
2504          * @param {JXG.Curve} c2 Curve, Line or Circle
2505          * @param {Number} low Lower bound of the search domain (between [0, 1])
2506          * @param {Number} up Upper bound of the search domain (between [0, 1])
2507          * @param {Number} i Return a list of the first i intersection points
2508          * @param {Boolean} testSegment If true require that t1 and t2 are inside of the allowed bounds.
2509          * @returns {Array} List of the first i intersection points, given by the parameter t.
2510          * @private
2511          * @see JXG.Math.Geometry.meetCurveCurveNewton
2512          * @see JXG.Math.Geometry.meetCurveCurve
2513          */
2514         _meetCurveCurveIterative: function(c1, c2, range1, range2, i, testSegment) {
2515             var ret,
2516                 t1,// t2,
2517                 low1, low2, up1, up2,
2518                 eps = Mat.eps * 100, // Minimum difference between zeros
2519                 j1, j2,
2520                 steps = 20,
2521                 d1, d2,
2522                 zeros = [];
2523 
2524             low1 = range1[0];
2525             up1 = range1[1];
2526             low2 = range2[0];
2527             up2 = range2[1];
2528             if (up1 < low1 || up2 < low2) {
2529                 return [];
2530             }
2531 
2532             // console.log('DO iterative', [low1, up1], [low2, up2])
2533 
2534             d1 = (up1 - low1) / steps;
2535             d2 = (up2 - low2) / steps;
2536             for (j1 = 0; j1 < steps; j1++) {
2537                 for (j2 = 0; j2 < steps; j2++) {
2538 
2539                     ret = this.meetCurveCurveNewton(c1, c2,
2540                         [low1 + j1 * d1, low1 + (j1 + 1) * d1],
2541                         [low2 + j2 * d2, low2 + (j2 + 1) * d2],
2542                         testSegment);
2543 
2544                     if (ret[3] < Mat.eps) {
2545                         t1 = ret[1];
2546                         // t2 = ret[2];
2547                         // console.log("\tFOUND", t1, t2, c1.Ft(t1)[2])
2548                         zeros = zeros.concat([t1]);
2549                         zeros = Type.toUniqueArrayFloat(zeros, eps);
2550                         // console.log(zeros, i)
2551                         if (zeros.length > i) {
2552                             return zeros;
2553                         }
2554                     }
2555                 }
2556             }
2557 
2558             return zeros;
2559         },
2560 
2561         /**
2562          * Compute an intersection of the curves c1 and c2.
2563          * We want to find values t1, t2 such that
2564          * c1(t1) = c2(t2), i.e. (c1_x(t1) - c2_x(t2), c1_y(t1) - c2_y(t2)) = (0, 0).
2565          *
2566          * Available methods:
2567          * <ul>
2568          *  <li> discrete, segment-wise intersections
2569          *  <li> generalized damped Newton-Raphson
2570          * </ul>
2571          *
2572          * Segment-wise intersection is more stable, but has problems with tangent points.
2573          * Damped Newton-Raphson converges very rapidly but sometimes behaves chaotic.
2574          *
2575          * @param {JXG.Curve} c1 Curve, Line or Circle
2576          * @param {JXG.Curve} c2 Curve, Line or Circle
2577          * @param {Number|Function} nr the nr-th intersection point will be returned. For backwards compatibility:
2578          * if method='newton' and nr is not an integer, {@link JXG.Math.Numerics.generalizedNewton} is called
2579          * directly with nr as start value (not recommended).
2580          * @param {Number} t2ini not longer used. Must be supplied and is ignored.
2581          * @param {JXG.Board} [board=c1.board] Reference to a board object.
2582          * @param {String} [method] Intersection method, possible values are 'newton' and 'segment'.
2583          * If both curves are given by functions (assumed to be continuous), 'newton' is the default, otherwise
2584          * 'segment' is the default.
2585          * @returns {JXG.Coords} intersection point
2586          *
2587          * @see JXG.Math.Geometry.meetCurveCurveDiscrete
2588          * @see JXG.Math.Geometry._meetCurveCurveRecursive
2589          * @see JXG.Math.Geometry.meetCurveCurveIterative
2590          */
2591         meetCurveCurve: function (c1, c2, nr, t2ini, board, method, testSegment) {
2592             var co,
2593                 zeros,
2594                 mi1, ma1, mi2, ma2,
2595                 i = Type.evaluate(nr);
2596 
2597             board = board || c1.board;
2598             if (method === 'segment' || Type.exists(c1.dataX) || Type.exists(c2.dataX)) {
2599                 // Discrete data points, i.e. x-coordinates of c1 or c2 are given in an array)
2600                 return this.meetCurveCurveDiscrete(c1, c2, i, board);
2601             }
2602 
2603             // Outdated:
2604             // Backwards compatibility if nr is not a positive integer then
2605             // generalizedNewton is still used.
2606             if (Type.exists(method) && method === 'newton' && i < 0 || parseInt(i) !== i) {
2607                 co = Numerics.generalizedNewton(c1, c2, i, t2ini);
2608                 return new Coords(Const.COORDS_BY_USER, co, board);
2609             }
2610 
2611             // Method 'newton'
2612             mi1 = c1.minX();
2613             ma1 = c1.maxX();
2614             mi2 = c2.minX();
2615             ma2 = c2.maxX();
2616 
2617             // console.time('curvecurve')
2618             zeros = this._meetCurveCurveIterative(c1, c2, [mi1, ma1], [mi2, ma2], i, testSegment);
2619             // console.timeEnd('curvecurve')
2620 
2621             if (zeros.length > i) {
2622                 co = c1.Ft(zeros[i]);
2623             } else {
2624                 return [0, NaN, NaN];
2625             }
2626 
2627             return new Coords(Const.COORDS_BY_USER, co, board);
2628         },
2629 
2630         /**
2631          * Intersection of curve with line,
2632          * Order of input does not matter for el1 and el2.
2633          * From version 0.99.7 on this method calls
2634          * {@link JXG.Math.Geometry.meetCurveLineDiscrete}.
2635          * If higher precision is needed, {@link JXG.Math.Geometry.meetCurveLineContinuous}
2636          * has to be used.
2637          *
2638          * @param {JXG.Curve|JXG.Line} el1 Curve or Line
2639          * @param {JXG.Curve|JXG.Line} el2 Curve or Line
2640          * @param {Number|Function} nr the nr-th intersection point will be returned.
2641          * @param {JXG.Board} [board=el1.board] Reference to a board object.
2642          * @param {Boolean} alwaysIntersect If false just the segment between the two defining points are tested for intersection
2643          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
2644          * the ideal point [0,1,0] is returned.
2645          */
2646         meetCurveLine: function (el1, el2, nr, board, alwaysIntersect) {
2647             var v = [0, NaN, NaN],
2648                 cu,
2649                 li;
2650 
2651             if (!Type.exists(board)) {
2652                 board = el1.board;
2653             }
2654 
2655             if (el1.elementClass === Const.OBJECT_CLASS_CURVE) {
2656                 cu = el1;
2657                 li = el2;
2658             } else {
2659                 cu = el2;
2660                 li = el1;
2661             }
2662 
2663             if (Type.exists(cu.dataX)) {
2664                 // We use the discrete version if
2665                 //   the curve is not a parametric curve, e.g. implicit plots
2666                 v = this.meetCurveLineDiscrete(cu, li, nr, board, !alwaysIntersect);
2667             } else {
2668                 v = this.meetCurveCurve(cu, li, nr, 0, board, 'newton', !alwaysIntersect);
2669             }
2670 
2671             return v;
2672         },
2673 
2674         /**
2675          * Intersection of line and curve, continuous case.
2676          * Finds the nr-th intersection point
2677          * Uses {@link JXG.Math.Geometry.meetCurveLineDiscrete} as a first approximation.
2678          * A more exact solution is then found with {@link JXG.Math.Numerics.root}.
2679          *
2680          * @param {JXG.Curve} cu Curve
2681          * @param {JXG.Line} li Line
2682          * @param {NumberFunction} nr Will return the nr-th intersection point.
2683          * @param {JXG.Board} board
2684          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the
2685          * line defined by the segment
2686          * @returns {JXG.Coords} Coords object containing the intersection.
2687          */
2688         meetCurveLineContinuous: function (cu, li, nr, board, testSegment) {
2689             var func0, func1,
2690                 t, v, x, y, z,
2691                 eps = Mat.eps,
2692                 epsLow = Mat.eps,
2693                 steps,
2694                 delta,
2695                 tnew, tmin, fmin,
2696                 i, ft;
2697 
2698             v = this.meetCurveLineDiscrete(cu, li, nr, board, testSegment);
2699             x = v.usrCoords[1];
2700             y = v.usrCoords[2];
2701 
2702             func0 = function (t) {
2703                 var c1, c2;
2704 
2705                 if (t > cu.maxX() || t < cu.minX()) {
2706                     return Infinity;
2707                 }
2708                 c1 = cu.X(t) - x;
2709                 c2 = cu.Y(t) - y;
2710                 return c1 * c1 + c2 * c2;
2711                 // return c1 * (cu.X(t + h) - cu.X(t - h)) + c2 * (cu.Y(t + h) - cu.Y(t - h)) / h;
2712             };
2713 
2714             func1 = function (t) {
2715                 var v = li.stdform[0] + li.stdform[1] * cu.X(t) + li.stdform[2] * cu.Y(t);
2716                 return v * v;
2717             };
2718 
2719             // Find t
2720             steps = 50;
2721             delta = (cu.maxX() - cu.minX()) / steps;
2722             tnew = cu.minX();
2723             fmin = 0.0001; //eps;
2724             tmin = NaN;
2725             for (i = 0; i < steps; i++) {
2726                 t = Numerics.root(func0, [
2727                     Math.max(tnew, cu.minX()),
2728                     Math.min(tnew + delta, cu.maxX())
2729                 ]);
2730                 ft = Math.abs(func0(t));
2731                 if (ft <= fmin) {
2732                     fmin = ft;
2733                     tmin = t;
2734                     if (fmin < eps) {
2735                         break;
2736                     }
2737                 }
2738 
2739                 tnew += delta;
2740             }
2741             t = tmin;
2742             // Compute "exact" t
2743             t = Numerics.root(func1, [
2744                 Math.max(t - delta, cu.minX()),
2745                 Math.min(t + delta, cu.maxX())
2746             ]);
2747 
2748             ft = func1(t);
2749             // Is the point on the line?
2750             if (isNaN(ft) || Math.abs(ft) > epsLow) {
2751                 z = 0.0; //NaN;
2752             } else {
2753                 z = 1.0;
2754             }
2755 
2756             return new Coords(Const.COORDS_BY_USER, [z, cu.X(t), cu.Y(t)], board);
2757         },
2758 
2759         /**
2760          * Intersection of line and curve, discrete case.
2761          * Segments are treated as lines.
2762          * Finding the nr-th intersection point should work for all nr.
2763          * @param {JXG.Curve} cu
2764          * @param {JXG.Line} li
2765          * @param {Number|Function} nr
2766          * @param {JXG.Board} board
2767          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the
2768          * line defined by the segment
2769          *
2770          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
2771          * the ideal point [0,1,0] is returned.
2772          */
2773         meetCurveLineDiscrete: function (cu, li, nr, board, testSegment) {
2774             var i, j,
2775                 n = Type.evaluate(nr),
2776                 p1, p2,
2777                 p, q,
2778                 lip1 = li.point1.coords.usrCoords,
2779                 lip2 = li.point2.coords.usrCoords,
2780                 d, res,
2781                 cnt = 0,
2782                 len = cu.numberPoints,
2783                 ev_sf = li.evalVisProp('straightfirst'),
2784                 ev_sl = li.evalVisProp('straightlast');
2785 
2786             // In case, no intersection will be found we will take this
2787             q = new Coords(Const.COORDS_BY_USER, [0, NaN, NaN], board);
2788 
2789             if (lip1[0] === 0.0) {
2790                 lip1 = [1, lip2[1] + li.stdform[2], lip2[2] - li.stdform[1]];
2791             } else if (lip2[0] === 0.0) {
2792                 lip2 = [1, lip1[1] + li.stdform[2], lip1[2] - li.stdform[1]];
2793             }
2794 
2795             p2 = cu.points[0].usrCoords;
2796             for (i = 1; i < len; i += cu.bezierDegree) {
2797                 p1 = p2.slice(0);
2798                 p2 = cu.points[i].usrCoords;
2799                 d = this.distance(p1, p2);
2800 
2801                 // The defining points are not identical
2802                 if (d > Mat.eps) {
2803                     if (cu.bezierDegree === 3) {
2804                         res = this.meetBeziersegmentBeziersegment(
2805                             [
2806                                 cu.points[i - 1].usrCoords.slice(1),
2807                                 cu.points[i].usrCoords.slice(1),
2808                                 cu.points[i + 1].usrCoords.slice(1),
2809                                 cu.points[i + 2].usrCoords.slice(1)
2810                             ],
2811                             [lip1.slice(1), lip2.slice(1)],
2812                             testSegment
2813                         );
2814                     } else {
2815                         res = [this.meetSegmentSegment(p1, p2, lip1, lip2)];
2816                     }
2817 
2818                     for (j = 0; j < res.length; j++) {
2819                         p = res[j];
2820                         if (0 <= p[1] && p[1] <= 1) {
2821                             if (cnt === n) {
2822                                 /**
2823                                  * If the intersection point is not part of the segment,
2824                                  * this intersection point is set to non-existent.
2825                                  * This prevents jumping behavior of the intersection points.
2826                                  * But it may be discussed if it is the desired behavior.
2827                                  */
2828                                 if (
2829                                     testSegment &&
2830                                     ((!ev_sf && p[2] < 0) || (!ev_sl && p[2] > 1))
2831                                 ) {
2832                                     return q; // break;
2833                                 }
2834 
2835                                 q = new Coords(Const.COORDS_BY_USER, p[0], board);
2836                                 return q; // break;
2837                             }
2838                             cnt += 1;
2839                         }
2840                     }
2841                 }
2842             }
2843 
2844             return q;
2845         },
2846 
2847         /**
2848          * Find the n-th intersection point of two curves named red (first parameter) and blue (second parameter).
2849          * We go through each segment of the red curve and search if there is an intersection with a segment of the blue curve.
2850          * This double loop, i.e. the outer loop runs along the red curve and the inner loop runs along the blue curve, defines
2851          * the n-th intersection point. The segments are either line segments or Bezier curves of degree 3. This depends on
2852          * the property bezierDegree of the curves.
2853          * <p>
2854          * This method works also for transformed curves, since only the already
2855          * transformed points are used.
2856          *
2857          * @param {JXG.Curve} red
2858          * @param {JXG.Curve} blue
2859          * @param {Number|Function} nr
2860          */
2861         meetCurveRedBlueSegments: function (red, blue, nr) {
2862             var i,
2863                 j,
2864                 n = Type.evaluate(nr),
2865                 red1,
2866                 red2,
2867                 blue1,
2868                 blue2,
2869                 m,
2870                 minX,
2871                 maxX,
2872                 iFound = 0,
2873                 lenBlue = blue.numberPoints,
2874                 lenRed = red.numberPoints;
2875 
2876             if (lenBlue <= 1 || lenRed <= 1) {
2877                 return [0, NaN, NaN];
2878             }
2879 
2880             for (i = 1; i < lenRed; i++) {
2881                 red1 = red.points[i - 1].usrCoords;
2882                 red2 = red.points[i].usrCoords;
2883                 minX = Math.min(red1[1], red2[1]);
2884                 maxX = Math.max(red1[1], red2[1]);
2885 
2886                 blue2 = blue.points[0].usrCoords;
2887                 for (j = 1; j < lenBlue; j++) {
2888                     blue1 = blue2;
2889                     blue2 = blue.points[j].usrCoords;
2890                     if (
2891                         Math.min(blue1[1], blue2[1]) < maxX &&
2892                         Math.max(blue1[1], blue2[1]) > minX
2893                     ) {
2894                         m = this.meetSegmentSegment(red1, red2, blue1, blue2);
2895                         if (
2896                             m[1] >= 0.0 && m[2] >= 0.0 &&
2897                             // The two segments meet in the interior or at the start points
2898                             ((m[1] < 1.0 && m[2] < 1.0) ||
2899                               // One of the curve is intersected in the very last point
2900                                 (i === lenRed - 1 && m[1] === 1.0) ||
2901                                 (j === lenBlue - 1 && m[2] === 1.0))
2902                         ) {
2903                             if (iFound === n) {
2904                                 return m[0];
2905                             }
2906 
2907                             iFound++;
2908                         }
2909                     }
2910                 }
2911             }
2912 
2913             return [0, NaN, NaN];
2914         },
2915 
2916         /**
2917          * (Virtual) Intersection of two segments.
2918          * @param {Array} p1 First point of segment 1 using normalized homogeneous coordinates [1,x,y]
2919          * @param {Array} p2 Second point or direction of segment 1 using normalized homogeneous coordinates [1,x,y] or point at infinity [0,x,y], respectively
2920          * @param {Array} q1 First point of segment 2 using normalized homogeneous coordinates [1,x,y]
2921          * @param {Array} q2 Second point or direction of segment 2 using normalized homogeneous coordinates [1,x,y] or point at infinity [0,x,y], respectively
2922          * @returns {Array} [Intersection point, t, u] The first entry contains the homogeneous coordinates
2923          * of the intersection point. The second and third entry give the position of the intersection with respect
2924          * to the definiting parameters. For example, the second entry t is defined by: intersection point = p1 + t * deltaP, where
2925          * deltaP = (p2 - p1) when both parameters are coordinates, and deltaP = p2 if p2 is a point at infinity.
2926          * If the two segments are collinear, [[0,0,0], Infinity, Infinity] is returned.
2927          **/
2928         meetSegmentSegment: function (p1, p2, q1, q2) {
2929             var t,
2930                 u,
2931                 i,
2932                 d,
2933                 li1 = Mat.crossProduct(p1, p2),
2934                 li2 = Mat.crossProduct(q1, q2),
2935                 c = Mat.crossProduct(li1, li2);
2936 
2937             if (Math.abs(c[0]) < Mat.eps) {
2938                 return [c, Infinity, Infinity];
2939             }
2940 
2941             // Normalize the intersection coordinates
2942             c[1] /= c[0];
2943             c[2] /= c[0];
2944             c[0] /= c[0];
2945 
2946             // Now compute in principle:
2947             //    t = dist(c - p1) / dist(p2 - p1) and
2948             //    u = dist(c - q1) / dist(q2 - q1)
2949             // However: the points q1, q2, p1, p2 might be ideal points - or in general - the
2950             // coordinates might be not normalized.
2951             // Note that the z-coordinates of p2 and q2 are used to determine whether it should be interpreted
2952             // as a segment coordinate or a direction.
2953             i = Math.abs(p2[1] - p2[0] * p1[1]) < Mat.eps ? 2 : 1;
2954             d = p1[i] / p1[0];
2955             t = (c[i] - d) / (p2[0] !== 0 ? p2[i] / p2[0] - d : p2[i]);
2956 
2957             i = Math.abs(q2[1] - q2[0] * q1[1]) < Mat.eps ? 2 : 1;
2958             d = q1[i] / q1[0];
2959             u = (c[i] - d) / (q2[0] !== 0 ? q2[i] / q2[0] - d : q2[i]);
2960 
2961             return [c, t, u];
2962         },
2963 
2964         /**
2965          * Find the n-th intersection point of two pathes, usually given by polygons. Uses parts of the
2966          * Greiner-Hormann algorithm in JXG.Math.Clip.
2967          *
2968          * @param {JXG.Circle|JXG.Curve|JXG.Polygon} path1
2969          * @param {JXG.Circle|JXG.Curve|JXG.Polygon} path2
2970          * @param {Number|Function} n
2971          * @param {JXG.Board} board
2972          *
2973          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
2974          * the ideal point [0,0,0] is returned.
2975          *
2976          */
2977         meetPathPath: function (path1, path2, nr, board) {
2978             var S, C, len, intersections,
2979                 n = Type.evaluate(nr);
2980 
2981             S = JXG.Math.Clip._getPath(path1, board);
2982             len = S.length;
2983             if (
2984                 len > 0 &&
2985                 this.distance(S[0].coords.usrCoords, S[len - 1].coords.usrCoords, 3) < Mat.eps
2986             ) {
2987                 S.pop();
2988             }
2989 
2990             C = JXG.Math.Clip._getPath(path2, board);
2991             len = C.length;
2992             if (
2993                 len > 0 &&
2994                 this.distance(C[0].coords.usrCoords, C[len - 1].coords.usrCoords, 3) <
2995                 Mat.eps * Mat.eps
2996             ) {
2997                 C.pop();
2998             }
2999 
3000             // Handle cases where at least one of the paths is empty
3001             if (nr < 0 || JXG.Math.Clip.isEmptyCase(S, C, 'intersection')) {
3002                 return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
3003             }
3004 
3005             JXG.Math.Clip.makeDoublyLinkedList(S);
3006             JXG.Math.Clip.makeDoublyLinkedList(C);
3007 
3008             intersections = JXG.Math.Clip.findIntersections(S, C, board)[0];
3009             if (n < intersections.length) {
3010                 return intersections[n].coords;
3011             }
3012             return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
3013         },
3014 
3015         /**
3016          * Find the n-th intersection point between a polygon and a line.
3017          * @param {JXG.Polygon} path
3018          * @param {JXG.Line} line
3019          * @param {Number|Function} nr
3020          * @param {JXG.Board} board
3021          * @param {Boolean} alwaysIntersect If false just the segment between the two defining points of the line are tested for intersection.
3022          *
3023          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
3024          * the ideal point [0,0,0] is returned.
3025          */
3026         meetPolygonLine: function (path, line, nr, board, alwaysIntersect) {
3027             var i,
3028                 n = Type.evaluate(nr),
3029                 res,
3030                 border,
3031                 crds = [0, 0, 0],
3032                 len = path.borders.length,
3033                 intersections = [];
3034 
3035             for (i = 0; i < len; i++) {
3036                 border = path.borders[i];
3037                 res = this.meetSegmentSegment(
3038                     border.point1.coords.usrCoords,
3039                     border.point2.coords.usrCoords,
3040                     line.point1.coords.usrCoords,
3041                     line.point2.coords.usrCoords
3042                 );
3043 
3044                 if (
3045                     (!alwaysIntersect || (res[2] >= 0 && res[2] < 1)) &&
3046                     res[1] >= 0 &&
3047                     res[1] < 1
3048                 ) {
3049                     intersections.push(res[0]);
3050                 }
3051             }
3052 
3053             if (n >= 0 && n < intersections.length) {
3054                 crds = intersections[n];
3055             }
3056             return new Coords(Const.COORDS_BY_USER, crds, board);
3057         },
3058 
3059         /****************************************/
3060         /****   BEZIER CURVE ALGORITHMS      ****/
3061         /****************************************/
3062 
3063         /**
3064          * Splits a Bezier curve segment defined by four points into
3065          * two Bezier curve segments. Dissection point is t=1/2.
3066          * @param {Array} curve Array of four coordinate arrays of length 2 defining a
3067          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3068          * @returns {Array} Array consisting of two coordinate arrays for Bezier curves.
3069          */
3070         _bezierSplit: function (curve) {
3071             var p0, p1, p2, p00, p22, p000;
3072 
3073             p0 = [(curve[0][0] + curve[1][0]) * 0.5, (curve[0][1] + curve[1][1]) * 0.5];
3074             p1 = [(curve[1][0] + curve[2][0]) * 0.5, (curve[1][1] + curve[2][1]) * 0.5];
3075             p2 = [(curve[2][0] + curve[3][0]) * 0.5, (curve[2][1] + curve[3][1]) * 0.5];
3076 
3077             p00 = [(p0[0] + p1[0]) * 0.5, (p0[1] + p1[1]) * 0.5];
3078             p22 = [(p1[0] + p2[0]) * 0.5, (p1[1] + p2[1]) * 0.5];
3079 
3080             p000 = [(p00[0] + p22[0]) * 0.5, (p00[1] + p22[1]) * 0.5];
3081 
3082             return [
3083                 [curve[0], p0, p00, p000],
3084                 [p000, p22, p2, curve[3]]
3085             ];
3086         },
3087 
3088         /**
3089          * Computes the bounding box [minX, maxY, maxX, minY] of a Bezier curve segment
3090          * from its control points.
3091          * @param {Array} curve Array of four coordinate arrays of length 2 defining a
3092          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3093          * @returns {Array} Bounding box [minX, maxY, maxX, minY]
3094          */
3095         _bezierBbox: function (curve) {
3096             var bb = [];
3097 
3098             if (curve.length === 4) {
3099                 // bezierDegree == 3
3100                 bb[0] = Math.min(curve[0][0], curve[1][0], curve[2][0], curve[3][0]); // minX
3101                 bb[1] = Math.max(curve[0][1], curve[1][1], curve[2][1], curve[3][1]); // maxY
3102                 bb[2] = Math.max(curve[0][0], curve[1][0], curve[2][0], curve[3][0]); // maxX
3103                 bb[3] = Math.min(curve[0][1], curve[1][1], curve[2][1], curve[3][1]); // minY
3104             } else {
3105                 // bezierDegree == 1
3106                 bb[0] = Math.min(curve[0][0], curve[1][0]); // minX
3107                 bb[1] = Math.max(curve[0][1], curve[1][1]); // maxY
3108                 bb[2] = Math.max(curve[0][0], curve[1][0]); // maxX
3109                 bb[3] = Math.min(curve[0][1], curve[1][1]); // minY
3110             }
3111 
3112             return bb;
3113         },
3114 
3115         /**
3116          * Decide if two Bezier curve segments overlap by comparing their bounding boxes.
3117          * @param {Array} bb1 Bounding box of the first Bezier curve segment
3118          * @param {Array} bb2 Bounding box of the second Bezier curve segment
3119          * @returns {Boolean} true if the bounding boxes overlap, false otherwise.
3120          */
3121         _bezierOverlap: function (bb1, bb2) {
3122             return bb1[2] >= bb2[0] && bb1[0] <= bb2[2] && bb1[1] >= bb2[3] && bb1[3] <= bb2[1];
3123         },
3124 
3125         /**
3126          * Append list of intersection points to a list.
3127          * @private
3128          */
3129         _bezierListConcat: function (L, Lnew, t1, t2) {
3130             var i,
3131                 t2exists = Type.exists(t2),
3132                 start = 0,
3133                 len = Lnew.length,
3134                 le = L.length;
3135 
3136             if (
3137                 le > 0 &&
3138                 len > 0 &&
3139                 ((L[le - 1][1] === 1 && Lnew[0][1] === 0) ||
3140                     (t2exists && L[le - 1][2] === 1 && Lnew[0][2] === 0))
3141             ) {
3142                 start = 1;
3143             }
3144 
3145             for (i = start; i < len; i++) {
3146                 if (t2exists) {
3147                     Lnew[i][2] *= 0.5;
3148                     Lnew[i][2] += t2;
3149                 }
3150 
3151                 Lnew[i][1] *= 0.5;
3152                 Lnew[i][1] += t1;
3153 
3154                 L.push(Lnew[i]);
3155             }
3156         },
3157 
3158         /**
3159          * Find intersections of two Bezier curve segments by recursive subdivision.
3160          * Below maxlevel determine intersections by intersection line segments.
3161          * @param {Array} red Array of four coordinate arrays of length 2 defining the first
3162          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3163          * @param {Array} blue Array of four coordinate arrays of length 2 defining the second
3164          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3165          * @param {Number} level Recursion level
3166          * @returns {Array} List of intersection points (up to nine). Each intersection point is an
3167          * array of length three (homogeneous coordinates) plus preimages.
3168          */
3169         _bezierMeetSubdivision: function (red, blue, level) {
3170             var bbb,
3171                 bbr,
3172                 ar,
3173                 b0,
3174                 b1,
3175                 r0,
3176                 r1,
3177                 m,
3178                 p0,
3179                 p1,
3180                 q0,
3181                 q1,
3182                 L = [],
3183                 maxLev = 5; // Maximum recursion level
3184 
3185             bbr = this._bezierBbox(blue);
3186             bbb = this._bezierBbox(red);
3187 
3188             if (!this._bezierOverlap(bbr, bbb)) {
3189                 return [];
3190             }
3191 
3192             if (level < maxLev) {
3193                 ar = this._bezierSplit(red);
3194                 r0 = ar[0];
3195                 r1 = ar[1];
3196 
3197                 ar = this._bezierSplit(blue);
3198                 b0 = ar[0];
3199                 b1 = ar[1];
3200 
3201                 this._bezierListConcat(
3202                     L,
3203                     this._bezierMeetSubdivision(r0, b0, level + 1),
3204                     0.0,
3205                     0.0
3206                 );
3207                 this._bezierListConcat(
3208                     L,
3209                     this._bezierMeetSubdivision(r0, b1, level + 1),
3210                     0,
3211                     0.5
3212                 );
3213                 this._bezierListConcat(
3214                     L,
3215                     this._bezierMeetSubdivision(r1, b0, level + 1),
3216                     0.5,
3217                     0.0
3218                 );
3219                 this._bezierListConcat(
3220                     L,
3221                     this._bezierMeetSubdivision(r1, b1, level + 1),
3222                     0.5,
3223                     0.5
3224                 );
3225 
3226                 return L;
3227             }
3228 
3229             // Make homogeneous coordinates
3230             q0 = [1].concat(red[0]);
3231             q1 = [1].concat(red[3]);
3232             p0 = [1].concat(blue[0]);
3233             p1 = [1].concat(blue[3]);
3234 
3235             m = this.meetSegmentSegment(q0, q1, p0, p1);
3236 
3237             if (m[1] >= 0.0 && m[2] >= 0.0 && m[1] <= 1.0 && m[2] <= 1.0) {
3238                 return [m];
3239             }
3240 
3241             return [];
3242         },
3243 
3244         /**
3245          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the line defined by the segment
3246          */
3247         _bezierLineMeetSubdivision: function (red, blue, level, testSegment) {
3248             var bbb, bbr, ar,
3249                 r0, r1,
3250                 m,
3251                 p0, p1, q0, q1,
3252                 L = [],
3253                 maxLev = 5; // Maximum recursion level
3254 
3255             bbb = this._bezierBbox(blue);
3256             bbr = this._bezierBbox(red);
3257 
3258             if (testSegment && !this._bezierOverlap(bbr, bbb)) {
3259                 return [];
3260             }
3261 
3262             if (level < maxLev) {
3263                 ar = this._bezierSplit(red);
3264                 r0 = ar[0];
3265                 r1 = ar[1];
3266 
3267                 this._bezierListConcat(
3268                     L,
3269                     this._bezierLineMeetSubdivision(r0, blue, level + 1),
3270                     0.0
3271                 );
3272                 this._bezierListConcat(
3273                     L,
3274                     this._bezierLineMeetSubdivision(r1, blue, level + 1),
3275                     0.5
3276                 );
3277 
3278                 return L;
3279             }
3280 
3281             // Make homogeneous coordinates
3282             q0 = [1].concat(red[0]);
3283             q1 = [1].concat(red[3]);
3284             p0 = [1].concat(blue[0]);
3285             p1 = [1].concat(blue[1]);
3286 
3287             m = this.meetSegmentSegment(q0, q1, p0, p1);
3288 
3289             if (m[1] >= 0.0 && m[1] <= 1.0) {
3290                 if (!testSegment || (m[2] >= 0.0 && m[2] <= 1.0)) {
3291                     return [m];
3292                 }
3293             }
3294 
3295             return [];
3296         },
3297 
3298         /**
3299          * Find the nr-th intersection point of two Bezier curve segments.
3300          * @param {Array} red Array of four coordinate arrays of length 2 defining the first
3301          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3302          * @param {Array} blue Array of four coordinate arrays of length 2 defining the second
3303          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
3304          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the line defined by the segment
3305          * @returns {Array} Array containing the list of all intersection points as homogeneous coordinate arrays plus
3306          * preimages [x,y], t_1, t_2] of the two Bezier curve segments.
3307          *
3308          */
3309         meetBeziersegmentBeziersegment: function (red, blue, testSegment) {
3310             var L, L2, i;
3311 
3312             if (red.length === 4 && blue.length === 4) {
3313                 L = this._bezierMeetSubdivision(red, blue, 0);
3314             } else {
3315                 L = this._bezierLineMeetSubdivision(red, blue, 0, testSegment);
3316             }
3317 
3318             L.sort(function (a, b) {
3319                 return (a[1] - b[1]) * 10000000.0 + (a[2] - b[2]);
3320             });
3321 
3322             L2 = [];
3323             for (i = 0; i < L.length; i++) {
3324                 // Only push entries different from their predecessor
3325                 if (i === 0 || L[i][1] !== L[i - 1][1] || L[i][2] !== L[i - 1][2]) {
3326                     L2.push(L[i]);
3327                 }
3328             }
3329             return L2;
3330         },
3331 
3332         /**
3333          * Find the nr-th intersection point of two Bezier curves, i.e. curves with bezierDegree == 3.
3334          * @param {JXG.Curve} red Curve with bezierDegree == 3
3335          * @param {JXG.Curve} blue Curve with bezierDegree == 3
3336          * @param {Number|Function} nr The number of the intersection point which should be returned.
3337          * @returns {Array} The homogeneous coordinates of the nr-th intersection point.
3338          */
3339         meetBezierCurveRedBlueSegments: function (red, blue, nr) {
3340             var p, i, j, k,
3341                 n = Type.evaluate(nr),
3342                 po, tmp,
3343                 redArr,
3344                 blueArr,
3345                 bbr,
3346                 bbb,
3347                 intersections,
3348                 startRed = 0,
3349                 startBlue = 0,
3350                 lenBlue, lenRed,
3351                 L = [];
3352 
3353             if (blue.numberPoints < blue.bezierDegree + 1 || red.numberPoints < red.bezierDegree + 1) {
3354                 return [0, NaN, NaN];
3355             }
3356             if (red.bezierDegree === 1 && blue.bezierDegree === 3) {
3357                 tmp = red;
3358                 red = blue;
3359                 blue = tmp;
3360             }
3361 
3362             lenBlue = blue.numberPoints - blue.bezierDegree;
3363             lenRed = red.numberPoints - red.bezierDegree;
3364 
3365             // For sectors, we ignore the "legs"
3366             if (red.type === Const.OBJECT_TYPE_SECTOR) {
3367                 startRed = 3;
3368                 lenRed -= 3;
3369             }
3370             if (blue.type === Const.OBJECT_TYPE_SECTOR) {
3371                 startBlue = 3;
3372                 lenBlue -= 3;
3373             }
3374 
3375             for (i = startRed; i < lenRed; i += red.bezierDegree) {
3376                 p = red.points;
3377                 redArr = [p[i].usrCoords.slice(1), p[i + 1].usrCoords.slice(1)];
3378                 if (red.bezierDegree === 3) {
3379                     redArr[2] = p[i + 2].usrCoords.slice(1);
3380                     redArr[3] = p[i + 3].usrCoords.slice(1);
3381                 }
3382 
3383                 bbr = this._bezierBbox(redArr);
3384 
3385                 for (j = startBlue; j < lenBlue; j += blue.bezierDegree) {
3386                     p = blue.points;
3387                     blueArr = [p[j].usrCoords.slice(1), p[j + 1].usrCoords.slice(1)];
3388                     if (blue.bezierDegree === 3) {
3389                         blueArr[2] = p[j + 2].usrCoords.slice(1);
3390                         blueArr[3] = p[j + 3].usrCoords.slice(1);
3391                     }
3392 
3393                     bbb = this._bezierBbox(blueArr);
3394                     if (this._bezierOverlap(bbr, bbb)) {
3395                         intersections = this.meetBeziersegmentBeziersegment(redArr, blueArr);
3396                         if (intersections.length === 0) {
3397                             continue;
3398                         }
3399                         for (k = 0; k < intersections.length; k++) {
3400                             po = intersections[k];
3401                             if (
3402                                 po[1] < -Mat.eps ||
3403                                 po[1] > 1 + Mat.eps ||
3404                                 po[2] < -Mat.eps ||
3405                                 po[2] > 1 + Mat.eps
3406                             ) {
3407                                 continue;
3408                             }
3409                             L.push(po);
3410                         }
3411                         if (L.length > n) {
3412                             return L[n][0];
3413                         }
3414                     }
3415                 }
3416             }
3417             if (L.length > n) {
3418                 return L[n][0];
3419             }
3420 
3421             return [0, NaN, NaN];
3422         },
3423 
3424         bezierSegmentEval: function (t, curve) {
3425             var f,
3426                 x,
3427                 y,
3428                 t1 = 1.0 - t;
3429 
3430             x = 0;
3431             y = 0;
3432 
3433             f = t1 * t1 * t1;
3434             x += f * curve[0][0];
3435             y += f * curve[0][1];
3436 
3437             f = 3.0 * t * t1 * t1;
3438             x += f * curve[1][0];
3439             y += f * curve[1][1];
3440 
3441             f = 3.0 * t * t * t1;
3442             x += f * curve[2][0];
3443             y += f * curve[2][1];
3444 
3445             f = t * t * t;
3446             x += f * curve[3][0];
3447             y += f * curve[3][1];
3448 
3449             return [1.0, x, y];
3450         },
3451 
3452         /**
3453          * Generate the defining points of a 3rd degree bezier curve that approximates
3454          * a circle sector defined by three coordinate points A, B, C, each defined by an array of length three.
3455          * The coordinate arrays are given in homogeneous coordinates.
3456          * @param {Array} A First point
3457          * @param {Array} B Second point (intersection point)
3458          * @param {Array} C Third point
3459          * @param {Boolean} withLegs Flag. If true the legs to the intersection point are part of the curve.
3460          * @param {Number} sgn Wither 1 or -1. Needed for minor and major arcs. In case of doubt, use 1.
3461          */
3462         bezierArc: function (A, B, C, withLegs, sgn) {
3463             var p1, p2, p3, p4,
3464                 r,
3465                 phi, beta, delta,
3466                 // PI2 = Math.PI * 0.5,
3467                 x = B[1],
3468                 y = B[2],
3469                 z = B[0],
3470                 dataX = [],
3471                 dataY = [],
3472                 co, si,
3473                 ax, ay,
3474                 bx, by,
3475                 k, v, d,
3476                 matrix;
3477 
3478             r = this.distance(B, A);
3479 
3480             // x,y, z is intersection point. Normalize it.
3481             x /= z;
3482             y /= z;
3483 
3484             phi = this.rad(A.slice(1), B.slice(1), C.slice(1));
3485             if (sgn === -1) {
3486                 phi = 2 * Math.PI - phi;
3487             }
3488 
3489             // Always divide the arc into four Bezier arcs.
3490             // Otherwise, the position of gliders on this arc
3491             // will be wrong.
3492             delta = phi / 4;
3493 
3494 
3495             p1 = A;
3496             p1[1] /= p1[0];
3497             p1[2] /= p1[0];
3498             p1[0] /= p1[0];
3499 
3500             p4 = p1.slice(0);
3501 
3502             if (withLegs) {
3503                 dataX = [x, x + 0.333 * (p1[1] - x), x + 0.666 * (p1[1] - x), p1[1]];
3504                 dataY = [y, y + 0.333 * (p1[2] - y), y + 0.666 * (p1[2] - y), p1[2]];
3505             } else {
3506                 dataX = [p1[1]];
3507                 dataY = [p1[2]];
3508             }
3509 
3510             while (phi > Mat.eps) {
3511                 // if (phi > PI2) {
3512                 //     beta = PI2;
3513                 //     phi -= PI2;
3514                 // } else {
3515                 //     beta = phi;
3516                 //     phi = 0;
3517                 // }
3518                 if (phi > delta) {
3519                     beta = delta;
3520                     phi -= delta;
3521                 } else {
3522                     beta = phi;
3523                     phi = 0;
3524                 }
3525 
3526                 co = Math.cos(sgn * beta);
3527                 si = Math.sin(sgn * beta);
3528 
3529                 matrix = [
3530                     [1, 0, 0],
3531                     [x * (1 - co) + y * si, co, -si],
3532                     [y * (1 - co) - x * si, si, co]
3533                 ];
3534                 v = Mat.matVecMult(matrix, p1);
3535                 p4 = [v[0] / v[0], v[1] / v[0], v[2] / v[0]];
3536 
3537                 ax = p1[1] - x;
3538                 ay = p1[2] - y;
3539                 bx = p4[1] - x;
3540                 by = p4[2] - y;
3541                 d = Mat.hypot(ax + bx, ay + by);
3542 
3543                 if (Math.abs(by - ay) > Mat.eps) {
3544                     k = ((((ax + bx) * (r / d - 0.5)) / (by - ay)) * 8) / 3;
3545                 } else {
3546                     k = ((((ay + by) * (r / d - 0.5)) / (ax - bx)) * 8) / 3;
3547                 }
3548 
3549                 p2 = [1, p1[1] - k * ay, p1[2] + k * ax];
3550                 p3 = [1, p4[1] + k * by, p4[2] - k * bx];
3551 
3552                 Type.concat(dataX, [p2[1], p3[1], p4[1]]);
3553                 Type.concat(dataY, [p2[2], p3[2], p4[2]]);
3554                 p1 = p4.slice(0);
3555             }
3556 
3557             if (withLegs) {
3558                 Type.concat(dataX, [
3559                     p4[1] + 0.333 * (x - p4[1]),
3560                     p4[1] + 0.666 * (x - p4[1]),
3561                     x
3562                 ]);
3563                 Type.concat(dataY, [
3564                     p4[2] + 0.333 * (y - p4[2]),
3565                     p4[2] + 0.666 * (y - p4[2]),
3566                     y
3567                 ]);
3568             }
3569 
3570             return [dataX, dataY];
3571         },
3572 
3573         /****************************************/
3574         /****           PROJECTIONS          ****/
3575         /****************************************/
3576 
3577         /**
3578          * Calculates the coordinates of the projection of a given point on a given circle. I.o.w. the
3579          * nearest one of the two intersection points of the line through the given point and the circles
3580          * center.
3581          * @param {JXG.Point|JXG.Coords} point Point to project or coords object to project.
3582          * @param {JXG.Circle} circle Circle on that the point is projected.
3583          * @param {JXG.Board} [board=point.board] Reference to the board
3584          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given circle.
3585          */
3586         projectPointToCircle: function (point, circle, board) {
3587             var dist,
3588                 P,
3589                 x,
3590                 y,
3591                 factor,
3592                 M = circle.center.coords.usrCoords;
3593 
3594             if (!Type.exists(board)) {
3595                 board = point.board;
3596             }
3597 
3598             // gave us a point
3599             if (Type.isPoint(point)) {
3600                 dist = point.coords.distance(Const.COORDS_BY_USER, circle.center.coords);
3601                 P = point.coords.usrCoords;
3602                 // gave us coords
3603             } else {
3604                 dist = point.distance(Const.COORDS_BY_USER, circle.center.coords);
3605                 P = point.usrCoords;
3606             }
3607 
3608             if (Math.abs(dist) < Mat.eps) {
3609                 dist = Mat.eps;
3610             }
3611 
3612             factor = circle.Radius() / dist;
3613             x = M[1] + factor * (P[1] - M[1]);
3614             y = M[2] + factor * (P[2] - M[2]);
3615 
3616             return new Coords(Const.COORDS_BY_USER, [x, y], board);
3617         },
3618 
3619         /**
3620          * Calculates the coordinates of the orthogonal projection of a given point on a given line. I.o.w. the
3621          * intersection point of the given line and its perpendicular through the given point.
3622          * @param {JXG.Point|JXG.Coords} point Point to project.
3623          * @param {JXG.Line} line Line on that the point is projected.
3624          * @param {JXG.Board} [board=point.board|board=line.board] Reference to a board.
3625          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given line.
3626          */
3627         projectPointToLine: function (point, line, board) {
3628             var v = [0, line.stdform[1], line.stdform[2]],
3629                 coords;
3630 
3631             if (!Type.exists(board)) {
3632                 if (Type.exists(point.coords)) {
3633                     board = point.board;
3634                 } else {
3635                     board = line.board;
3636                 }
3637             }
3638 
3639             if (Type.exists(point.coords)) {
3640                 coords = point.coords.usrCoords;
3641             } else {
3642                 coords = point.usrCoords;
3643             }
3644 
3645             v = Mat.crossProduct(v, coords);
3646             return new Coords(Const.COORDS_BY_USER, Mat.crossProduct(v, line.stdform), board);
3647         },
3648 
3649         /**
3650          * Calculates the coordinates of the orthogonal projection of a given coordinate array on a given line
3651          * segment defined by two coordinate arrays.
3652          * @param {Array} p Point to project.
3653          * @param {Array} q1 Start point of the line segment on that the point is projected.
3654          * @param {Array} q2 End point of the line segment on that the point is projected.
3655          * @returns {Array} The coordinates of the projection of the given point on the given segment
3656          * and the factor that determines the projected point as a convex combination of the
3657          * two endpoints q1 and q2 of the segment.
3658          */
3659         projectCoordsToSegment: function (p, q1, q2) {
3660             var t,
3661                 denom,
3662                 s = [q2[1] - q1[1], q2[2] - q1[2]],
3663                 v = [p[1] - q1[1], p[2] - q1[2]];
3664 
3665             /**
3666              * If the segment has length 0, i.e. is a point,
3667              * the projection is equal to that point.
3668              */
3669             if (Math.abs(s[0]) < Mat.eps && Math.abs(s[1]) < Mat.eps) {
3670                 return [q1, 0];
3671             }
3672 
3673             t = Mat.innerProduct(v, s);
3674             denom = Mat.innerProduct(s, s);
3675             t /= denom;
3676 
3677             return [[1, t * s[0] + q1[1], t * s[1] + q1[2]], t];
3678         },
3679 
3680         /**
3681          * Finds the coordinates of the closest point on a Bezier segment of a
3682          * {@link JXG.Curve} to a given coordinate array.
3683          * @param {Array} pos Point to project in homogeneous coordinates.
3684          * @param {JXG.Curve} curve Curve of type "plot" having Bezier degree 3.
3685          * @param {Number} start Number of the Bezier segment of the curve.
3686          * @returns {Array} The coordinates of the projection of the given point
3687          * on the given Bezier segment and the preimage of the curve which
3688          * determines the closest point.
3689          */
3690         projectCoordsToBeziersegment: function (pos, curve, start) {
3691             var t0,
3692                 /** @ignore */
3693                 minfunc = function (t) {
3694                     var z = [1, curve.X(start + t), curve.Y(start + t)];
3695 
3696                     z[1] -= pos[1];
3697                     z[2] -= pos[2];
3698 
3699                     return z[1] * z[1] + z[2] * z[2];
3700                 };
3701 
3702             t0 = JXG.Math.Numerics.fminbr(minfunc, [0.0, 1.0]);
3703 
3704             return [[1, curve.X(t0 + start), curve.Y(t0 + start)], t0];
3705         },
3706 
3707         /**
3708          * Calculates the coordinates of the projection of a given point on a given curve.
3709          * Uses {@link JXG.Math.Geometry.projectCoordsToCurve}.
3710          *
3711          * @param {JXG.Point} point Point to project.
3712          * @param {JXG.Curve} curve Curve on that the point is projected.
3713          * @param {JXG.Board} [board=point.board] Reference to a board.
3714          * @see JXG.Math.Geometry.projectCoordsToCurve
3715          * @returns {Array} [JXG.Coords, position] The coordinates of the projection of the given
3716          * point on the given graph and the relative position on the curve (real number).
3717          */
3718         projectPointToCurve: function (point, curve, board) {
3719             if (!Type.exists(board)) {
3720                 board = point.board;
3721             }
3722 
3723             var x = point.X(),
3724                 y = point.Y(),
3725                 t = point.position,
3726                 result;
3727 
3728             if (!Type.exists(t)) {
3729                 t = curve.evalVisProp('curvetype') === 'functiongraph' ? x : 0.0;
3730             }
3731             result = this.projectCoordsToCurve(x, y, t, curve, board);
3732             // point.position = result[1];
3733 
3734             return result;
3735         },
3736 
3737         /**
3738          * Calculates the coordinates of the projection of a coordinates pair on a given curve. In case of
3739          * function graphs this is the
3740          * intersection point of the curve and the parallel to y-axis through the given point.
3741          * @param {Number} x coordinate to project.
3742          * @param {Number} y coordinate to project.
3743          * @param {Number} t start value for newtons method
3744          * @param {JXG.Curve} curve Curve on that the point is projected.
3745          * @param {JXG.Board} [board=curve.board] Reference to a board.
3746          * @see JXG.Math.Geometry.projectPointToCurve
3747          * @returns {JXG.Coords} Array containing the coordinates of the projection of the given point on the given curve and
3748          * the position on the curve.
3749          */
3750         projectCoordsToCurve: function (x, y, t, curve, board) {
3751             var newCoords, newCoordsObj,
3752                 i, j, mindist, dist, lbda,
3753                 v, coords, d, p1, p2, res, minfunc,
3754                 t_new, f_new, f_old, dy,
3755                 delta, delta1, delta2, steps,
3756                 minX, maxX, minX_glob, maxX_glob,
3757                 infty = Number.POSITIVE_INFINITY;
3758 
3759             if (!Type.exists(board)) {
3760                 board = curve.board;
3761             }
3762 
3763             if (curve.evalVisProp('curvetype') === 'plot') {
3764                 t = 0;
3765                 mindist = infty;
3766                 if (curve.numberPoints === 0) {
3767                     newCoords = [0, 1, 1];
3768                 } else {
3769                     newCoords = [curve.Z(0), curve.X(0), curve.Y(0)];
3770                 }
3771 
3772                 if (curve.numberPoints > 1) {
3773                     v = [1, x, y];
3774                     if (curve.bezierDegree === 3) {
3775                         j = 0;
3776                     } else {
3777                         p1 = [curve.Z(0), curve.X(0), curve.Y(0)];
3778                     }
3779                     for (i = 0; i < curve.numberPoints - 1; i++) {
3780                         if (curve.bezierDegree === 3) {
3781                             res = this.projectCoordsToBeziersegment(v, curve, j);
3782                         } else {
3783                             p2 = [curve.Z(i + 1), curve.X(i + 1), curve.Y(i + 1)];
3784                             res = this.projectCoordsToSegment(v, p1, p2);
3785                         }
3786                         lbda = res[1];
3787                         coords = res[0];
3788 
3789                         if (0.0 <= lbda && lbda <= 1.0) {
3790                             dist = this.distance(coords, v);
3791                             d = i + lbda;
3792                         } else if (lbda < 0.0) {
3793                             coords = p1;
3794                             dist = this.distance(p1, v);
3795                             d = i;
3796                         } else if (lbda > 1.0 && i === curve.numberPoints - 2) {
3797                             coords = p2;
3798                             dist = this.distance(coords, v);
3799                             d = curve.numberPoints - 1;
3800                         }
3801 
3802                         if (dist < mindist) {
3803                             mindist = dist;
3804                             t = d;
3805                             newCoords = coords;
3806                         }
3807 
3808                         if (curve.bezierDegree === 3) {
3809                             j++;
3810                             i += 2;
3811                         } else {
3812                             p1 = p2;
3813                         }
3814                     }
3815                 }
3816 
3817                 newCoordsObj = new Coords(Const.COORDS_BY_USER, newCoords, board);
3818             } else {
3819                 // 'parameter', 'polar', 'functiongraph'
3820 
3821                 minX_glob = curve.minX();
3822                 maxX_glob = curve.maxX();
3823                 minX = minX_glob;
3824                 maxX = maxX_glob;
3825 
3826                 if (curve.evalVisProp('curvetype') === 'functiongraph') {
3827                     // Restrict the possible position of t
3828                     // to the projection of a circle to the x-axis (= t-axis)
3829                     dy = Math.abs(y - curve.Y(x));
3830                     if (!isNaN(dy)) {
3831                         minX = x - dy;
3832                         maxX = x + dy;
3833                     }
3834                 }
3835 
3836                 /**
3837                  * @ignore
3838                  * Find t such that the Euclidean distance between
3839                  * [x, y] and [curve.X(t), curve.Y(t)]
3840                  * is minimized.
3841                  */
3842                 minfunc = function (t) {
3843                     var dx, dy;
3844 
3845                     if (t < minX_glob || t > curve.maxX_glob) {
3846                         return Infinity;
3847                     }
3848                     dx = x - curve.X(t);
3849                     dy = y - curve.Y(t);
3850                     return dx * dx + dy * dy;
3851                 };
3852 
3853                 // Search t which minimizes minfunc(t)
3854                 // in discrete steps
3855                 f_old = minfunc(t);
3856                 steps = 50;
3857                 delta = (maxX - minX) / steps;
3858                 t_new = minX;
3859                 for (i = 0; i < steps; i++) {
3860                     f_new = minfunc(t_new);
3861 
3862                     if (f_new < f_old || f_old === Infinity || isNaN(f_old)) {
3863                         t = t_new;
3864                         f_old = f_new;
3865                     }
3866 
3867                     t_new += delta;
3868                 }
3869 
3870                 // t = Numerics.root(Numerics.D(minfunc), t);
3871 
3872                 // Ensure that minfunc is defined on the
3873                 // enclosing interval [t-delta1, t+delta2]
3874                 delta1 = delta;
3875                 for (i = 0; i < 20 && isNaN(minfunc(t - delta1)); i++, delta1 *= 0.5);
3876                 if (isNaN(minfunc(t - delta1))) {
3877                     delta1 = 0.0;
3878                 }
3879                 delta2 = delta;
3880                 for (i = 0; i < 20 && isNaN(minfunc(t + delta2)); i++, delta2 *= 0.5);
3881                 if (isNaN(minfunc(t + delta2))) {
3882                     delta2 = 0.0;
3883                 }
3884 
3885                 // Finally, apply mathemetical optimization in the determined interval
3886                 t = Numerics.fminbr(minfunc, [
3887                     Math.max(t - delta1, minX),
3888                     Math.min(t + delta2, maxX)
3889                 ]);
3890 
3891                 // Distinction between closed and open curves is not necessary.
3892                 // If closed, the cyclic projection shift will work anyhow
3893                 // if (Math.abs(curve.X(minX) - curve.X(maxX)) < Mat.eps &&
3894                 //     Math.abs(curve.Y(minX) - curve.Y(maxX)) < Mat.eps) {
3895                 //     // Cyclically
3896                 //     if (t < minX) {console.log(t)
3897                 //         t = maxX + t - minX;
3898                 //     }
3899                 //     if (t > maxX) {
3900                 //         t = minX + t - maxX;
3901                 //     }
3902                 // } else {
3903 
3904                 t = t < minX_glob ? minX_glob : t;
3905                 t = t > maxX_glob ? maxX_glob : t;
3906                 // }
3907 
3908                 newCoordsObj = new Coords(
3909                     Const.COORDS_BY_USER,
3910                     [curve.X(t), curve.Y(t)],
3911                     board
3912                 );
3913             }
3914 
3915             return [curve.updateTransform(newCoordsObj), t];
3916         },
3917 
3918         /**
3919          * Calculates the coordinates of the closest orthogonal projection of a given coordinate array onto the
3920          * border of a polygon.
3921          * @param {Array} p Point to project.
3922          * @param {JXG.Polygon} pol Polygon element
3923          * @returns {Array} The coordinates of the closest projection of the given point to the border of the polygon.
3924          */
3925         projectCoordsToPolygon: function (p, pol) {
3926             var i,
3927                 len = pol.vertices.length,
3928                 d_best = Infinity,
3929                 d,
3930                 projection,
3931                 proj,
3932                 bestprojection;
3933 
3934             for (i = 0; i < len - 1; i++) {
3935                 projection = JXG.Math.Geometry.projectCoordsToSegment(
3936                     p,
3937                     pol.vertices[i].coords.usrCoords,
3938                     pol.vertices[i + 1].coords.usrCoords
3939                 );
3940 
3941                 if (0 <= projection[1] && projection[1] <= 1) {
3942                     d = JXG.Math.Geometry.distance(projection[0], p, 3);
3943                     proj = projection[0];
3944                 } else if (projection[1] < 0) {
3945                     d = JXG.Math.Geometry.distance(pol.vertices[i].coords.usrCoords, p, 3);
3946                     proj = pol.vertices[i].coords.usrCoords;
3947                 } else {
3948                     d = JXG.Math.Geometry.distance(pol.vertices[i + 1].coords.usrCoords, p, 3);
3949                     proj = pol.vertices[i + 1].coords.usrCoords;
3950                 }
3951                 if (d < d_best) {
3952                     bestprojection = proj.slice(0);
3953                     d_best = d;
3954                 }
3955             }
3956             return bestprojection;
3957         },
3958 
3959         /**
3960          * Calculates the coordinates of the projection of a given point on a given turtle. A turtle consists of
3961          * one or more curves of curveType 'plot'. Uses {@link JXG.Math.Geometry.projectPointToCurve}.
3962          * @param {JXG.Point} point Point to project.
3963          * @param {JXG.Turtle} turtle on that the point is projected.
3964          * @param {JXG.Board} [board=point.board] Reference to a board.
3965          * @returns {Array} [JXG.Coords, position] Array containing the coordinates of the projection of the given point on the turtle and
3966          * the position on the turtle.
3967          */
3968         projectPointToTurtle: function (point, turtle, board) {
3969             var newCoords,
3970                 t,
3971                 x,
3972                 y,
3973                 i,
3974                 dist,
3975                 el,
3976                 minEl,
3977                 res,
3978                 newPos,
3979                 np = 0,
3980                 npmin = 0,
3981                 mindist = Number.POSITIVE_INFINITY,
3982                 len = turtle.objects.length;
3983 
3984             if (!Type.exists(board)) {
3985                 board = point.board;
3986             }
3987 
3988             // run through all curves of this turtle
3989             for (i = 0; i < len; i++) {
3990                 el = turtle.objects[i];
3991 
3992                 if (el.elementClass === Const.OBJECT_CLASS_CURVE) {
3993                     res = this.projectPointToCurve(point, el);
3994                     newCoords = res[0];
3995                     newPos = res[1];
3996                     dist = this.distance(newCoords.usrCoords, point.coords.usrCoords);
3997 
3998                     if (dist < mindist) {
3999                         x = newCoords.usrCoords[1];
4000                         y = newCoords.usrCoords[2];
4001                         t = newPos;
4002                         mindist = dist;
4003                         minEl = el;
4004                         npmin = np;
4005                     }
4006                     np += el.numberPoints;
4007                 }
4008             }
4009 
4010             newCoords = new Coords(Const.COORDS_BY_USER, [x, y], board);
4011             // point.position = t + npmin;
4012             // return minEl.updateTransform(newCoords);
4013             return [minEl.updateTransform(newCoords), t + npmin];
4014         },
4015 
4016         /**
4017          * Trivial projection of a point to another point.
4018          * @param {JXG.Point} point Point to project (not used).
4019          * @param {JXG.Point} dest Point on that the point is projected.
4020          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given circle.
4021          */
4022         projectPointToPoint: function (point, dest) {
4023             return dest.coords;
4024         },
4025 
4026         /**
4027          *
4028          * @param {JXG.Point|JXG.Coords} point
4029          * @param {JXG.Board} [board]
4030          */
4031         projectPointToBoard: function (point, board) {
4032             var i,
4033                 l,
4034                 c,
4035                 brd = board || point.board,
4036                 // comparison factor, point coord idx, bbox idx, 1st bbox corner x & y idx, 2nd bbox corner x & y idx
4037                 config = [
4038                     // left
4039                     [1, 1, 0, 0, 3, 0, 1],
4040                     // top
4041                     [-1, 2, 1, 0, 1, 2, 1],
4042                     // right
4043                     [-1, 1, 2, 2, 1, 2, 3],
4044                     // bottom
4045                     [1, 2, 3, 0, 3, 2, 3]
4046                 ],
4047                 coords = point.coords || point,
4048                 bbox = brd.getBoundingBox();
4049 
4050             for (i = 0; i < 4; i++) {
4051                 c = config[i];
4052                 if (c[0] * coords.usrCoords[c[1]] < c[0] * bbox[c[2]]) {
4053                     // define border
4054                     l = Mat.crossProduct(
4055                         [1, bbox[c[3]], bbox[c[4]]],
4056                         [1, bbox[c[5]], bbox[c[6]]]
4057                     );
4058                     l[3] = 0;
4059                     l = Mat.normalize(l);
4060 
4061                     // project point
4062                     coords = this.projectPointToLine({ coords: coords }, { stdform: l }, brd);
4063                 }
4064             }
4065 
4066             return coords;
4067         },
4068 
4069         /**
4070          * Calculates the distance of a point to a line. The point and the line are given by homogeneous
4071          * coordinates. For lines this can be line.stdform.
4072          * @param {Array} point Homogeneous coordinates of a point.
4073          * @param {Array} line Homogeneous coordinates of a line ([C,A,B] where A*x+B*y+C*z=0).
4074          * @returns {Number} Distance of the point to the line.
4075          */
4076         distPointLine: function (point, line) {
4077             var a = line[1],
4078                 b = line[2],
4079                 c = line[0],
4080                 nom;
4081 
4082             if (Math.abs(a) + Math.abs(b) < Mat.eps) {
4083                 return Number.POSITIVE_INFINITY;
4084             }
4085 
4086             nom = a * point[1] + b * point[2] + c;
4087             a *= a;
4088             b *= b;
4089 
4090             return Math.abs(nom) / Math.sqrt(a + b);
4091         },
4092 
4093         /**
4094          * Determine the (Euclidean) distance between a point q and a line segment
4095          * defined by two points p1 and p2. In case p1 equals p2, the distance to this
4096          * point is returned.
4097          *
4098          * @param {Array} q Homogeneous coordinates of q
4099          * @param {Array} p1 Homogeneous coordinates of p1
4100          * @param {Array} p2 Homogeneous coordinates of p2
4101          * @returns {Number} Distance of q to line segment [p1, p2]
4102          */
4103         distPointSegment: function (q, p1, p2) {
4104             var x, y, dx, dy,
4105                 den, lbda,
4106                 eps = Mat.eps * Mat.eps,
4107                 huge = 1000000;
4108 
4109             // Difference q - p1
4110             x = q[1] - p1[1];
4111             y = q[2] - p1[2];
4112             x = (x === Infinity) ? huge : (x === -Infinity) ? -huge : x;
4113             y = (y === Infinity) ? huge : (y === -Infinity) ? -huge : y;
4114 
4115             // Difference p2 - p1
4116             dx = p2[1] - p1[1];
4117             dy = p2[2] - p1[2];
4118             dx = (dx === Infinity) ? huge : (dx === -Infinity) ? -huge : dx;
4119             dy = (dy === Infinity) ? huge : (dy === -Infinity) ? -huge : dy;
4120 
4121             // If den==0 then p1 and p2 are identical
4122             // In this case the distance to p1 is returned
4123             den = dx * dx + dy * dy;
4124             if (den > eps) {
4125                 lbda = (x * dx + y * dy) / den;
4126                 if (lbda < 0.0) {
4127                     lbda = 0.0;
4128                 } else if (lbda > 1.0) {
4129                     lbda = 1.0;
4130                 }
4131                 x -= lbda * dx;
4132                 y -= lbda * dy;
4133             }
4134 
4135             return Mat.hypot(x, y);
4136         },
4137 
4138         /* ***************************************/
4139         /* *** 3D CALCULATIONS ****/
4140         /* ***************************************/
4141 
4142         /**
4143          * Generate the function which computes the data of the intersection between
4144          * <ul>
4145          * <li> plane3d, plane3d,
4146          * <li> plane3d, sphere3d,
4147          * <li> sphere3d, plane3d,
4148          * <li> sphere3d, sphere3d
4149          * </ul>
4150          *
4151          * @param {JXG.GeometryElement3D} el1 Plane or sphere element
4152          * @param {JXG.GeometryElement3D} el2 Plane or sphere element
4153          * @returns {Array} of functions needed as input to create the intersecting line or circle.
4154          *
4155          */
4156         intersectionFunction3D: function (view, el1, el2) {
4157             var func,
4158                 that = this;
4159 
4160             if (el1.type === Const.OBJECT_TYPE_PLANE3D) {
4161                 if (el2.type === Const.OBJECT_TYPE_PLANE3D) {
4162                     // func = () => view.intersectionPlanePlane(el1, el2)[i];
4163                     func = view.intersectionPlanePlane(el1, el2);
4164                 } else if (el2.type === Const.OBJECT_TYPE_SPHERE3D) {
4165                     func = that.meetPlaneSphere(el1, el2);
4166                 }
4167             } else if (el1.type === Const.OBJECT_TYPE_SPHERE3D) {
4168                 if (el2.type === Const.OBJECT_TYPE_PLANE3D) {
4169                     func = that.meetPlaneSphere(el2, el1);
4170                 } else if (el2.type === Const.OBJECT_TYPE_SPHERE3D) {
4171                     func = that.meetSphereSphere(el1, el2);
4172                 }
4173             }
4174 
4175             return func;
4176         },
4177 
4178         /**
4179          * Intersecting point of three planes in 3D. The planes
4180          * are given in Hesse normal form.
4181          *
4182          * @param {Array} n1 Hesse normal form vector of plane 1
4183          * @param {Number} d1 Hesse normal form right hand side of plane 1
4184          * @param {Array} n2 Hesse normal form vector of plane 2
4185          * @param {Number} d2 Hesse normal form right hand side of plane 2
4186          * @param {Array} n3 Hesse normal form vector of plane 1
4187          * @param {Number} d3 Hesse normal form right hand side of plane 3
4188          * @returns {Array} Coordinates array of length 4 of the intersecting point
4189          */
4190         meet3Planes: function (n1, d1, n2, d2, n3, d3) {
4191             var p = [1, 0, 0, 0],
4192                 n31, n12, n23,
4193                 denom,
4194                 i;
4195 
4196             n31 = Mat.crossProduct(n3.slice(1), n1.slice(1));
4197             n12 = Mat.crossProduct(n1.slice(1), n2.slice(1));
4198             n23 = Mat.crossProduct(n2.slice(1), n3.slice(1));
4199 
4200             denom = Mat.innerProduct(n1.slice(1), n23, 3);
4201             for (i = 0; i < 3; i++) {
4202                 p[i + 1] = (d1 * n23[i] + d2 * n31[i] + d3 * n12[i]) / denom;
4203             }
4204 
4205             return p;
4206         },
4207 
4208         /**
4209          * Direction of intersecting line of two planes in 3D.
4210          *
4211          * @param {Array} v11 First vector spanning plane 1 (homogeneous coordinates)
4212          * @param {Array} v12 Second vector spanning plane 1 (homogeneous coordinates)
4213          * @param {Array} v21 First vector spanning plane 2 (homogeneous coordinates)
4214          * @param {Array} v22 Second vector spanning plane 2 (homogeneous coordinates)
4215          * @returns {Array} Coordinates array of length 4 of the direction  (homogeneous coordinates)
4216          */
4217         meetPlanePlane: function (v11, v12, v21, v22) {
4218             var no1,
4219                 no2,
4220                 v, w;
4221 
4222             v = v11.slice(1);
4223             w = v12.slice(1);
4224             no1 = Mat.crossProduct(v, w);
4225 
4226             v = v21.slice(1);
4227             w = v22.slice(1);
4228             no2 = Mat.crossProduct(v, w);
4229 
4230             w = Mat.crossProduct(no1, no2);
4231             w.unshift(0);
4232             return w;
4233         },
4234 
4235         meetPlaneSphere: function (el1, el2) {
4236             var dis = function () {
4237                     return Mat.innerProduct(el1.normal, el2.center.coords, 4) - el1.d;
4238                 };
4239 
4240             return [
4241                 // Center
4242                 function() {
4243                     return Mat.axpy(-dis(), el1.normal, el2.center.coords);
4244                 },
4245                 // Normal
4246                 el1.normal,
4247                 // Radius
4248                 function () {
4249                     // Radius (returns NaN if spheres don't touch)
4250                     var r = el2.Radius(),
4251                         s = dis();
4252                     return Math.sqrt(r * r - s * s);
4253                 }
4254             ];
4255         },
4256 
4257         meetSphereSphere: function (el1, el2) {
4258             var skew = function () {
4259                     var dist = el1.center.distance(el2.center),
4260                         r1 = el1.Radius(),
4261                         r2 = el2.Radius();
4262                     return (r1 - r2) * (r1 + r2) / (dist * dist);
4263                 };
4264             return [
4265                 // Center
4266                 function () {
4267                     var s = skew();
4268                     return [
4269                         1,
4270                         0.5 * ((1 - s) * el1.center.coords[1] + (1 + s) * el2.center.coords[1]),
4271                         0.5 * ((1 - s) * el1.center.coords[2] + (1 + s) * el2.center.coords[2]),
4272                         0.5 * ((1 - s) * el1.center.coords[3] + (1 + s) * el2.center.coords[3])
4273                     ];
4274                 },
4275                 // Normal
4276                 function() {
4277                     return Stat.subtract(el2.center.coords, el1.center.coords);
4278                 },
4279                 // Radius
4280                 function () {
4281                     // Radius (returns NaN if spheres don't touch)
4282                     var dist = el1.center.distance(el2.center),
4283                         r1 = el1.Radius(),
4284                         r2 = el2.Radius(),
4285                         s = skew(),
4286                         rIxnSq = 0.5 * (r1 * r1 + r2 * r2 - 0.5 * dist * dist * (1 + s * s));
4287                     return Math.sqrt(rIxnSq);
4288                 }
4289             ];
4290         },
4291 
4292         /**
4293          * Test if parameters are inside of allowed ranges
4294          *
4295          * @param {Array} params Array of length 1 or 2
4296          * @param {Array} r_u First range
4297          * @param {Array} [r_v] Second range
4298          * @returns Boolean
4299          * @private
4300          */
4301         _paramsOutOfRange: function(params, r_u, r_v) {
4302             return params[0] < r_u[0] || params[0] > r_u[1] ||
4303                 (params.length > 1 && (params[1] < r_v[0] || params[1] > r_v[1]));
4304         },
4305 
4306         /**
4307          * Given the 2D screen coordinates of a point, finds the nearest point on the given
4308          * parametric curve or surface, and returns its view-space coordinates.
4309          * @param {Array} p 3D coordinates for which the closest point on the curve point is searched.
4310          * @param {JXG.Curve3D|JXG.Surface3D} target Parametric curve or surface to project to.
4311          * @param {Array} params New position of point on the target (i.e. it is a return value),
4312          * modified in place during the search, ending up at the nearest point.
4313          * Usually, point.position is supplied for params.
4314          *
4315          * @returns {Array} Array of length 4 containing the coordinates of the nearest point on the curve or surface.
4316          */
4317         projectCoordsToParametric: function (p, target, n, params) {
4318             // The variables and parameters for the Cobyla constrained
4319             // minimization algorithm are explained in the Cobyla.js comments
4320             var rhobeg,                // initial size of simplex (Cobyla)
4321                 rhoend,                // finial size of simplex (Cobyla)
4322                 iprint = 0,            // no console output (Cobyla)
4323                 maxfun = 200,          // call objective function at most 200 times (Cobyla)
4324                 _minFunc,              // Objective function for Cobyla
4325                 f = Math.random() * 0.01 + 0.5,
4326                 r_u, r_v,
4327                 m = 2 * n;
4328 
4329             // adapt simplex size to parameter range
4330             if (n === 1) {
4331                 r_u = [Type.evaluate(target.range[0]), Type.evaluate(target.range[1])];
4332 
4333                 rhobeg = 0.1 * (r_u[1] - r_u[0]);
4334             } else if (n === 2) {
4335                 r_u = [Type.evaluate(target.range_u[0]), Type.evaluate(target.range_u[1])];
4336                 r_v = [Type.evaluate(target.range_v[0]), Type.evaluate(target.range_v[1])];
4337 
4338                 rhobeg = 0.1 * Math.min(
4339                     r_u[1] - r_u[0],
4340                     r_v[1] - r_v[0]
4341                 );
4342             }
4343             rhoend = rhobeg / 5e6;
4344 
4345             // Minimize distance of the new position to the original position
4346             _minFunc = function (n, m, w, con) {
4347                 var p_new = [
4348                         target.X.apply(target, w),
4349                         target.Y.apply(target, w),
4350                         target.Z.apply(target, w)
4351                     ],
4352                     xDiff = p[0] - p_new[0],
4353                     yDiff = p[1] - p_new[1],
4354                     zDiff = p[2] - p_new[2];
4355 
4356                 if (m >= 2) {
4357                     con[0] =  w[0] - r_u[0];
4358                     con[1] = -w[0] + r_u[1];
4359                 }
4360                 if (m >= 4) {
4361                     con[2] =  w[1] - r_v[0];
4362                     con[3] = -w[1] + r_v[1];
4363                 }
4364 
4365                 return xDiff * xDiff + yDiff * yDiff + zDiff * zDiff;
4366             };
4367 
4368             // First optimization without range constraints to give a smooth draag experience on
4369             // cyclic structures.
4370 
4371             // Set the start values
4372             if (params.length === 0) {
4373                 // If length > 0: take the previous position as start values for the optimization
4374                 params[0] = f * (r_u[0] + r_u[1]);
4375                 if (n === 2) { params[1] = f * (r_v[0] + r_v[1]); }
4376             }
4377             Mat.Nlp.FindMinimum(_minFunc, n, 0, params, rhobeg, rhoend, iprint, maxfun);
4378             // Update p which is used subsequently in _minFunc
4379             p = [target.X.apply(target, params),
4380                 target.Y.apply(target, params),
4381                 target.Z.apply(target, params)
4382             ];
4383 
4384             // If the optimal params are outside of the rang
4385             // Second optimization to obey the range constraints
4386 
4387             if (this._paramsOutOfRange(params, r_u, r_v)) {
4388                 // Set the start values again
4389                 params[0] = f * (r_u[0] + r_u[1]);
4390                 if (n === 2) { params[1] = f * (r_v[0] + r_v[1]); }
4391 
4392                 Mat.Nlp.FindMinimum(_minFunc, n, m, params, rhobeg, rhoend, iprint, maxfun);
4393             }
4394 
4395             return [1,
4396                 target.X.apply(target, params),
4397                 target.Y.apply(target, params),
4398                 target.Z.apply(target, params)
4399             ];
4400         },
4401 
4402         // /**
4403         //  * Given a the screen coordinates of a point, finds the point on the
4404         //  * given parametric curve or surface which is nearest in screen space,
4405         //  * and returns its view-space coordinates.
4406         //  * @param {Array} pScr Screen coordinates to project.
4407         //  * @param {JXG.Curve3D|JXG.Surface3D} target Parametric curve or surface to project to.
4408         //  * @param {Array} params Parameters of point on the target, initially specifying the starting point of
4409         //  * the search. The parameters are modified in place during the search, ending up at the nearest point.
4410         //  * @returns {Array} Array of length 4 containing the coordinates of the nearest point on the curve or surface.
4411         //  */
4412         // projectScreenCoordsToParametric: function (pScr, target, params) {
4413         //     // The variables and parameters for the Cobyla constrained
4414         //     // minimization algorithm are explained in the Cobyla.js comments
4415         //     var rhobeg, // initial size of simplex (Cobyla)
4416         //         rhoend, // finial size of simplex (Cobyla)
4417         //         iprint = 0, // no console output (Cobyla)
4418         //         maxfun = 200, // call objective function at most 200 times (Cobyla)
4419         //         dim = params.length,
4420         //         _minFunc; // objective function (Cobyla)
4421 
4422         //     // adapt simplex size to parameter range
4423         //     if (dim === 1) {
4424         //         rhobeg = 0.1 * (target.range[1] - target.range[0]);
4425         //     } else if (dim === 2) {
4426         //         rhobeg = 0.1 * Math.min(
4427         //             target.range_u[1] - target.range_u[0],
4428         //             target.range_v[1] - target.range_v[0]
4429         //         );
4430         //     }
4431         //     rhoend = rhobeg / 5e6;
4432 
4433         //     // minimize screen distance to cursor
4434         //     _minFunc = function (n, m, w, con) {
4435         //         var c3d = [
4436         //             1,
4437         //             target.X.apply(target, w),
4438         //             target.Y.apply(target, w),
4439         //             target.Z.apply(target, w)
4440         //         ],
4441         //         c2d = target.view.project3DTo2D(c3d),
4442         //         xDiff = pScr[0] - c2d[1],
4443         //         yDiff = pScr[1] - c2d[2];
4444 
4445         //         if (n === 1) {
4446         //             con[0] = w[0] - target.range[0];
4447         //             con[1] = -w[0] + target.range[1];
4448         //         } else if (n === 2) {
4449         //             con[0] = w[0] - target.range_u[0];
4450         //             con[1] = -w[0] + target.range_u[1];
4451         //             con[2] = w[1] - target.range_v[0];
4452         //             con[3] = -w[1] + target.range_v[1];
4453         //         }
4454 
4455         //         return xDiff * xDiff + yDiff * yDiff;
4456         //     };
4457 
4458         //     Mat.Nlp.FindMinimum(_minFunc, dim, 2 * dim, params, rhobeg, rhoend, iprint, maxfun);
4459 
4460         //     return [1, target.X.apply(target, params), target.Y.apply(target, params), target.Z.apply(target, params)];
4461         // },
4462 
4463         project3DTo3DPlane: function (point, normal, foot) {
4464             // TODO: homogeneous 3D coordinates
4465             var sol = [0, 0, 0],
4466                 le,
4467                 d1,
4468                 d2,
4469                 lbda;
4470 
4471             foot = foot || [0, 0, 0];
4472 
4473             le = Mat.norm(normal);
4474             d1 = Mat.innerProduct(point, normal, 3);
4475             d2 = Mat.innerProduct(foot, normal, 3);
4476             // (point - lbda * normal / le) * normal / le == foot * normal / le
4477             // => (point * normal - foot * normal) ==  lbda * le
4478             lbda = (d1 - d2) / le;
4479             sol = Mat.axpy(-lbda, normal, point);
4480 
4481             return sol;
4482         },
4483 
4484         getPlaneBounds: function (v1, v2, q, s, e) {
4485             var s1, s2, e1, e2, mat, rhs, sol;
4486 
4487             if (v1[2] + v2[0] !== 0) {
4488                 mat = [
4489                     [v1[0], v2[0]],
4490                     [v1[1], v2[1]]
4491                 ];
4492                 rhs = [s - q[0], s - q[1]];
4493 
4494                 sol = Numerics.Gauss(mat, rhs);
4495                 s1 = sol[0];
4496                 s2 = sol[1];
4497 
4498                 rhs = [e - q[0], e - q[1]];
4499                 sol = Numerics.Gauss(mat, rhs);
4500                 e1 = sol[0];
4501                 e2 = sol[1];
4502                 return [s1, e1, s2, e2];
4503             }
4504             return null;
4505         },
4506 
4507         /* ***************************************/
4508         /* *** Various ****/
4509         /* ***************************************/
4510 
4511         /**
4512          * Helper function to create curve which displays a Reuleaux polygons.
4513          * @param {Array} points Array of points which should be the vertices of the Reuleaux polygon. Typically,
4514          * these point list is the array vertices of a regular polygon.
4515          * @param {Number} nr Number of vertices
4516          * @returns {Array} An array containing the two functions defining the Reuleaux polygon and the two values
4517          * for the start and the end of the paramtric curve. array may be used as parent array of a
4518          * {@link JXG.Curve}.
4519          *
4520          * @example
4521          * var A = brd.create('point',[-2,-2]);
4522          * var B = brd.create('point',[0,1]);
4523          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0});
4524          * var reuleauxTriangle = brd.create('curve', JXG.Math.Geometry.reuleauxPolygon(pol.vertices, 3),
4525          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
4526          *
4527          * </pre><div class="jxgbox" id="JXG2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div>
4528          * <script type="text/javascript">
4529          * var brd = JXG.JSXGraph.initBoard('JXG2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false});
4530          * var A = brd.create('point',[-2,-2]);
4531          * var B = brd.create('point',[0,1]);
4532          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0});
4533          * var reuleauxTriangle = brd.create('curve', JXG.Math.Geometry.reuleauxPolygon(pol.vertices, 3),
4534          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
4535          * </script><pre>
4536          */
4537         reuleauxPolygon: function (points, nr) {
4538             var beta,
4539                 pi2 = Math.PI * 2,
4540                 pi2_n = pi2 / nr,
4541                 diag = (nr - 1) / 2,
4542                 d = 0,
4543                 makeFct = function (which, trig) {
4544                     return function (t, suspendUpdate) {
4545                         var t1 = ((t % pi2) + pi2) % pi2,
4546                             j = Math.floor(t1 / pi2_n) % nr;
4547 
4548                         if (!suspendUpdate) {
4549                             d = points[0].Dist(points[diag]);
4550                             beta = Mat.Geometry.rad(
4551                                 [points[0].X() + 1, points[0].Y()],
4552                                 points[0],
4553                                 points[diag % nr]
4554                             );
4555                         }
4556 
4557                         if (isNaN(j)) {
4558                             return j;
4559                         }
4560 
4561                         t1 = t1 * 0.5 + j * pi2_n * 0.5 + beta;
4562 
4563                         return points[j][which]() + d * Math[trig](t1);
4564                     };
4565                 };
4566 
4567             return [makeFct("X", 'cos'), makeFct("Y", 'sin'), 0, pi2];
4568         }
4569 
4570     }
4571 );
4572 
4573 export default Mat.Geometry;
4574