1 /*
  2     Copyright 2008-2022
  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 <http://www.gnu.org/licenses/>
 30     and <http://opensource.org/licenses/MIT/>.
 31  */
 32 
 33 
 34 /*global JXG: true, define: true*/
 35 /*jslint nomen: true, plusplus: true*/
 36 
 37 /* depends:
 38  jxg
 39  base/constants
 40  base/coords
 41  math/math
 42  math/numerics
 43  utils/type
 44  */
 45 
 46 /**
 47  * @fileoverview This file contains the Math.Geometry namespace for calculating algebraic/geometric
 48  * stuff like intersection points, angles, midpoint, and so on.
 49  */
 50 
 51 define([
 52     'jxg', 'base/constants', 'base/coords', 'math/math', 'math/numerics', 'utils/type', 'utils/expect'
 53 ], function (JXG, Const, Coords, Mat, Numerics, Type, Expect) {
 54 
 55     "use strict";
 56 
 57     /**
 58      * Math.Geometry namespace definition. This namespace holds geometrical algorithms,
 59      * especially intersection algorithms.
 60      * @name JXG.Math.Geometry
 61      * @namespace
 62      */
 63     Mat.Geometry = {};
 64 
 65 // the splitting is necessary due to the shortcut for the circumcircleMidpoint method to circumcenter.
 66 
 67     JXG.extend(Mat.Geometry, /** @lends JXG.Math.Geometry */ {
 68         /* ***************************************/
 69         /* *** GENERAL GEOMETRIC CALCULATIONS ****/
 70         /* ***************************************/
 71 
 72         /**
 73          * Calculates the angle defined by the points A, B, C.
 74          * @param {JXG.Point,Array} A A point  or [x,y] array.
 75          * @param {JXG.Point,Array} B Another point or [x,y] array.
 76          * @param {JXG.Point,Array} C A circle - no, of course the third point or [x,y] array.
 77          * @deprecated Use {@link JXG.Math.Geometry.rad} instead.
 78          * @see #rad
 79          * @see #trueAngle
 80          * @returns {Number} The angle in radian measure.
 81          */
 82         angle: function (A, B, C) {
 83             var u, v, s, t,
 84                 a = [],
 85                 b = [],
 86                 c = [];
 87 
 88             JXG.deprecated('Geometry.angle()', 'Geometry.rad()');
 89             if (A.coords) {
 90                 a[0] = A.coords.usrCoords[1];
 91                 a[1] = A.coords.usrCoords[2];
 92             } else {
 93                 a[0] = A[0];
 94                 a[1] = A[1];
 95             }
 96 
 97             if (B.coords) {
 98                 b[0] = B.coords.usrCoords[1];
 99                 b[1] = B.coords.usrCoords[2];
100             } else {
101                 b[0] = B[0];
102                 b[1] = B[1];
103             }
104 
105             if (C.coords) {
106                 c[0] = C.coords.usrCoords[1];
107                 c[1] = C.coords.usrCoords[2];
108             } else {
109                 c[0] = C[0];
110                 c[1] = C[1];
111             }
112 
113             u = a[0] - b[0];
114             v = a[1] - b[1];
115             s = c[0] - b[0];
116             t = c[1] - b[1];
117 
118             return Math.atan2(u * t - v * s, u * s + v * t);
119         },
120 
121         /**
122          * Calculates the angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
123          * @param {JXG.Point,Array} A Point or [x,y] array
124          * @param {JXG.Point,Array} B Point or [x,y] array
125          * @param {JXG.Point,Array} C Point or [x,y] array
126          * @see #rad
127          * @returns {Number} The angle in degrees.
128          */
129         trueAngle: function (A, B, C) {
130             return this.rad(A, B, C) * 57.295779513082323; // *180.0/Math.PI;
131         },
132 
133         /**
134          * Calculates the internal angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
135          * @param {JXG.Point,Array} A Point or [x,y] array
136          * @param {JXG.Point,Array} B Point or [x,y] array
137          * @param {JXG.Point,Array} C Point or [x,y] array
138          * @see #trueAngle
139          * @returns {Number} Angle in radians.
140          */
141         rad: function (A, B, C) {
142             var ax, ay, bx, by, cx, cy, phi;
143 
144             if (A.coords) {
145                 ax = A.coords.usrCoords[1];
146                 ay = A.coords.usrCoords[2];
147             } else {
148                 ax = A[0];
149                 ay = A[1];
150             }
151 
152             if (B.coords) {
153                 bx = B.coords.usrCoords[1];
154                 by = B.coords.usrCoords[2];
155             } else {
156                 bx = B[0];
157                 by = B[1];
158             }
159 
160             if (C.coords) {
161                 cx = C.coords.usrCoords[1];
162                 cy = C.coords.usrCoords[2];
163             } else {
164                 cx = C[0];
165                 cy = C[1];
166             }
167 
168             phi = Math.atan2(cy - by, cx - bx) - Math.atan2(ay - by, ax - bx);
169 
170             if (phi < 0) {
171                 phi += 6.2831853071795862;
172             }
173 
174             return phi;
175         },
176 
177         /**
178          * Calculates a point on the bisection line between the three points A, B, C.
179          * As a result, the bisection line is defined by two points:
180          * Parameter B and the point with the coordinates calculated in this function.
181          * Does not work for ideal points.
182          * @param {JXG.Point} A Point
183          * @param {JXG.Point} B Point
184          * @param {JXG.Point} C Point
185          * @param [board=A.board] Reference to the board
186          * @returns {JXG.Coords} Coordinates of the second point defining the bisection.
187          */
188         angleBisector: function (A, B, C, board) {
189             var phiA, phiC, phi,
190                 Ac = A.coords.usrCoords,
191                 Bc = B.coords.usrCoords,
192                 Cc = C.coords.usrCoords,
193                 x, y;
194 
195             if (!Type.exists(board)) {
196                 board = A.board;
197             }
198 
199             // Parallel lines
200             if (Bc[0] === 0) {
201                 return new Coords(Const.COORDS_BY_USER,
202                     [1, (Ac[1] + Cc[1]) * 0.5, (Ac[2] + Cc[2]) * 0.5], board);
203             }
204 
205             // Non-parallel lines
206             x = Ac[1] - Bc[1];
207             y = Ac[2] - Bc[2];
208             phiA =  Math.atan2(y, x);
209 
210             x = Cc[1] - Bc[1];
211             y = Cc[2] - Bc[2];
212             phiC =  Math.atan2(y, x);
213 
214             phi = (phiA + phiC) * 0.5;
215 
216             if (phiA > phiC) {
217                 phi += Math.PI;
218             }
219 
220             x = Math.cos(phi) + Bc[1];
221             y = Math.sin(phi) + Bc[2];
222 
223             return new Coords(Const.COORDS_BY_USER, [1, x, y], board);
224         },
225 
226         // /**
227         //  * Calculates a point on the m-section line between the three points A, B, C.
228         //  * As a result, the m-section line is defined by two points:
229         //  * Parameter B and the point with the coordinates calculated in this function.
230         //  * The m-section generalizes the bisector to any real number.
231         //  * For example, the trisectors of an angle are simply the 1/3-sector and the 2/3-sector.
232         //  * Does not work for ideal points.
233         //  * @param {JXG.Point} A Point
234         //  * @param {JXG.Point} B Point
235         //  * @param {JXG.Point} C Point
236         //  * @param {Number} m Number
237         //  * @param [board=A.board] Reference to the board
238         //  * @returns {JXG.Coords} Coordinates of the second point defining the bisection.
239         //  */
240         // angleMsector: function (A, B, C, m, board) {
241         //     var phiA, phiC, phi,
242         //         Ac = A.coords.usrCoords,
243         //         Bc = B.coords.usrCoords,
244         //         Cc = C.coords.usrCoords,
245         //         x, y;
246 
247         //     if (!Type.exists(board)) {
248         //         board = A.board;
249         //     }
250 
251         //     // Parallel lines
252         //     if (Bc[0] === 0) {
253         //         return new Coords(Const.COORDS_BY_USER,
254         //             [1, (Ac[1] + Cc[1]) * m, (Ac[2] + Cc[2]) * m], board);
255         //     }
256 
257         //     // Non-parallel lines
258         //     x = Ac[1] - Bc[1];
259         //     y = Ac[2] - Bc[2];
260         //     phiA =  Math.atan2(y, x);
261 
262         //     x = Cc[1] - Bc[1];
263         //     y = Cc[2] - Bc[2];
264         //     phiC =  Math.atan2(y, x);
265 
266         //     phi = phiA + ((phiC - phiA) * m);
267 
268         //     if (phiA - phiC > Math.PI) {
269         //         phi += 2*m*Math.PI;
270         //     }
271 
272         //     x = Math.cos(phi) + Bc[1];
273         //     y = Math.sin(phi) + Bc[2];
274 
275         //     return new Coords(Const.COORDS_BY_USER, [1, x, y], board);
276         // },
277 
278         /**
279          * Reflects the point along the line.
280          * @param {JXG.Line} line Axis of reflection.
281          * @param {JXG.Point} point Point to reflect.
282          * @param [board=point.board] Reference to the board
283          * @returns {JXG.Coords} Coordinates of the reflected point.
284          */
285         reflection: function (line, point, board) {
286             // (v,w) defines the slope of the line
287             var x0, y0, x1, y1, v, w, mu,
288                 pc = point.coords.usrCoords,
289                 p1c = line.point1.coords.usrCoords,
290                 p2c = line.point2.coords.usrCoords;
291 
292             if (!Type.exists(board)) {
293                 board = point.board;
294             }
295 
296             v = p2c[1] - p1c[1];
297             w = p2c[2] - p1c[2];
298 
299             x0 = pc[1] - p1c[1];
300             y0 = pc[2] - p1c[2];
301 
302             mu = (v * y0 - w * x0) / (v * v + w * w);
303 
304             // point + mu*(-y,x) is the perpendicular foot
305             x1 = pc[1] + 2 * mu * w;
306             y1 = pc[2] - 2 * mu * v;
307 
308             return new Coords(Const.COORDS_BY_USER, [x1, y1], board);
309         },
310 
311         /**
312          * Computes the new position of a point which is rotated
313          * around a second point (called rotpoint) by the angle phi.
314          * @param {JXG.Point} rotpoint Center of the rotation
315          * @param {JXG.Point} point point to be rotated
316          * @param {Number} phi rotation angle in arc length
317          * @param {JXG.Board} [board=point.board] Reference to the board
318          * @returns {JXG.Coords} Coordinates of the new position.
319          */
320         rotation: function (rotpoint, point, phi, board) {
321             var x0, y0, c, s, x1, y1,
322                 pc = point.coords.usrCoords,
323                 rotpc = rotpoint.coords.usrCoords;
324 
325             if (!Type.exists(board)) {
326                 board = point.board;
327             }
328 
329             x0 = pc[1] - rotpc[1];
330             y0 = pc[2] - rotpc[2];
331 
332             c = Math.cos(phi);
333             s = Math.sin(phi);
334 
335             x1 = x0 * c - y0 * s + rotpc[1];
336             y1 = x0 * s + y0 * c + rotpc[2];
337 
338             return new Coords(Const.COORDS_BY_USER, [x1, y1], board);
339         },
340 
341         /**
342          * Calculates the coordinates of a point on the perpendicular to the given line through
343          * the given point.
344          * @param {JXG.Line} line A line.
345          * @param {JXG.Point} point Point which is projected to the line.
346          * @param {JXG.Board} [board=point.board] Reference to the board
347          * @returns {Array} Array of length two containing coordinates of a point on the perpendicular to the given line
348          *                  through the given point and boolean flag "change".
349          */
350         perpendicular: function (line, point, board) {
351             var x, y, change,
352                 c, z,
353                 A = line.point1.coords.usrCoords,
354                 B = line.point2.coords.usrCoords,
355                 C = point.coords.usrCoords;
356 
357             if (!Type.exists(board)) {
358                 board = point.board;
359             }
360 
361             // special case: point is the first point of the line
362             if (point === line.point1) {
363                 x = A[1] + B[2] - A[2];
364                 y = A[2] - B[1] + A[1];
365                 z = A[0] * B[0];
366 
367                 if (Math.abs(z) < Mat.eps) {
368                     x =  B[2];
369                     y = -B[1];
370                 }
371                 c = [z, x, y];
372                 change = true;
373 
374             // special case: point is the second point of the line
375             } else if (point === line.point2) {
376                 x = B[1] + A[2] - B[2];
377                 y = B[2] - A[1] + B[1];
378                 z = A[0] * B[0];
379 
380                 if (Math.abs(z) < Mat.eps) {
381                     x =  A[2];
382                     y = -A[1];
383                 }
384                 c = [z, x, y];
385                 change = false;
386 
387             // special case: point lies somewhere else on the line
388             } else if (Math.abs(Mat.innerProduct(C, line.stdform, 3)) < Mat.eps) {
389                 x = C[1] + B[2] - C[2];
390                 y = C[2] - B[1] + C[1];
391                 z = B[0];
392 
393                 if (Math.abs(z) < Mat.eps) {
394                     x =  B[2];
395                     y = -B[1];
396                 }
397 
398                 change = true;
399                 if (Math.abs(z) > Mat.eps && Math.abs(x - C[1]) < Mat.eps && Math.abs(y - C[2]) < Mat.eps) {
400                     x = C[1] + A[2] - C[2];
401                     y = C[2] - A[1] + C[1];
402                     change = false;
403                 }
404                 c = [z, x, y];
405 
406             // general case: point does not lie on the line
407             // -> calculate the foot of the dropped perpendicular
408             } else {
409                 c = [0, line.stdform[1], line.stdform[2]];
410                 c = Mat.crossProduct(c, C);                  // perpendicuar to line
411                 c = Mat.crossProduct(c, line.stdform);       // intersection of line and perpendicular
412                 change = true;
413             }
414 
415             return [new Coords(Const.COORDS_BY_USER, c, board), change];
416         },
417 
418         /**
419          * @deprecated Please use {@link JXG.Math.Geometry.circumcenter} instead.
420          */
421         circumcenterMidpoint: function () {
422             JXG.deprecated('Geometry.circumcenterMidpoint()', 'Geometry.circumcenter()');
423             this.circumcenter.apply(this, arguments);
424         },
425 
426         /**
427          * Calculates the center of the circumcircle of the three given points.
428          * @param {JXG.Point} point1 Point
429          * @param {JXG.Point} point2 Point
430          * @param {JXG.Point} point3 Point
431          * @param {JXG.Board} [board=point1.board] Reference to the board
432          * @returns {JXG.Coords} Coordinates of the center of the circumcircle of the given points.
433          */
434         circumcenter: function (point1, point2, point3, board) {
435             var u, v, m1, m2,
436                 A = point1.coords.usrCoords,
437                 B = point2.coords.usrCoords,
438                 C = point3.coords.usrCoords;
439 
440             if (!Type.exists(board)) {
441                 board = point1.board;
442             }
443 
444             u = [B[0] - A[0], -B[2] + A[2], B[1] - A[1]];
445             v = [(A[0] + B[0])  * 0.5, (A[1] + B[1]) * 0.5, (A[2] + B[2]) * 0.5];
446             m1 = Mat.crossProduct(u, v);
447 
448             u = [C[0] - B[0], -C[2] + B[2], C[1] - B[1]];
449             v = [(B[0] + C[0]) * 0.5, (B[1] + C[1]) * 0.5, (B[2] + C[2]) * 0.5];
450             m2 = Mat.crossProduct(u, v);
451 
452             return new Coords(Const.COORDS_BY_USER, Mat.crossProduct(m1, m2), board);
453         },
454 
455         /**
456          * Calculates the Euclidean distance for two given arrays of the same length.
457          * @param {Array} array1 Array of Number
458          * @param {Array} array2 Array of Number
459          * @param {Number} [n] Length of the arrays. Default is the minimum length of the given arrays.
460          * @returns {Number} Euclidean distance of the given vectors.
461          */
462         distance: function (array1, array2, n) {
463             var i,
464                 sum = 0;
465 
466             if (!n) {
467                 n = Math.min(array1.length, array2.length);
468             }
469 
470             for (i = 0; i < n; i++) {
471                 sum += (array1[i] - array2[i]) * (array1[i] - array2[i]);
472             }
473 
474             return Math.sqrt(sum);
475         },
476 
477         /**
478          * Calculates Euclidean distance for two given arrays of the same length.
479          * If one of the arrays contains a zero in the first coordinate, and the Euclidean distance
480          * is different from zero it is a point at infinity and we return Infinity.
481          * @param {Array} array1 Array containing elements of type number.
482          * @param {Array} array2 Array containing elements of type number.
483          * @param {Number} [n] Length of the arrays. Default is the minimum length of the given arrays.
484          * @returns {Number} Euclidean (affine) distance of the given vectors.
485          */
486         affineDistance: function (array1, array2, n) {
487             var d;
488 
489             d = this.distance(array1, array2, n);
490 
491             if (d > Mat.eps && (Math.abs(array1[0]) < Mat.eps || Math.abs(array2[0]) < Mat.eps)) {
492                 return Infinity;
493             }
494 
495             return d;
496         },
497 
498         /**
499          * Affine ratio of three collinear points a, b, c: (c - a) / (b - a).
500          * If r > 1 or r < 0 then c is outside of the segment ab.
501          *
502          * @param {Array|JXG.Coords} a
503          * @param {Array|JXG.Coords} b
504          * @param {Array|JXG.Coords} c
505          * @returns {Number} affine ratio (c - a) / (b - a)
506          */
507         affineRatio: function(a, b, c) {
508             var r = 0.0, dx;
509 
510             if (Type.exists(a.usrCoords)) {
511                 a = a.usrCoords;
512             }
513             if (Type.exists(b.usrCoords)) {
514                 b = b.usrCoords;
515             }
516             if (Type.exists(c.usrCoords)) {
517                 c = c.usrCoords;
518             }
519 
520             dx =  b[1] - a[1];
521 
522             if (Math.abs(dx) > Mat.eps) {
523                 r = (c[1] - a[1]) / dx;
524             } else {
525                 r = (c[2] - a[2]) / (b[2] - a[2]);
526             }
527             return r;
528         },
529 
530         /**
531          * Sort vertices counter clockwise starting with the first point.
532          *
533          * @param {Array} p An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
534          *
535          * @returns {Array}
536          */
537         sortVertices: function (p) {
538             var ll,
539                 ps = Expect.each(p, Expect.coordsArray),
540                 N = ps.length,
541                 lastPoint = null;
542 
543             // If the last point equals the first point, we take the last point out of the array.
544             // It may be that the several points at the end of the array are equal to the first point.
545             // The polygonal chain is been closed by JSXGraph, but this may also have been done by the user.
546             // Therefore, we use a while lopp to pop the last points.
547             while (ps[0][0] === ps[N - 1][0] && ps[0][1] === ps[N - 1][1] && ps[0][2] === ps[N - 1][2]) {
548                 lastPoint = ps.pop();
549                 N--;
550             }
551             // Find the point with the lowest y value
552             // for (i = 1; i < N; i++) {
553             //     if ((ps[i][2] < ps[0][2]) ||
554             //         // if the current and the lowest point have the same y value, pick the one with
555             //         // the lowest x value.
556             //         (Math.abs(ps[i][2] - ps[0][2]) < Mat.eps && ps[i][1] < ps[0][1])) {
557             //         console.log(i, 0);
558             //         ps = Type.swap(ps, i, 0);
559             //     }
560             // }
561 
562             ll = ps[0];
563             // Sort ps in increasing order of the angle between a point and the first point ll.
564             // If a point is equal to the first point ll, the angle is defined to be -Infinity.
565             // Otherwise, atan2 would return zero, which is a value which also attained by points
566             // on the same horizontal line.
567             ps.sort(function (a, b) {
568                 var rad1 = (a[2] === ll[2] && a[1] === ll[1]) ? -Infinity : Math.atan2(a[2] - ll[2], a[1] - ll[1]),
569                     rad2 = (b[2] === ll[2] && b[1] === ll[1]) ? -Infinity : Math.atan2(b[2] - ll[2], b[1] - ll[1]);
570 
571                 return rad1 - rad2;
572             });
573 
574             // If the last point has been taken out of the array, we put it in again.
575             if (lastPoint !== null) {
576                 ps.push(lastPoint);
577             }
578 
579             return ps;
580         },
581 
582         /**
583          * Signed triangle area of the three points given.
584          *
585          * @param {JXG.Point|JXG.Coords|Array} p1
586          * @param {JXG.Point|JXG.Coords|Array} p2
587          * @param {JXG.Point|JXG.Coords|Array} p3
588          *
589          * @returns {Number}
590          */
591         signedTriangle: function (p1, p2, p3) {
592             var A = Expect.coordsArray(p1),
593                 B = Expect.coordsArray(p2),
594                 C = Expect.coordsArray(p3);
595 
596             return 0.5 * ((B[1] - A[1]) * (C[2] - A[2]) - (B[2] - A[2]) * (C[1] - A[1]));
597         },
598 
599         /**
600          * Determine the signed area of a non-selfintersecting polygon.
601          * Surveyor's Formula
602          *
603          * @param {Array} p An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
604          * @param {Boolean} [sort=true]
605          *
606          * @returns {Number}
607          */
608         signedPolygon: function (p, sort) {
609             var i, N,
610                 A = 0,
611                 ps = Expect.each(p, Expect.coordsArray);
612 
613             if (sort === undefined) {
614                 sort = true;
615             }
616 
617             if (!sort) {
618                 ps = this.sortVertices(ps);
619             } else {
620                 // Make sure the polygon is closed. If it is already closed this won't change the sum because the last
621                 // summand will be 0.
622                 ps.unshift(ps[ps.length - 1]);
623             }
624 
625             N = ps.length;
626 
627             for (i = 1; i < N; i++) {
628                 A += ps[i - 1][1] * ps[i][2] - ps[i][1] * ps[i - 1][2];
629             }
630 
631             return 0.5 * A;
632         },
633 
634         /**
635          * Calculate the complex hull of a point cloud.
636          *
637          * @param {Array} points An array containing {@link JXG.Point}, {@link JXG.Coords}, and/or arrays.
638          *
639          * @returns {Array}
640          */
641         GrahamScan: function (points) {
642             var i,
643                 M = 1,
644                 ps = Expect.each(points, Expect.coordsArray),
645                 N = ps.length;
646 
647             ps = this.sortVertices(ps);
648             N = ps.length;
649 
650             for (i = 2; i < N; i++) {
651                 while (this.signedTriangle(ps[M - 1], ps[M], ps[i]) <= 0) {
652                     if (M > 1) {
653                         M -= 1;
654                     } else if (i === N - 1) {
655                         break;
656                     }
657                     i += 1;
658                 }
659 
660                 M += 1;
661                 ps = Type.swap(ps, M, i);
662             }
663 
664             return ps.slice(0, M);
665         },
666 
667         /**
668          * A line can be a segment, a straight, or a ray. So it is not always delimited by point1 and point2
669          * calcStraight determines the visual start point and end point of the line. A segment is only drawn
670          * from start to end point, a straight line is drawn until it meets the boards boundaries.
671          * @param {JXG.Line} el Reference to a line object, that needs calculation of start and end point.
672          * @param {JXG.Coords} point1 Coordinates of the point where line drawing begins. This value is calculated and
673          * set by this method.
674          * @param {JXG.Coords} point2 Coordinates of the point where line drawing ends. This value is calculated and set
675          * by this method.
676          * @param {Number} margin Optional margin, to avoid the display of the small sides of lines.
677          * @returns null
678          * @see Line
679          * @see JXG.Line
680          */
681         calcStraight: function (el, point1, point2, margin) {
682             var takePoint1, takePoint2, intersection, intersect1, intersect2, straightFirst, straightLast,
683                 c, p1, p2;
684 
685             if (!Type.exists(margin)) {
686                 // Enlarge the drawable region slightly. This hides the small sides
687                 // of thick lines in most cases.
688                 margin = 10;
689             }
690 
691             straightFirst = Type.evaluate(el.visProp.straightfirst);
692             straightLast = Type.evaluate(el.visProp.straightlast);
693 
694             // If one of the point is an ideal point in homogeneous coordinates
695             // drawing of line segments or rays are not possible.
696             if (Math.abs(point1.scrCoords[0]) < Mat.eps) {
697                 straightFirst = true;
698             }
699             if (Math.abs(point2.scrCoords[0]) < Mat.eps) {
700                 straightLast = true;
701             }
702 
703             // Do nothing in case of line segments (inside or outside of the board)
704             if (!straightFirst && !straightLast) {
705                 return;
706             }
707 
708             // Compute the stdform of the line in screen coordinates.
709             c = [];
710             c[0] = el.stdform[0] -
711                 el.stdform[1] * el.board.origin.scrCoords[1] / el.board.unitX +
712                 el.stdform[2] * el.board.origin.scrCoords[2] / el.board.unitY;
713             c[1] =  el.stdform[1] / el.board.unitX;
714             c[2] = -el.stdform[2] / el.board.unitY;
715 
716             // p1=p2
717             if (isNaN(c[0] + c[1] + c[2])) {
718                 return;
719             }
720 
721             takePoint1 = false;
722             takePoint2 = false;
723 
724             // Line starts at point1 and point1 is inside the board
725             takePoint1 = !straightFirst &&
726                 Math.abs(point1.usrCoords[0]) >= Mat.eps &&
727                 point1.scrCoords[1] >= 0.0 && point1.scrCoords[1] <= el.board.canvasWidth &&
728                 point1.scrCoords[2] >= 0.0 && point1.scrCoords[2] <= el.board.canvasHeight;
729 
730             // Line ends at point2 and point2 is inside the board
731             takePoint2 = !straightLast &&
732                 Math.abs(point2.usrCoords[0]) >= Mat.eps &&
733                 point2.scrCoords[1] >= 0.0 && point2.scrCoords[1] <= el.board.canvasWidth &&
734                 point2.scrCoords[2] >= 0.0 && point2.scrCoords[2] <= el.board.canvasHeight;
735 
736             // Intersect the line with the four borders of the board.
737             intersection = this.meetLineBoard(c, el.board, margin);
738             intersect1 = intersection[0];
739             intersect2 = intersection[1];
740 
741             /**
742              * At this point we have four points:
743              * point1 and point2 are the first and the second defining point on the line,
744              * intersect1, intersect2 are the intersections of the line with border around the board.
745              */
746 
747             /*
748              * Here we handle rays where both defining points are outside of the board.
749              */
750             // If both points are outside and the complete ray is outside we do nothing
751             if (!takePoint1 && !takePoint2) {
752                 // Ray starting at point 1
753                 if (!straightFirst && straightLast &&
754                         !this.isSameDirection(point1, point2, intersect1) && !this.isSameDirection(point1, point2, intersect2)) {
755                     return;
756                 }
757 
758                 // Ray starting at point 2
759                 if (straightFirst && !straightLast &&
760                         !this.isSameDirection(point2, point1, intersect1) && !this.isSameDirection(point2, point1, intersect2)) {
761                     return;
762                 }
763             }
764 
765             /*
766              * If at least one of the defining points is outside of the board
767              * we take intersect1 or intersect2 as one of the end points
768              * The order is also important for arrows of axes
769              */
770             if (!takePoint1) {
771                 if (!takePoint2) {
772                     // Two border intersection points are used
773                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
774                         p1 = intersect1;
775                         p2 = intersect2;
776                     } else {
777                         p2 = intersect1;
778                         p1 = intersect2;
779                     }
780                 } else {
781                     // One border intersection points is used
782                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
783                         p1 = intersect1;
784                     } else {
785                         p1 = intersect2;
786                     }
787                 }
788             } else {
789                 if (!takePoint2) {
790                     // One border intersection points is used
791                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
792                         p2 = intersect2;
793                     } else {
794                         p2 = intersect1;
795                     }
796                 }
797             }
798 
799             if (p1) {
800                 //point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords.slice(1));
801                 point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords);
802             }
803 
804             if (p2) {
805                 //point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords.slice(1));
806                 point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords);
807             }
808         },
809 
810         /**
811          * A line can be a segment, a straight, or a ray. so it is not always delimited by point1 and point2.
812          *
813          * This method adjusts the line's delimiting points taking into account its nature, the viewport defined
814          * by the board.
815          *
816          * A segment is delimited by start and end point, a straight line or ray is delimited until it meets the
817          * boards boundaries. However, if the line has infinite ticks, it will be delimited by the projection of
818          * the boards vertices onto itself.
819          *
820          * @param {JXG.Line} el Reference to a line object, that needs calculation of start and end point.
821          * @param {JXG.Coords} point1 Coordinates of the point where line drawing begins. This value is calculated and
822          * set by this method.
823          * @param {JXG.Coords} point2 Coordinates of the point where line drawing ends. This value is calculated and set
824          * by this method.
825          * @see Line
826          * @see JXG.Line
827          */
828         calcLineDelimitingPoints: function (el, point1, point2) {
829             var distP1P2, boundingBox, lineSlope,
830                 intersect1, intersect2, straightFirst, straightLast,
831                 c, p1, p2,
832                 takePoint1 = false,
833                 takePoint2 = false;
834 
835             straightFirst = Type.evaluate(el.visProp.straightfirst);
836             straightLast = Type.evaluate(el.visProp.straightlast);
837 
838             // If one of the point is an ideal point in homogeneous coordinates
839             // drawing of line segments or rays are not possible.
840             if (Math.abs(point1.scrCoords[0]) < Mat.eps) {
841                 straightFirst = true;
842             }
843             if (Math.abs(point2.scrCoords[0]) < Mat.eps) {
844                 straightLast = true;
845             }
846 
847             // Compute the stdform of the line in screen coordinates.
848             c = [];
849             c[0] = el.stdform[0] -
850                 el.stdform[1] * el.board.origin.scrCoords[1] / el.board.unitX +
851                 el.stdform[2] * el.board.origin.scrCoords[2] / el.board.unitY;
852             c[1] =  el.stdform[1] / el.board.unitX;
853             c[2] = -el.stdform[2] / el.board.unitY;
854 
855             // p1=p2
856             if (isNaN(c[0] + c[1] + c[2])) {
857                 return;
858             }
859 
860             takePoint1 = !straightFirst;
861             takePoint2 = !straightLast;
862             // Intersect the board vertices on the line to establish the available visual space for the infinite ticks
863             // Based on the slope of the line we can optimise and only project the two outer vertices
864 
865             // boundingBox = [x1, y1, x2, y2] upper left, lower right vertices
866             boundingBox = el.board.getBoundingBox();
867             lineSlope = el.getSlope();
868             if (lineSlope >= 0) {
869                 // project vertices (x2,y1) (x1, y2)
870                 intersect1 = this.projectPointToLine({ coords: { usrCoords: [1, boundingBox[2], boundingBox[1]] } }, el, el.board);
871                 intersect2 = this.projectPointToLine({ coords: { usrCoords: [1, boundingBox[0], boundingBox[3]] } }, el, el.board);
872             } else {
873                 // project vertices (x1, y1) (x2, y2)
874                 intersect1 = this.projectPointToLine({ coords: { usrCoords: [1, boundingBox[0], boundingBox[1]] } }, el, el.board);
875                 intersect2 = this.projectPointToLine({ coords: { usrCoords: [1, boundingBox[2], boundingBox[3]] } }, el, el.board);
876             }
877 
878             /**
879              * we have four points:
880              * point1 and point2 are the first and the second defining point on the line,
881              * intersect1, intersect2 are the intersections of the line with border around the board.
882              */
883 
884             /*
885              * Here we handle rays/segments where both defining points are outside of the board.
886              */
887             if (!takePoint1 && !takePoint2) {
888                 // Segment, if segment does not cross the board, do nothing
889                 if (!straightFirst && !straightLast) {
890                     distP1P2 = point1.distance(Const.COORDS_BY_USER, point2);
891                     // if  intersect1 not between point1 and point2
892                     if (Math.abs(point1.distance(Const.COORDS_BY_USER, intersect1) +
893                             intersect1.distance(Const.COORDS_BY_USER, point2) - distP1P2) > Mat.eps) {
894                         return;
895                     }
896                     // if insersect2 not between point1 and point2
897                     if (Math.abs(point1.distance(Const.COORDS_BY_USER, intersect2) +
898                             intersect2.distance(Const.COORDS_BY_USER, point2) - distP1P2) > Mat.eps) {
899                         return;
900                     }
901                 }
902 
903                 // If both points are outside and the complete ray is outside we do nothing
904                 // Ray starting at point 1
905                 if (!straightFirst && straightLast &&
906                         !this.isSameDirection(point1, point2, intersect1) && !this.isSameDirection(point1, point2, intersect2)) {
907                     return;
908                 }
909 
910                 // Ray starting at point 2
911                 if (straightFirst && !straightLast &&
912                         !this.isSameDirection(point2, point1, intersect1) && !this.isSameDirection(point2, point1, intersect2)) {
913                     return;
914                 }
915             }
916 
917             /*
918              * If at least one of the defining points is outside of the board
919              * we take intersect1 or intersect2 as one of the end points
920              * The order is also important for arrows of axes
921              */
922             if (!takePoint1) {
923                 if (!takePoint2) {
924                     // Two border intersection points are used
925                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
926                         p1 = intersect1;
927                         p2 = intersect2;
928                     } else {
929                         p2 = intersect1;
930                         p1 = intersect2;
931                     }
932                 } else {
933                     // One border intersection points is used
934                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
935                         p1 = intersect1;
936                     } else {
937                         p1 = intersect2;
938                     }
939                 }
940             } else {
941                 if (!takePoint2) {
942                     // One border intersection points is used
943                     if (this.isSameDir(point1, point2, intersect1, intersect2)) {
944                         p2 = intersect2;
945                     } else {
946                         p2 = intersect1;
947                     }
948                 }
949             }
950 
951             if (p1) {
952                 //point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords.slice(1));
953                 point1.setCoordinates(Const.COORDS_BY_USER, p1.usrCoords);
954             }
955 
956             if (p2) {
957                 //point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords.slice(1));
958                 point2.setCoordinates(Const.COORDS_BY_USER, p2.usrCoords);
959             }
960         },
961 
962         /**
963          * Calculates the visProp.position corresponding to a given angle.
964          * @param {number} angle angle in radians. Must be in range (-2pi,2pi).
965          */
966         calcLabelQuadrant: function(angle) {
967             var q;
968             if (angle < 0) {
969                 angle += 2*Math.PI;
970             }
971             q = Math.floor((angle+Math.PI/8)/(Math.PI/4))%8;
972             return ['rt','urt','top','ulft','lft','llft','lrt'][q];
973         },
974 
975         /**
976          * The vectors <tt>p2-p1</tt> and <tt>i2-i1</tt> are supposed to be collinear. If their cosine is positive
977          * they point into the same direction otherwise they point in opposite direction.
978          * @param {JXG.Coords} p1
979          * @param {JXG.Coords} p2
980          * @param {JXG.Coords} i1
981          * @param {JXG.Coords} i2
982          * @returns {Boolean} True, if <tt>p2-p1</tt> and <tt>i2-i1</tt> point into the same direction
983          */
984         isSameDir: function (p1, p2, i1, i2) {
985             var dpx = p2.usrCoords[1] - p1.usrCoords[1],
986                 dpy = p2.usrCoords[2] - p1.usrCoords[2],
987                 dix = i2.usrCoords[1] - i1.usrCoords[1],
988                 diy = i2.usrCoords[2] - i1.usrCoords[2];
989 
990             if (Math.abs(p2.usrCoords[0]) < Mat.eps) {
991                 dpx = p2.usrCoords[1];
992                 dpy = p2.usrCoords[2];
993             }
994 
995             if (Math.abs(p1.usrCoords[0]) < Mat.eps) {
996                 dpx = -p1.usrCoords[1];
997                 dpy = -p1.usrCoords[2];
998             }
999 
1000             return dpx * dix + dpy * diy >= 0;
1001         },
1002 
1003         /**
1004          * If you're looking from point "start" towards point "s" and you can see the point "p", return true.
1005          * Otherwise return false.
1006          * @param {JXG.Coords} start The point you're standing on.
1007          * @param {JXG.Coords} p The point in which direction you're looking.
1008          * @param {JXG.Coords} s The point that should be visible.
1009          * @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.
1010          */
1011         isSameDirection: function (start, p, s) {
1012             var dx, dy, sx, sy, r = false;
1013 
1014             dx = p.usrCoords[1] - start.usrCoords[1];
1015             dy = p.usrCoords[2] - start.usrCoords[2];
1016 
1017             sx = s.usrCoords[1] - start.usrCoords[1];
1018             sy = s.usrCoords[2] - start.usrCoords[2];
1019 
1020             if (Math.abs(dx) < Mat.eps) {
1021                 dx = 0;
1022             }
1023 
1024             if (Math.abs(dy) < Mat.eps) {
1025                 dy = 0;
1026             }
1027 
1028             if (Math.abs(sx) < Mat.eps) {
1029                 sx = 0;
1030             }
1031 
1032             if (Math.abs(sy) < Mat.eps) {
1033                 sy = 0;
1034             }
1035 
1036             if (dx >= 0 && sx >= 0) {
1037                 r = (dy >= 0 && sy >= 0) || (dy <= 0 && sy <= 0);
1038             } else if (dx <= 0 && sx <= 0) {
1039                 r = (dy >= 0 && sy >= 0) || (dy <= 0 && sy <= 0);
1040             }
1041 
1042             return r;
1043         },
1044 
1045         /**
1046          * Determinant of three points in the Euclidean plane.
1047          * Zero, if the points are collinear. Used to determine of a point q is left or
1048          * right to a segment defined by points p1 and p2.
1049          *
1050          * @param  {Array} p1 Coordinates of the first point of the segment. Array of length 3. First coordinate is equal to 1.
1051          * @param  {Array} p2 Coordinates of the second point of the segment. Array of length 3. First coordinate is equal to 1.
1052          * @param  {Array} q Coordinates of the point. Array of length 3. First coordinate is equal to 1.
1053          * @return {Number} Signed area of the triangle formed by these three points.
1054          *
1055          * @see #windingNumber
1056          */
1057         det3p: function(p1, p2, q) {
1058             return (p1[1] - q[1]) * (p2[2] - q[2]) - (p2[1] - q[1]) * (p1[2] - q[2]);
1059         },
1060 
1061         /**
1062          * Winding number of a point in respect to a polygon path.
1063          *
1064          * The point is regarded outside if the winding number is zero,
1065          * inside otherwise. The algorithm tries to find degenerate cases, i.e.
1066          * if the point is on the path. This is regarded as "outside".
1067          * If the point is a vertex of the path, it is regarded as "inside".
1068          *
1069          * Implementation of algorithm 7 from "The point in polygon problem for
1070          * arbitrary polygons" by Kai Hormann and Alexander Agathos, Computational Geometry,
1071          * Volume 20, Issue 3, November 2001, Pages 131-144.
1072          *
1073          * @param  {Array} usrCoords Homogenous coordinates of the point
1074          * @param  {Array} path      Array of points / coords determining a path, i.e. the vertices of the polygon / path. The array elements
1075          * do not have to be full points, but have to have a subobject "coords" or should be of type JXG.Coords.
1076          * @param  {Boolean} [doNotClosePath=false] If true the last point of the path is not connected to the first point.
1077          * This is necessary if the path consists of two or more closed subpaths, e.g. if the figure has a hole.
1078          *
1079          * @return {Number}          Winding number of the point. The point is
1080          *                           regarded outside if the winding number is zero,
1081          *                           inside otherwise.
1082          */
1083         windingNumber: function(usrCoords, path, doNotClosePath) {
1084             var wn = 0,
1085                 le = path.length,
1086                 x = usrCoords[1],
1087                 y = usrCoords[2],
1088                 p0, p1, p2, d, sign, i, off = 0;
1089 
1090             if (le === 0) {
1091                 return 0;
1092             }
1093 
1094             doNotClosePath = doNotClosePath || false;
1095             if (doNotClosePath) {
1096                 off = 1;
1097             }
1098 
1099             // Infinite points are declared outside
1100             if (isNaN(x) || isNaN(y)) {
1101                 return 1;
1102             }
1103 
1104             if (Type.exists(path[0].coords)) {
1105                 p0 = path[0].coords;
1106                 p1 = path[le - 1].coords;
1107             } else {
1108                 p0 = path[0];
1109                 p1 = path[le - 1];
1110             }
1111             // Handle the case if the point is the first vertex of the path, i.e. inside.
1112             if (p0.usrCoords[1] === x && p0.usrCoords[2] === y) {
1113                 return 1;
1114             }
1115 
1116             for (i = 0; i < le - off; i++) {
1117                 // Consider the edge from p1 = path[i] to p2 = path[i+1]isClosedPath
1118                 if (Type.exists(path[i].coords)) {
1119                     p1 = path[i].coords.usrCoords;
1120                     p2 = path[(i + 1) % le].coords.usrCoords;
1121                 } else {
1122                     p1 = path[i].usrCoords;
1123                     p2 = path[(i + 1) % le].usrCoords;
1124                 }
1125 
1126                 // If one of the two points p1, p2 is undefined or infinite,
1127                 // move on.
1128                 if (p1[0] === 0 || p2[0] === 0 ||
1129                     isNaN(p1[1]) || isNaN(p2[1]) ||
1130                     isNaN(p1[2]) || isNaN(p2[2])) {
1131                     continue;
1132                 }
1133 
1134                 if (p2[2] === y) {
1135                     if (p2[1] === x) {
1136                         return 1;
1137                     }
1138                     if (p1[2] === y && ((p2[1] > x) === (p1[1] < x))) {
1139                         return 0;
1140                     }
1141                 }
1142 
1143                 if ((p1[2] < y) !== (p2[2] < y)) {                        // Crossing
1144                     sign = 2 * ((p2[2] > p1[2]) ? 1 : 0) - 1;
1145                     if (p1[1] >= x) {
1146                         if (p2[1] > x) {
1147                             wn += sign;
1148                         } else {
1149                             d = this.det3p(p1, p2, usrCoords);
1150                             if (d === 0) {  // Point is on line, i.e. outside
1151                                 return 0;
1152                             }
1153                             if ((d > 0 + Mat.eps) === (p2[2] > p1[2])) {  // Right crossing
1154                                 wn += sign;
1155                             }
1156                         }
1157                     } else {
1158                         if (p2[1] > x) {
1159                             d = this.det3p(p1, p2, usrCoords);
1160                             if ((d > 0 + Mat.eps) === (p2[2] > p1[2])) {  // Right crossing
1161                                 wn += sign;
1162                             }
1163                         }
1164                     }
1165                 }
1166             }
1167 
1168             return wn;
1169         },
1170 
1171         /**
1172          * Decides if a point (x,y) is inside of a path / polygon.
1173          * Does not work correct if the path has hole. In this case, windingNumber is the preferred method.
1174          * Implements W. Randolf Franklin's pnpoly method.
1175          *
1176          * See <a href="https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html">https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html</a>.
1177          *
1178          * @param {Number} x_in x-coordinate (screen or user coordinates)
1179          * @param {Number} y_in y-coordinate (screen or user coordinates)
1180          * @param  {Array} path  Array of points / coords determining a path, i.e. the vertices of the polygon / path. The array elements
1181          * do not have to be full points, but have to have a subobject "coords" or should be of type JXG.Coords.
1182          * @param {Number} [coord_type=JXG.COORDS_BY_SCREEN] Type of coordinates used here.
1183          *   Possible values are <b>JXG.COORDS_BY_USER</b> and <b>JXG.COORDS_BY_SCREEN</b>.
1184          *   Default value is JXG.COORDS_BY_SCREEN.
1185          *
1186          * @returns {Boolean} if (x_in, y_in) is inside of the polygon.
1187          * @see JXG.Polygon.hasPoint
1188          * @see JXG.Polygon.pnpoly
1189          * @see #windingNumber
1190          *
1191          * @example
1192          * var pol = board.create('polygon', [[-1,2], [2,2], [-1,4]]);
1193          * var p = board.create('point', [4, 3]);
1194          * var txt = board.create('text', [-1, 0.5, function() {
1195          *   return 'Point A is inside of the polygon = ' +
1196          *     JXG.Math.Geometry.pnpoly(p.X(), p.Y(), JXG.COORDS_BY_USER, pol.vertices);
1197          * }]);
1198          *
1199          * </pre><div id="JXG4656ed42-f965-4e35-bb66-c334a4529683" class="jxgbox" style="width: 300px; height: 300px;"></div>
1200          * <script type="text/javascript">
1201          *     (function() {
1202          *         var board = JXG.JSXGraph.initBoard('JXG4656ed42-f965-4e35-bb66-c334a4529683',
1203          *             {boundingbox: [-2, 5, 5,-2], axis: true, showcopyright: false, shownavigation: false});
1204          *     var pol = board.create('polygon', [[-1,2], [2,2], [-1,4]]);
1205          *     var p = board.create('point', [4, 3]);
1206          *     var txt = board.create('text', [-1, 0.5, function() {
1207          *     		return 'Point A is inside of the polygon = ' + JXG.Math.Geometry.pnpoly(p.X(), p.Y(), JXG.COORDS_BY_USER, pol.vertices);
1208          *     }]);
1209          *
1210          *     })();
1211          *
1212          * </script><pre>
1213          *
1214          */
1215         pnpoly: function(x_in, y_in, path, coord_type) {
1216             var i, j, len,
1217                 x, y, crds,
1218                 v = path,
1219                 vi, vj,
1220                 isIn = false;
1221 
1222             if (coord_type === Const.COORDS_BY_USER) {
1223                 crds = new Coords(Const.COORDS_BY_USER, [x_in, y_in], this.board);
1224                 x = crds.scrCoords[1];
1225                 y = crds.scrCoords[2];
1226             } else {
1227                 x = x_in;
1228                 y = y_in;
1229             }
1230 
1231             len = path.length;
1232             for (i = 0, j = len - 2; i < len - 1; j = i++) {
1233                 vi = (Type.exists(v[i].coords)) ? v[i].coords : v[i];
1234                 vj = (Type.exists(v[j].coords)) ? v[j].coords : v[j];
1235 
1236                 if (((vi.scrCoords[2] > y) !== (vj.scrCoords[2] > y)) &&
1237                     (x < (vj.scrCoords[1] - vi.scrCoords[1]) *
1238                     (y - vi.scrCoords[2]) / (vj.scrCoords[2] - vi.scrCoords[2]) + vi.scrCoords[1])
1239                    ) {
1240                     isIn = !isIn;
1241                 }
1242             }
1243 
1244             return isIn;
1245         },
1246 
1247         /****************************************/
1248         /****          INTERSECTIONS         ****/
1249         /****************************************/
1250 
1251         /**
1252          * Generate the function which computes the coordinates of the intersection point.
1253          * Primarily used in {@link JXG.Point#createIntersectionPoint}.
1254          * @param {JXG.Board} board object
1255          * @param {JXG.Line,JXG.Circle_JXG.Line,JXG.Circle_Number} el1,el2,i The result will be a intersection point on el1 and el2.
1256          * i determines the intersection point if two points are available: <ul>
1257          *   <li>i==0: use the positive square root,</li>
1258          *   <li>i==1: use the negative square root.</li></ul>
1259          * See further {@link JXG.Point#createIntersectionPoint}.
1260          * @param {Boolean} alwaysintersect. Flag that determines if segments and arc can have an outer intersection point
1261          * on their defining line or circle.
1262          * @returns {Function} Function returning a {@link JXG.Coords} object that determines
1263          * the intersection point.
1264          */
1265         intersectionFunction: function (board, el1, el2, i, j, alwaysintersect) {
1266             var func, that = this,
1267                 el1_isArcType = false,
1268                 el2_isArcType = false;
1269 
1270             el1_isArcType = (el1.elementClass === Const.OBJECT_CLASS_CURVE &&
1271                 (el1.type === Const.OBJECT_TYPE_ARC || el1.type === Const.OBJECT_TYPE_SECTOR)
1272                 ) ? true : false;
1273             el2_isArcType = (el2.elementClass === Const.OBJECT_CLASS_CURVE &&
1274                 (el2.type === Const.OBJECT_TYPE_ARC || el2.type === Const.OBJECT_TYPE_SECTOR)
1275                 ) ? true : false;
1276 
1277             if ((el1.elementClass === Const.OBJECT_CLASS_CURVE || el2.elementClass === Const.OBJECT_CLASS_CURVE) &&
1278                 (el1.elementClass === Const.OBJECT_CLASS_CURVE || el1.elementClass === Const.OBJECT_CLASS_CIRCLE) &&
1279                 (el2.elementClass === Const.OBJECT_CLASS_CURVE || el2.elementClass === Const.OBJECT_CLASS_CIRCLE) /*&&
1280                 !(el1_isArcType && el2_isArcType)*/ ) {
1281                 // curve - curve
1282                 // with the exception that both elements are arc types
1283                 /** @ignore */
1284                 func = function () {
1285                     return that.meetCurveCurve(el1, el2, i, j, el1.board);
1286                 };
1287 
1288             } else if ((
1289                         el1.elementClass === Const.OBJECT_CLASS_CURVE &&
1290                         !el1_isArcType &&
1291                         el2.elementClass === Const.OBJECT_CLASS_LINE
1292                        ) ||
1293                        (
1294                         el2.elementClass === Const.OBJECT_CLASS_CURVE &&
1295                         !el2_isArcType &&
1296                         el1.elementClass === Const.OBJECT_CLASS_LINE
1297                        )
1298                     ) {
1299                 // curve - line (this includes intersections between conic sections and lines)
1300                 // with the exception that the curve is of arc type
1301                 /** @ignore */
1302                 func = function () {
1303                     return that.meetCurveLine(el1, el2, i, el1.board, alwaysintersect);
1304                 };
1305 
1306             } else if (el1.type === Const.OBJECT_TYPE_POLYGON || el2.type === Const.OBJECT_TYPE_POLYGON) {
1307                 // polygon - other
1308                 // Uses the Greiner-Hormann clipping algorithm
1309                 // Not implemented: polygon - point
1310 
1311                 if (el1.elementClass === Const.OBJECT_CLASS_LINE) {
1312                     // line - path
1313                     /** @ignore */
1314                     func = function () {
1315                         return that.meetPolygonLine(el2, el1, i, el1.board, alwaysintersect);
1316                     };
1317                 } else if (el2.elementClass === Const.OBJECT_CLASS_LINE) {
1318                     // path - line
1319                     func = function () {
1320                         return that.meetPolygonLine(el1, el2, i, el1.board, alwaysintersect);
1321                     };
1322                 } else {
1323                     // path - path
1324                     /** @ignore */
1325                     func = function () {
1326                         return that.meetPathPath(el1, el2, i, el1.board);
1327                     };
1328                 }
1329 
1330             } else if (el1.elementClass === Const.OBJECT_CLASS_LINE && el2.elementClass === Const.OBJECT_CLASS_LINE) {
1331                 // line - line, lines may also be segments.
1332                 /** @ignore */
1333                 func = function () {
1334                     var res, c,
1335                         first1 = Type.evaluate(el1.visProp.straightfirst),
1336                         last1 = Type.evaluate(el1.visProp.straightlast),
1337                         first2 = Type.evaluate(el2.visProp.straightfirst),
1338                         last2 = Type.evaluate(el2.visProp.straightlast);
1339 
1340                     /**
1341                      * If one of the lines is a segment or ray and
1342                      * the intersection point should disappear if outside
1343                      * of the segment or ray we call
1344                      * meetSegmentSegment
1345                      */
1346                     if (!Type.evaluate(alwaysintersect) && (!first1 || !last1 || !first2 || !last2)) {
1347                         res = that.meetSegmentSegment(
1348                             el1.point1.coords.usrCoords,
1349                             el1.point2.coords.usrCoords,
1350                             el2.point1.coords.usrCoords,
1351                             el2.point2.coords.usrCoords
1352                         );
1353 
1354                         if ((!first1 && res[1] < 0) || (!last1 && res[1] > 1) ||
1355                                 (!first2 && res[2] < 0) || (!last2 && res[2] > 1)) {
1356                             // Non-existent
1357                             c = [0, NaN, NaN];
1358                         } else {
1359                             c = res[0];
1360                         }
1361 
1362                         return (new Coords(Const.COORDS_BY_USER, c, el1.board));
1363                     }
1364 
1365                     return that.meet(el1.stdform, el2.stdform, i, el1.board);
1366                 };
1367             } else {
1368                 // All other combinations of circles and lines,
1369                 // Arc types are treated as circles.
1370                 /** @ignore */
1371                 func = function () {
1372                     var res = that.meet(el1.stdform, el2.stdform, i, el1.board),
1373                         has = true,
1374                         first, last, r, dx;
1375 
1376                     if (alwaysintersect) {
1377                         return res;
1378                     }
1379                     if (el1.elementClass === Const.OBJECT_CLASS_LINE) {
1380                         first = Type.evaluate(el1.visProp.straightfirst);
1381                         last  = Type.evaluate(el1.visProp.straightlast);
1382                         if (!first || !last) {
1383                             r = that.affineRatio(el1.point1.coords, el1.point2.coords, res);
1384                             if ( (!last && r > 1 + Mat.eps) || (!first && r < 0 - Mat.eps) ) {
1385                                 return (new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board));
1386                             }
1387                         }
1388                     }
1389                     if (el2.elementClass === Const.OBJECT_CLASS_LINE) {
1390                         first = Type.evaluate(el2.visProp.straightfirst);
1391                         last  = Type.evaluate(el2.visProp.straightlast);
1392                         if (!first || !last) {
1393                             r = that.affineRatio(el2.point1.coords, el2.point2.coords, res);
1394                             if ( (!last && r > 1 + Mat.eps) || (!first && r < 0 - Mat.eps) ) {
1395                                 return (new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board));
1396                             }
1397                         }
1398                     }
1399                     if (el1_isArcType) {
1400                         has = that.coordsOnArc(el1, res);
1401                         if (has && el2_isArcType) {
1402                             has = that.coordsOnArc(el2, res);
1403                         }
1404                         if (!has) {
1405                             return (new Coords(JXG.COORDS_BY_USER, [0, NaN, NaN], el1.board));
1406                         }
1407                     }
1408                     return res;
1409                 };
1410             }
1411 
1412             return func;
1413         },
1414 
1415         /**
1416          * Returns true if the coordinates are on the arc element,
1417          * false otherwise. Usually, coords is an intersection
1418          * on the circle line. Now it is decided if coords are on the
1419          * circle restricted to the arc line.
1420          * @param  {Arc} arc arc or sector element
1421          * @param  {JXG.Coords} coords Coords object of an intersection
1422          * @returns {Boolean}
1423          * @private
1424          */
1425         coordsOnArc: function(arc, coords) {
1426             var angle = this.rad(arc.radiuspoint, arc.center, coords.usrCoords.slice(1)),
1427                 alpha = 0.0,
1428                 beta = this.rad(arc.radiuspoint, arc.center, arc.anglepoint),
1429                 ev_s = Type.evaluate(arc.visProp.selection);
1430 
1431             if ((ev_s === 'minor' && beta > Math.PI) ||
1432                 (ev_s === 'major' && beta < Math.PI)) {
1433                 alpha = beta;
1434                 beta = 2 * Math.PI;
1435             }
1436             if (angle < alpha || angle > beta) {
1437                 return false;
1438             }
1439             return true;
1440         },
1441 
1442         /**
1443          * Computes the intersection of a pair of lines, circles or both.
1444          * It uses the internal data array stdform of these elements.
1445          * @param {Array} el1 stdform of the first element (line or circle)
1446          * @param {Array} el2 stdform of the second element (line or circle)
1447          * @param {Number} i Index of the intersection point that should be returned.
1448          * @param board Reference to the board.
1449          * @returns {JXG.Coords} Coordinates of one of the possible two or more intersection points.
1450          * Which point will be returned is determined by i.
1451          */
1452         meet: function (el1, el2, i, board) {
1453             var result,
1454                 eps = Mat.eps;
1455 
1456             // line line
1457             if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) < eps) {
1458                 result = this.meetLineLine(el1, el2, i, board);
1459             // circle line
1460             } else if (Math.abs(el1[3]) >= eps && Math.abs(el2[3]) < eps) {
1461                 result = this.meetLineCircle(el2, el1, i, board);
1462             // line circle
1463             } else if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) >= eps) {
1464                 result = this.meetLineCircle(el1, el2, i, board);
1465             // circle circle
1466             } else {
1467                 result = this.meetCircleCircle(el1, el2, i, board);
1468             }
1469 
1470             return result;
1471         },
1472 
1473         /**
1474          * Intersection of the line with the board
1475          * @param  {Array}     line   stdform of the line in screen coordinates
1476          * @param  {JXG.Board} board  reference to a board.
1477          * @param  {Number}    margin optional margin, to avoid the display of the small sides of lines.
1478          * @returns {Array}            [intersection coords 1, intersection coords 2]
1479          */
1480         meetLineBoard: function (line, board, margin) {
1481              // Intersect the line with the four borders of the board.
1482             var s = [], intersect1, intersect2, i, j;
1483 
1484             if (!Type.exists(margin)) {
1485                 margin = 0;
1486             }
1487 
1488             // top
1489             s[0] = Mat.crossProduct(line, [margin, 0, 1]);
1490             // left
1491             s[1] = Mat.crossProduct(line, [margin, 1, 0]);
1492             // bottom
1493             s[2] = Mat.crossProduct(line, [-margin - board.canvasHeight, 0, 1]);
1494             // right
1495             s[3] = Mat.crossProduct(line, [-margin - board.canvasWidth, 1, 0]);
1496 
1497             // Normalize the intersections
1498             for (i = 0; i < 4; i++) {
1499                 if (Math.abs(s[i][0]) > Mat.eps) {
1500                     for (j = 2; j > 0; j--) {
1501                         s[i][j] /= s[i][0];
1502                     }
1503                     s[i][0] = 1.0;
1504                 }
1505             }
1506 
1507             // line is parallel to "left", take "top" and "bottom"
1508             if (Math.abs(s[1][0]) < Mat.eps) {
1509                 intersect1 = s[0];                          // top
1510                 intersect2 = s[2];                          // bottom
1511             // line is parallel to "top", take "left" and "right"
1512             } else if (Math.abs(s[0][0]) < Mat.eps) {
1513                 intersect1 = s[1];                          // left
1514                 intersect2 = s[3];                          // right
1515             // left intersection out of board (above)
1516             } else if (s[1][2] < 0) {
1517                 intersect1 = s[0];                          // top
1518 
1519                 // right intersection out of board (below)
1520                 if (s[3][2] > board.canvasHeight) {
1521                     intersect2 = s[2];                      // bottom
1522                 } else {
1523                     intersect2 = s[3];                      // right
1524                 }
1525             // left intersection out of board (below)
1526             } else if (s[1][2] > board.canvasHeight) {
1527                 intersect1 = s[2];                          // bottom
1528 
1529                 // right intersection out of board (above)
1530                 if (s[3][2] < 0) {
1531                     intersect2 = s[0];                      // top
1532                 } else {
1533                     intersect2 = s[3];                      // right
1534                 }
1535             } else {
1536                 intersect1 = s[1];                          // left
1537 
1538                 // right intersection out of board (above)
1539                 if (s[3][2] < 0) {
1540                     intersect2 = s[0];                      // top
1541                 // right intersection out of board (below)
1542                 } else if (s[3][2] > board.canvasHeight) {
1543                     intersect2 = s[2];                      // bottom
1544                 } else {
1545                     intersect2 = s[3];                      // right
1546                 }
1547             }
1548 
1549             intersect1 = new Coords(Const.COORDS_BY_SCREEN, intersect1.slice(1), board);
1550             intersect2 = new Coords(Const.COORDS_BY_SCREEN, intersect2.slice(1), board);
1551             return [intersect1, intersect2];
1552         },
1553 
1554         /**
1555          * Intersection of two lines.
1556          * @param {Array} l1 stdform of the first line
1557          * @param {Array} l2 stdform of the second line
1558          * @param {number} i unused
1559          * @param {JXG.Board} board Reference to the board.
1560          * @returns {JXG.Coords} Coordinates of the intersection point.
1561          */
1562         meetLineLine: function (l1, l2, i, board) {
1563             /*
1564             var s = Mat.crossProduct(l1, l2);
1565 
1566             if (Math.abs(s[0]) > Mat.eps) {
1567                 s[1] /= s[0];
1568                 s[2] /= s[0];
1569                 s[0] = 1.0;
1570             }
1571             */
1572             var s = isNaN(l1[5] + l2[5]) ? [0, 0, 0] : Mat.crossProduct(l1, l2);
1573             return new Coords(Const.COORDS_BY_USER, s, board);
1574         },
1575 
1576         /**
1577          * Intersection of line and circle.
1578          * @param {Array} lin stdform of the line
1579          * @param {Array} circ stdform of the circle
1580          * @param {number} i number of the returned intersection point.
1581          *   i==0: use the positive square root,
1582          *   i==1: use the negative square root.
1583          * @param {JXG.Board} board Reference to a board.
1584          * @returns {JXG.Coords} Coordinates of the intersection point
1585          */
1586         meetLineCircle: function (lin, circ, i, board) {
1587             var a, b, c, d, n,
1588                 A, B, C, k, t;
1589 
1590             // Radius is zero, return center of circle
1591             if (circ[4] < Mat.eps) {
1592                 if (Math.abs(Mat.innerProduct([1, circ[6], circ[7]], lin, 3)) < Mat.eps) {
1593                     return new Coords(Const.COORDS_BY_USER, circ.slice(6, 8), board);
1594                 }
1595 
1596                 return new Coords(Const.COORDS_BY_USER, [NaN, NaN], board);
1597             }
1598             c = circ[0];
1599             b = circ.slice(1, 3);
1600             a = circ[3];
1601             d = lin[0];
1602             n = lin.slice(1, 3);
1603 
1604             // Line is assumed to be normalized. Therefore, nn==1 and we can skip some operations:
1605             /*
1606              var nn = n[0]*n[0]+n[1]*n[1];
1607              A = a*nn;
1608              B = (b[0]*n[1]-b[1]*n[0])*nn;
1609              C = a*d*d - (b[0]*n[0]+b[1]*n[1])*d + c*nn;
1610              */
1611             A = a;
1612             B = (b[0] * n[1] - b[1] * n[0]);
1613             C = a * d * d - (b[0] * n[0] + b[1] * n[1]) * d + c;
1614 
1615             k = B * B - 4 * A * C;
1616             if (k > -Mat.eps * Mat.eps) {
1617                 k = Math.sqrt(Math.abs(k));
1618                 t = [(-B + k) / (2 * A), (-B - k) / (2 * A)];
1619 
1620                 return ((i === 0) ?
1621                         new Coords(Const.COORDS_BY_USER, [-t[0] * (-n[1]) - d * n[0], -t[0] * n[0] - d * n[1]], board) :
1622                         new Coords(Const.COORDS_BY_USER, [-t[1] * (-n[1]) - d * n[0], -t[1] * n[0] - d * n[1]], board)
1623                     );
1624             }
1625 
1626             return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
1627         },
1628 
1629         /**
1630          * Intersection of two circles.
1631          * @param {Array} circ1 stdform of the first circle
1632          * @param {Array} circ2 stdform of the second circle
1633          * @param {number} i number of the returned intersection point.
1634          *   i==0: use the positive square root,
1635          *   i==1: use the negative square root.
1636          * @param {JXG.Board} board Reference to the board.
1637          * @returns {JXG.Coords} Coordinates of the intersection point
1638          */
1639         meetCircleCircle: function (circ1, circ2, i, board) {
1640             var radicalAxis;
1641 
1642             // Radius is zero, return center of circle, if on other circle
1643             if (circ1[4] < Mat.eps) {
1644                 if (Math.abs(this.distance(circ1.slice(6, 2), circ2.slice(6, 8)) - circ2[4]) < Mat.eps) {
1645                     return new Coords(Const.COORDS_BY_USER, circ1.slice(6, 8), board);
1646                 }
1647 
1648                 return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
1649             }
1650 
1651             // Radius is zero, return center of circle, if on other circle
1652             if (circ2[4] < Mat.eps) {
1653                 if (Math.abs(this.distance(circ2.slice(6, 2), circ1.slice(6, 8)) - circ1[4]) < Mat.eps) {
1654                     return new Coords(Const.COORDS_BY_USER, circ2.slice(6, 8), board);
1655                 }
1656 
1657                 return new Coords(Const.COORDS_BY_USER, [0, 0, 0], board);
1658             }
1659 
1660             radicalAxis = [circ2[3] * circ1[0] - circ1[3] * circ2[0],
1661                 circ2[3] * circ1[1] - circ1[3] * circ2[1],
1662                 circ2[3] * circ1[2] - circ1[3] * circ2[2],
1663                 0, 1, Infinity, Infinity, Infinity];
1664             radicalAxis = Mat.normalize(radicalAxis);
1665 
1666             return this.meetLineCircle(radicalAxis, circ1, i, board);
1667         },
1668 
1669         /**
1670          * Compute an intersection of the curves c1 and c2.
1671          * We want to find values t1, t2 such that
1672          * c1(t1) = c2(t2), i.e. (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
1673          *
1674          * Methods: segment-wise intersections (default) or generalized Newton method.
1675          * @param {JXG.Curve} c1 Curve, Line or Circle
1676          * @param {JXG.Curve} c2 Curve, Line or Circle
1677          * @param {Number} nr the nr-th intersection point will be returned.
1678          * @param {Number} t2ini not longer used.
1679          * @param {JXG.Board} [board=c1.board] Reference to a board object.
1680          * @param {String} [method='segment'] Intersection method, possible values are 'newton' and 'segment'.
1681          * @returns {JXG.Coords} intersection point
1682          */
1683         meetCurveCurve: function (c1, c2, nr, t2ini, board, method) {
1684             var co;
1685 
1686             if (Type.exists(method) && method === 'newton') {
1687                 co = Numerics.generalizedNewton(c1, c2, nr, t2ini);
1688             } else {
1689                 if (c1.bezierDegree === 3 || c2.bezierDegree === 3) {
1690                     co = this.meetBezierCurveRedBlueSegments(c1, c2, nr);
1691                 } else {
1692                     co = this.meetCurveRedBlueSegments(c1, c2, nr);
1693                 }
1694             }
1695 
1696             return (new Coords(Const.COORDS_BY_USER, co, board));
1697         },
1698 
1699         /**
1700          * Intersection of curve with line,
1701          * Order of input does not matter for el1 and el2.
1702          * From version 0.99.7 on this method calls
1703          * {@link JXG.Math.Geometry.meetCurveLineDiscrete}.
1704          * If higher precision is needed, {@link JXG.Math.Geometry.meetCurveLineContinuous}
1705          * has to be used.
1706          *
1707          * @param {JXG.Curve,JXG.Line} el1 Curve or Line
1708          * @param {JXG.Curve,JXG.Line} el2 Curve or Line
1709          * @param {Number} nr the nr-th intersection point will be returned.
1710          * @param {JXG.Board} [board=el1.board] Reference to a board object.
1711          * @param {Boolean} alwaysIntersect If false just the segment between the two defining points are tested for intersection
1712          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
1713          * the ideal point [0,1,0] is returned.
1714          */
1715         meetCurveLine: function (el1, el2, nr, board, alwaysIntersect) {
1716             var v = [0, NaN, NaN], cu, li;
1717 
1718             if (!Type.exists(board)) {
1719                 board = el1.board;
1720             }
1721 
1722             if (el1.elementClass === Const.OBJECT_CLASS_CURVE) {
1723                 cu = el1;
1724                 li = el2;
1725             } else {
1726                 cu = el2;
1727                 li = el1;
1728             }
1729 
1730             v = this.meetCurveLineDiscrete(cu, li, nr, board, !alwaysIntersect);
1731 
1732             return v;
1733         },
1734 
1735         /**
1736          * Intersection of line and curve, continuous case.
1737          * Finds the nr-the intersection point
1738          * Uses {@link JXG.Math.Geometry.meetCurveLineDiscrete} as a first approximation.
1739          * A more exact solution is then found with {@link JXG.Math.Numerics.root}.
1740          *
1741          * @param {JXG.Curve} cu Curve
1742          * @param {JXG.Line} li Line
1743          * @param {Number} nr Will return the nr-th intersection point.
1744          * @param {JXG.Board} board
1745          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the
1746          * line defined by the segment
1747          * @returns {JXG.Coords} Coords object containing the intersection.
1748          */
1749         meetCurveLineContinuous: function (cu, li, nr, board, testSegment) {
1750             var t, func0, func1, v, x, y, z,
1751                 eps = Mat.eps,
1752                 epsLow = Mat.eps,
1753                 steps, delta, tnew, i,
1754                 tmin, fmin, ft;
1755 
1756             v = this.meetCurveLineDiscrete(cu, li, nr, board, testSegment);
1757             x = v.usrCoords[1];
1758             y = v.usrCoords[2];
1759 
1760             func0 = function (t) {
1761                 var c1, c2;
1762 
1763                 if (t > cu.maxX() || t < cu.minX()) {
1764                     return Infinity;
1765                 }
1766                 c1 = x - cu.X(t);
1767                 c2 = y - cu.Y(t);
1768                 return c1 * c1 + c2 * c2;
1769             };
1770 
1771             func1 = function (t) {
1772                 var v = li.stdform[0] + li.stdform[1] * cu.X(t) + li.stdform[2] * cu.Y(t);
1773                 return v * v;
1774             };
1775 
1776             // Find t
1777             steps = 50;
1778             delta = (cu.maxX() - cu.minX()) / steps;
1779             tnew = cu.minX();
1780 
1781             fmin = 0.0001; //eps;
1782             tmin = NaN;
1783             for (i = 0; i < steps; i++) {
1784                 t = Numerics.root(func0, [Math.max(tnew, cu.minX()), Math.min(tnew + delta, cu.maxX())]);
1785                 ft = Math.abs(func0(t));
1786                 if (ft <= fmin) {
1787                     fmin = ft;
1788                     tmin = t;
1789                     if (fmin < eps) {
1790                         break;
1791                     }
1792                 }
1793 
1794                 tnew += delta;
1795             }
1796             t = tmin;
1797             // Compute "exact" t
1798             t = Numerics.root(func1, [Math.max(t - delta, cu.minX()), Math.min(t + delta, cu.maxX())]);
1799 
1800             ft = func1(t);
1801             // Is the point on the line?
1802             if (isNaN(ft) || Math.abs(ft) > epsLow) {
1803                 z = 0.0; //NaN;
1804             } else {
1805                 z = 1.0;
1806             }
1807 
1808             return (new Coords(Const.COORDS_BY_USER, [z, cu.X(t), cu.Y(t)], board));
1809         },
1810 
1811         /**
1812          * Intersection of line and curve, discrete case.
1813          * Segments are treated as lines.
1814          * Finding the nr-th intersection point should work for all nr.
1815          * @param {JXG.Curve} cu
1816          * @param {JXG.Line} li
1817          * @param {Number} nr
1818          * @param {JXG.Board} board
1819          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the
1820          * line defined by the segment
1821          *
1822          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
1823          * the ideal point [0,1,0] is returned.
1824          */
1825         meetCurveLineDiscrete: function (cu, li, nr, board, testSegment) {
1826             var i, j,
1827                 p1, p2, p, q,
1828                 lip1 = li.point1.coords.usrCoords,
1829                 lip2 = li.point2.coords.usrCoords,
1830                 d, res,
1831                 cnt = 0,
1832                 len = cu.numberPoints,
1833                 ev_sf = Type.evaluate(li.visProp.straightfirst),
1834                 ev_sl = Type.evaluate(li.visProp.straightlast);
1835 
1836             // In case, no intersection will be found we will take this
1837             q = new Coords(Const.COORDS_BY_USER, [0, NaN, NaN], board);
1838 
1839             if (lip1[0] === 0.0) {
1840                 lip1 = [1, lip2[1] + li.stdform[2], lip2[2] - li.stdform[1]];
1841             } else if (lip2[0] === 0.0) {
1842                 lip2 = [1, lip1[1] + li.stdform[2], lip1[2] - li.stdform[1]];
1843             }
1844 
1845             p2 = cu.points[0].usrCoords;
1846             for (i = 1; i < len; i += cu.bezierDegree) {
1847                 p1 = p2.slice(0);
1848                 p2 = cu.points[i].usrCoords;
1849                 d = this.distance(p1, p2);
1850 
1851                 // The defining points are not identical
1852                 if (d > Mat.eps) {
1853                     if (cu.bezierDegree === 3) {
1854                         res = this.meetBeziersegmentBeziersegment([
1855                             cu.points[i - 1].usrCoords.slice(1),
1856                             cu.points[i].usrCoords.slice(1),
1857                             cu.points[i + 1].usrCoords.slice(1),
1858                             cu.points[i + 2].usrCoords.slice(1)
1859                         ], [
1860                             lip1.slice(1),
1861                             lip2.slice(1)
1862                         ], testSegment);
1863                     } else {
1864                         res = [this.meetSegmentSegment(p1, p2, lip1, lip2)];
1865                     }
1866 
1867                     for (j = 0; j < res.length; j++) {
1868                         p = res[j];
1869                         if (0 <= p[1] && p[1] <= 1) {
1870                             if (cnt === nr) {
1871                                 /**
1872                                 * If the intersection point is not part of the segment,
1873                                 * this intersection point is set to non-existent.
1874                                 * This prevents jumping behavior of the intersection points.
1875                                 * But it may be discussed if it is the desired behavior.
1876                                 */
1877                                 if (testSegment &&
1878                                         ((!ev_sf && p[2] < 0) || (!ev_sl && p[2] > 1))) {
1879                                     return q;  // break;
1880                                 }
1881 
1882                                 q = new Coords(Const.COORDS_BY_USER, p[0], board);
1883                                 return q;      // break;
1884                             }
1885                             cnt += 1;
1886                         }
1887                     }
1888                 }
1889             }
1890 
1891             return q;
1892         },
1893 
1894         /**
1895          * Find the n-th intersection point of two curves named red (first parameter) and blue (second parameter).
1896          * We go through each segment of the red curve and search if there is an intersection with a segemnt of the blue curve.
1897          * This double loop, i.e. the outer loop runs along the red curve and the inner loop runs along the blue curve, defines
1898          * the n-th intersection point. The segments are either line segments or Bezier curves of degree 3. This depends on
1899          * the property bezierDegree of the curves.
1900          * <p>
1901          * This method works also for transformed curves, since only the already
1902          * transformed points are used.
1903          *
1904          * @param {JXG.Curve} red
1905          * @param {JXG.Curve} blue
1906          * @param {Number} nr
1907          */
1908         meetCurveRedBlueSegments: function (red, blue, nr) {
1909             var i, j,
1910                 red1, red2, blue1, blue2, m,
1911                 minX, maxX,
1912                 iFound = 0,
1913                 lenBlue = blue.numberPoints, //points.length,
1914                 lenRed = red.numberPoints; //points.length;
1915 
1916             if (lenBlue <= 1 || lenRed <= 1) {
1917                 return [0, NaN, NaN];
1918             }
1919 
1920             for (i = 1; i < lenRed; i++) {
1921                 red1 = red.points[i - 1].usrCoords;
1922                 red2 = red.points[i].usrCoords;
1923                 minX = Math.min(red1[1], red2[1]);
1924                 maxX = Math.max(red1[1], red2[1]);
1925 
1926                 blue2 = blue.points[0].usrCoords;
1927                 for (j = 1; j < lenBlue; j++) {
1928                     blue1 = blue2;
1929                     blue2 = blue.points[j].usrCoords;
1930 
1931                     if (Math.min(blue1[1], blue2[1]) < maxX && Math.max(blue1[1], blue2[1]) > minX) {
1932                         m = this.meetSegmentSegment(red1, red2, blue1, blue2);
1933                         if (m[1] >= 0.0 && m[2] >= 0.0 &&
1934                                 // The two segments meet in the interior or at the start points
1935                                 ((m[1] < 1.0 && m[2] < 1.0) ||
1936                                 // One of the curve is intersected in the very last point
1937                                 (i === lenRed - 1 && m[1] === 1.0) ||
1938                                 (j === lenBlue - 1 && m[2] === 1.0))) {
1939                             if (iFound === nr) {
1940                                 return m[0];
1941                             }
1942 
1943                             iFound++;
1944                         }
1945                     }
1946                 }
1947             }
1948 
1949             return [0, NaN, NaN];
1950         },
1951 
1952         /**
1953          * (Virtual) Intersection of two segments.
1954          * @param {Array} p1 First point of segment 1 using normalized homogeneous coordinates [1,x,y]
1955          * @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
1956          * @param {Array} q1 First point of segment 2 using normalized homogeneous coordinates [1,x,y]
1957          * @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
1958          * @returns {Array} [Intersection point, t, u] The first entry contains the homogeneous coordinates
1959          * of the intersection point. The second and third entry give the position of the intersection with respect
1960          * to the definiting parameters. For example, the second entry t is defined by: intersection point = p1 + t * deltaP, where
1961          * deltaP = (p2 - p1) when both parameters are coordinates, and deltaP = p2 if p2 is a point at infinity.
1962          * If the two segments are collinear, [[0,0,0], Infinity, Infinity] is returned.
1963          **/
1964         meetSegmentSegment: function (p1, p2, q1, q2) {
1965             var t, u, i, d,
1966                 li1 = Mat.crossProduct(p1, p2),
1967                 li2 = Mat.crossProduct(q1, q2),
1968                 c = Mat.crossProduct(li1, li2);
1969 
1970             if (Math.abs(c[0]) < Mat.eps) {
1971                 return [c, Infinity, Infinity];
1972             }
1973 
1974             // Normalize the intersection coordinates
1975             c[1] /= c[0];
1976             c[2] /= c[0];
1977             c[0] /= c[0];
1978 
1979             // Now compute in principle:
1980             //    t = dist(c - p1) / dist(p2 - p1) and
1981             //    u = dist(c - q1) / dist(q2 - q1)
1982             // However: the points q1, q2, p1, p2 might be ideal points - or in general - the
1983             // coordinates might be not normalized.
1984             // Note that the z-coordinates of p2 and q2 are used to determine whether it should be interpreted
1985             // as a segment coordinate or a direction.
1986             i = (Math.abs(p2[1] - p2[0] * p1[1]) < Mat.eps) ? 2 : 1;
1987             d = p1[i] / p1[0];
1988             t = (c[i] - d) / ( (p2[0] !== 0) ? (p2[i] / p2[0] - d) : p2[i] );
1989 
1990             i = (Math.abs(q2[1] - q2[0] * q1[1]) < Mat.eps) ? 2 : 1;
1991             d = q1[i] / q1[0];
1992             u = (c[i] - d) / ( (q2[0] !== 0) ? (q2[i] / q2[0] - d) : q2[i] );
1993 
1994             return [c, t, u];
1995         },
1996 
1997         /**
1998          * Find the n-th intersection point of two pathes, usually given by polygons. Uses parts of the
1999          * Greiner-Hormann algorithm in JXG.Math.Clip.
2000          *
2001          * @param {JXG.Circle|JXG.Curve|JXG.Polygon} path1
2002          * @param {JXG.Circle|JXG.Curve|JXG.Polygon} path2
2003          * @param {Number} n
2004          * @param {JXG.Board} board
2005          *
2006          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
2007          * the ideal point [0,0,0] is returned.
2008          *
2009          */
2010         meetPathPath: function(path1, path2, nr, board) {
2011             var S, C, len, intersections;
2012 
2013             S = JXG.Math.Clip._getPath(path1, board);
2014             len = S.length;
2015             if (len > 0 && this.distance(S[0].coords.usrCoords, S[len - 1].coords.usrCoords, 3) < Mat.eps) {
2016                 S.pop();
2017             }
2018 
2019             C = JXG.Math.Clip._getPath(path2, board);
2020             len = C.length;
2021             if (len > 0 && this.distance(C[0].coords.usrCoords, C[len - 1].coords.usrCoords, 3) < Mat.eps * Mat.eps) {
2022                 C.pop();
2023             }
2024 
2025             // Handle cases where at least one of the paths is empty
2026             if (nr < 0 || JXG.Math.Clip.isEmptyCase(S, C, 'intersection')) {
2027                 return (new Coords(Const.COORDS_BY_USER, [0, 0, 0], board));
2028             }
2029 
2030             JXG.Math.Clip.makeDoublyLinkedList(S);
2031             JXG.Math.Clip.makeDoublyLinkedList(C);
2032 
2033             intersections = JXG.Math.Clip.findIntersections(S, C, board)[0];
2034             if (nr < intersections.length) {
2035                 return intersections[nr].coords;
2036             }
2037             return (new Coords(Const.COORDS_BY_USER, [0, 0, 0], board));
2038         },
2039 
2040         /**
2041          * Find the n-th intersection point between a polygon and a line.
2042          * @param {JXG.Polygon} path
2043          * @param {JXG.Line} line
2044          * @param {Number} nr
2045          * @param {JXG.Board} board
2046          * @param {Boolean} alwaysIntersect If false just the segment between the two defining points of the line are tested for intersection.
2047          *
2048          * @returns {JXG.Coords} Intersection point. In case no intersection point is detected,
2049          * the ideal point [0,0,0] is returned.
2050          */
2051         meetPolygonLine: function(path, line, nr, board, alwaysIntersect) {
2052             var i, res, border,
2053                 crds = [0,0,0],
2054                 len = path.borders.length,
2055                 intersections = [];
2056 
2057             for (i = 0; i < len; i++) {
2058                 border = path.borders[i];
2059                 res = this.meetSegmentSegment(
2060                     border.point1.coords.usrCoords,
2061                     border.point2.coords.usrCoords,
2062                     line.point1.coords.usrCoords,
2063                     line.point2.coords.usrCoords);
2064 
2065                 if (
2066                     (!alwaysIntersect || (res[2] >= 0 && res[2] < 1)) &&
2067                     res[1] >= 0 && res[1] < 1) {
2068                     intersections.push(res[0]);
2069                 }
2070             }
2071 
2072             if (nr >= 0 && nr < intersections.length) {
2073                 crds = intersections[nr];
2074             }
2075             return (new Coords(Const.COORDS_BY_USER, crds, board));
2076         },
2077 
2078         /****************************************/
2079         /****   BEZIER CURVE ALGORITHMS      ****/
2080         /****************************************/
2081 
2082         /**
2083          * Splits a Bezier curve segment defined by four points into
2084          * two Bezier curve segments. Dissection point is t=1/2.
2085          * @param {Array} curve Array of four coordinate arrays of length 2 defining a
2086          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2087          * @returns {Array} Array consisting of two coordinate arrays for Bezier curves.
2088          */
2089         _bezierSplit: function (curve) {
2090             var p0, p1, p2, p00, p22, p000;
2091 
2092             p0 = [(curve[0][0] + curve[1][0]) * 0.5, (curve[0][1] + curve[1][1]) * 0.5];
2093             p1 = [(curve[1][0] + curve[2][0]) * 0.5, (curve[1][1] + curve[2][1]) * 0.5];
2094             p2 = [(curve[2][0] + curve[3][0]) * 0.5, (curve[2][1] + curve[3][1]) * 0.5];
2095 
2096             p00 = [(p0[0] + p1[0]) * 0.5, (p0[1] + p1[1]) * 0.5];
2097             p22 = [(p1[0] + p2[0]) * 0.5, (p1[1] + p2[1]) * 0.5];
2098 
2099             p000 = [(p00[0] + p22[0]) * 0.5, (p00[1] + p22[1]) * 0.5];
2100 
2101             return [[curve[0], p0, p00, p000], [p000, p22, p2, curve[3]]];
2102         },
2103 
2104         /**
2105          * Computes the bounding box [minX, maxY, maxX, minY] of a Bezier curve segment
2106          * from its control points.
2107          * @param {Array} curve Array of four coordinate arrays of length 2 defining a
2108          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2109          * @returns {Array} Bounding box [minX, maxY, maxX, minY]
2110          */
2111         _bezierBbox: function (curve) {
2112             var bb = [];
2113 
2114             if (curve.length === 4) {   // bezierDegree == 3
2115                 bb[0] = Math.min(curve[0][0], curve[1][0], curve[2][0], curve[3][0]); // minX
2116                 bb[1] = Math.max(curve[0][1], curve[1][1], curve[2][1], curve[3][1]); // maxY
2117                 bb[2] = Math.max(curve[0][0], curve[1][0], curve[2][0], curve[3][0]); // maxX
2118                 bb[3] = Math.min(curve[0][1], curve[1][1], curve[2][1], curve[3][1]); // minY
2119             } else {                   // bezierDegree == 1
2120                 bb[0] = Math.min(curve[0][0], curve[1][0]); // minX
2121                 bb[1] = Math.max(curve[0][1], curve[1][1]); // maxY
2122                 bb[2] = Math.max(curve[0][0], curve[1][0]); // maxX
2123                 bb[3] = Math.min(curve[0][1], curve[1][1]); // minY
2124             }
2125 
2126             return bb;
2127         },
2128 
2129         /**
2130          * Decide if two Bezier curve segments overlap by comparing their bounding boxes.
2131          * @param {Array} bb1 Bounding box of the first Bezier curve segment
2132          * @param {Array} bb2 Bounding box of the second Bezier curve segment
2133          * @returns {Boolean} true if the bounding boxes overlap, false otherwise.
2134          */
2135         _bezierOverlap: function (bb1, bb2) {
2136             return bb1[2] >= bb2[0] && bb1[0] <= bb2[2] && bb1[1] >= bb2[3] && bb1[3] <= bb2[1];
2137         },
2138 
2139         /**
2140          * Append list of intersection points to a list.
2141          * @private
2142          */
2143         _bezierListConcat: function (L, Lnew, t1, t2) {
2144             var i,
2145                 t2exists = Type.exists(t2),
2146                 start = 0,
2147                 len = Lnew.length,
2148                 le = L.length;
2149 
2150             if (le > 0 && len > 0 &&
2151                     ((L[le - 1][1] === 1 && Lnew[0][1] === 0) ||
2152                     (t2exists && L[le - 1][2] === 1 && Lnew[0][2] === 0))) {
2153                 start = 1;
2154             }
2155 
2156             for (i = start; i < len; i++) {
2157                 if (t2exists) {
2158                     Lnew[i][2] *= 0.5;
2159                     Lnew[i][2] += t2;
2160                 }
2161 
2162                 Lnew[i][1] *= 0.5;
2163                 Lnew[i][1] += t1;
2164 
2165                 L.push(Lnew[i]);
2166             }
2167         },
2168 
2169         /**
2170          * Find intersections of two Bezier curve segments by recursive subdivision.
2171          * Below maxlevel determine intersections by intersection line segments.
2172          * @param {Array} red Array of four coordinate arrays of length 2 defining the first
2173          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2174          * @param {Array} blue Array of four coordinate arrays of length 2 defining the second
2175          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2176          * @param {Number} level Recursion level
2177          * @returns {Array} List of intersection points (up to nine). Each intersection point is an
2178          * array of length three (homogeneous coordinates) plus preimages.
2179          */
2180         _bezierMeetSubdivision: function (red, blue, level) {
2181             var bbb, bbr,
2182                 ar, b0, b1, r0, r1, m,
2183                 p0, p1, q0, q1,
2184                 L = [],
2185                 maxLev = 5;      // Maximum recursion level
2186 
2187             bbr = this._bezierBbox(blue);
2188             bbb = this._bezierBbox(red);
2189 
2190             if (!this._bezierOverlap(bbr, bbb)) {
2191                 return [];
2192             }
2193 
2194             if (level < maxLev) {
2195                 ar = this._bezierSplit(red);
2196                 r0 = ar[0];
2197                 r1 = ar[1];
2198 
2199                 ar = this._bezierSplit(blue);
2200                 b0 = ar[0];
2201                 b1 = ar[1];
2202 
2203                 this._bezierListConcat(L, this._bezierMeetSubdivision(r0, b0, level + 1), 0.0, 0.0);
2204                 this._bezierListConcat(L, this._bezierMeetSubdivision(r0, b1, level + 1), 0, 0.5);
2205                 this._bezierListConcat(L, this._bezierMeetSubdivision(r1, b0, level + 1), 0.5, 0.0);
2206                 this._bezierListConcat(L, this._bezierMeetSubdivision(r1, b1, level + 1), 0.5, 0.5);
2207 
2208                 return L;
2209             }
2210 
2211             // Make homogeneous coordinates
2212             q0 = [1].concat(red[0]);
2213             q1 = [1].concat(red[3]);
2214             p0 = [1].concat(blue[0]);
2215             p1 = [1].concat(blue[3]);
2216 
2217             m = this.meetSegmentSegment(q0, q1, p0, p1);
2218 
2219             if (m[1] >= 0.0 && m[2] >= 0.0 && m[1] <= 1.0 && m[2] <= 1.0) {
2220                 return [m];
2221             }
2222 
2223             return [];
2224         },
2225 
2226         /**
2227          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the line defined by the segment
2228          */
2229         _bezierLineMeetSubdivision: function (red, blue, level, testSegment) {
2230             var bbb, bbr,
2231                 ar, r0, r1, m,
2232                 p0, p1, q0, q1,
2233                 L = [],
2234                 maxLev = 5;      // Maximum recursion level
2235 
2236             bbb = this._bezierBbox(blue);
2237             bbr = this._bezierBbox(red);
2238 
2239             if (testSegment && !this._bezierOverlap(bbr, bbb)) {
2240                 return [];
2241             }
2242 
2243             if (level < maxLev) {
2244                 ar = this._bezierSplit(red);
2245                 r0 = ar[0];
2246                 r1 = ar[1];
2247 
2248                 this._bezierListConcat(L, this._bezierLineMeetSubdivision(r0, blue, level + 1), 0.0);
2249                 this._bezierListConcat(L, this._bezierLineMeetSubdivision(r1, blue, level + 1), 0.5);
2250 
2251                 return L;
2252             }
2253 
2254             // Make homogeneous coordinates
2255             q0 = [1].concat(red[0]);
2256             q1 = [1].concat(red[3]);
2257             p0 = [1].concat(blue[0]);
2258             p1 = [1].concat(blue[1]);
2259 
2260             m = this.meetSegmentSegment(q0, q1, p0, p1);
2261 
2262             if (m[1] >= 0.0 && m[1] <= 1.0) {
2263                 if (!testSegment || (m[2] >= 0.0 && m[2] <= 1.0)) {
2264                     return [m];
2265                 }
2266             }
2267 
2268             return [];
2269         },
2270 
2271         /**
2272          * Find the nr-th intersection point of two Bezier curve segments.
2273          * @param {Array} red Array of four coordinate arrays of length 2 defining the first
2274          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2275          * @param {Array} blue Array of four coordinate arrays of length 2 defining the second
2276          * Bezier curve segment, i.e. [[x0,y0], [x1,y1], [x2,y2], [x3,y3]].
2277          * @param {Boolean} testSegment Test if intersection has to be inside of the segment or somewhere on the line defined by the segment
2278          * @returns {Array} Array containing the list of all intersection points as homogeneous coordinate arrays plus
2279          * preimages [x,y], t_1, t_2] of the two Bezier curve segments.
2280          *
2281          */
2282         meetBeziersegmentBeziersegment: function (red, blue, testSegment) {
2283             var L, L2, i;
2284 
2285             if (red.length === 4 && blue.length === 4) {
2286                 L = this._bezierMeetSubdivision(red, blue, 0);
2287             } else {
2288                 L = this._bezierLineMeetSubdivision(red, blue, 0, testSegment);
2289             }
2290 
2291             L.sort(function (a, b) {
2292                 return (a[1] - b[1]) * 10000000.0 + (a[2] - b[2]);
2293             });
2294 
2295             L2 = [];
2296             for (i = 0; i < L.length; i++) {
2297                 // Only push entries different from their predecessor
2298                 if (i === 0 || (L[i][1] !== L[i - 1][1] || L[i][2] !== L[i - 1][2])) {
2299                     L2.push(L[i]);
2300                 }
2301             }
2302             return L2;
2303         },
2304 
2305         /**
2306          * Find the nr-th intersection point of two Bezier curves, i.e. curves with bezierDegree == 3.
2307          * @param {JXG.Curve} red Curve with bezierDegree == 3
2308          * @param {JXG.Curve} blue Curve with bezierDegree == 3
2309          * @param {Number} nr The number of the intersection point which should be returned.
2310          * @returns {Array} The homogeneous coordinates of the nr-th intersection point.
2311          */
2312         meetBezierCurveRedBlueSegments: function (red, blue, nr) {
2313             var p, i, j, k, po,
2314                 redArr, blueArr,
2315                 bbr, bbb, intersections,
2316                 startRed = 0,
2317                 startBlue = 0,
2318                 lenBlue = blue.numberPoints,
2319                 lenRed = red.numberPoints,
2320                 L = [];
2321 
2322             if (lenBlue < blue.bezierDegree + 1 || lenRed < red.bezierDegree + 1) {
2323                 return [0, NaN, NaN];
2324             }
2325             lenBlue -= blue.bezierDegree;
2326             lenRed  -= red.bezierDegree;
2327 
2328             // For sectors, we ignore the "legs"
2329             if (red.type === Const.OBJECT_TYPE_SECTOR) {
2330                 startRed = 3;
2331                 lenRed  -= 3;
2332             }
2333             if (blue.type === Const.OBJECT_TYPE_SECTOR) {
2334                 startBlue = 3;
2335                 lenBlue  -= 3;
2336             }
2337 
2338             for (i = startRed; i < lenRed; i += red.bezierDegree) {
2339                 p = red.points;
2340                 redArr = [
2341                     p[i].usrCoords.slice(1),
2342                     p[i + 1].usrCoords.slice(1)
2343                 ];
2344                 if (red.bezierDegree === 3) {
2345                     redArr[2] = p[i + 2].usrCoords.slice(1);
2346                     redArr[3] = p[i + 3].usrCoords.slice(1);
2347                 }
2348 
2349                 bbr = this._bezierBbox(redArr);
2350 
2351                 for (j = startBlue; j < lenBlue; j += blue.bezierDegree) {
2352                     p = blue.points;
2353                     blueArr = [
2354                         p[j].usrCoords.slice(1),
2355                         p[j + 1].usrCoords.slice(1)
2356                     ];
2357                     if (blue.bezierDegree === 3) {
2358                         blueArr[2] = p[j + 2].usrCoords.slice(1);
2359                         blueArr[3] = p[j + 3].usrCoords.slice(1);
2360                     }
2361 
2362                     bbb = this._bezierBbox(blueArr);
2363                     if (this._bezierOverlap(bbr, bbb)) {
2364                         intersections = this.meetBeziersegmentBeziersegment(redArr, blueArr);
2365                         if (intersections.length === 0) {
2366                             continue;
2367                         }
2368                         for (k = 0; k < intersections.length; k++) {
2369                             po = intersections[k];
2370                             if (po[1] < -Mat.eps ||
2371                                 po[1] > 1 + Mat.eps ||
2372                                 po[2] < -Mat.eps ||
2373                                 po[2] > 1 + Mat.eps) {
2374                                 continue;
2375                             }
2376                             L.push(po);
2377                         }
2378                         if (L.length > nr) {
2379                             return L[nr][0];
2380                         }
2381                     }
2382                 }
2383             }
2384             if (L.length > nr) {
2385                 return L[nr][0];
2386             }
2387 
2388             return [0, NaN, NaN];
2389         },
2390 
2391         bezierSegmentEval: function (t, curve) {
2392             var f, x, y,
2393                 t1 = 1.0 - t;
2394 
2395             x = 0;
2396             y = 0;
2397 
2398             f = t1 * t1 * t1;
2399             x += f * curve[0][0];
2400             y += f * curve[0][1];
2401 
2402             f = 3.0 * t * t1 * t1;
2403             x += f * curve[1][0];
2404             y += f * curve[1][1];
2405 
2406             f = 3.0 * t * t * t1;
2407             x += f * curve[2][0];
2408             y += f * curve[2][1];
2409 
2410             f = t * t * t;
2411             x += f * curve[3][0];
2412             y += f * curve[3][1];
2413 
2414             return [1.0, x, y];
2415         },
2416 
2417         /**
2418          * Generate the defining points of a 3rd degree bezier curve that approximates
2419          * a circle sector defined by three coordinate points A, B, C, each defined by an array of length three.
2420          * The coordinate arrays are given in homogeneous coordinates.
2421          * @param {Array} A First point
2422          * @param {Array} B Second point (intersection point)
2423          * @param {Array} C Third point
2424          * @param {Boolean} withLegs Flag. If true the legs to the intersection point are part of the curve.
2425          * @param {Number} sgn Wither 1 or -1. Needed for minor and major arcs. In case of doubt, use 1.
2426          */
2427         bezierArc: function (A, B, C, withLegs, sgn) {
2428             var p1, p2, p3, p4,
2429                 r, phi, beta,
2430                 PI2 = Math.PI * 0.5,
2431                 x = B[1],
2432                 y = B[2],
2433                 z = B[0],
2434                 dataX = [], dataY = [],
2435                 co, si, ax, ay, bx, by, k, v, d, matrix;
2436 
2437             r = this.distance(B, A);
2438 
2439             // x,y, z is intersection point. Normalize it.
2440             x /= z;
2441             y /= z;
2442 
2443             phi = this.rad(A.slice(1), B.slice(1), C.slice(1));
2444             if (sgn === -1) {
2445                 phi = 2 * Math.PI - phi;
2446             }
2447 
2448             p1 = A;
2449             p1[1] /= p1[0];
2450             p1[2] /= p1[0];
2451             p1[0] /= p1[0];
2452 
2453             p4 = p1.slice(0);
2454 
2455             if (withLegs) {
2456                 dataX = [x, x + 0.333 * (p1[1] - x), x + 0.666 * (p1[1] - x), p1[1]];
2457                 dataY = [y, y + 0.333 * (p1[2] - y), y + 0.666 * (p1[2] - y), p1[2]];
2458             } else {
2459                 dataX = [p1[1]];
2460                 dataY = [p1[2]];
2461             }
2462 
2463             while (phi > Mat.eps) {
2464                 if (phi > PI2) {
2465                     beta = PI2;
2466                     phi -= PI2;
2467                 } else {
2468                     beta = phi;
2469                     phi = 0;
2470                 }
2471 
2472                 co = Math.cos(sgn * beta);
2473                 si = Math.sin(sgn * beta);
2474 
2475                 matrix = [
2476                     [1, 0, 0],
2477                     [x * (1 - co) + y * si, co, -si],
2478                     [y * (1 - co) - x * si, si,  co]
2479                 ];
2480                 v = Mat.matVecMult(matrix, p1);
2481                 p4 = [v[0] / v[0], v[1] / v[0], v[2] / v[0]];
2482 
2483                 ax = p1[1] - x;
2484                 ay = p1[2] - y;
2485                 bx = p4[1] - x;
2486                 by = p4[2] - y;
2487 
2488                 d = Math.sqrt((ax + bx) * (ax + bx) + (ay + by) * (ay + by));
2489 
2490                 if (Math.abs(by - ay) > Mat.eps) {
2491                     k = (ax + bx) * (r / d - 0.5) / (by - ay) * 8 / 3;
2492                 } else {
2493                     k = (ay + by) * (r / d - 0.5) / (ax - bx) * 8 / 3;
2494                 }
2495 
2496                 p2 = [1, p1[1] - k * ay, p1[2] + k * ax];
2497                 p3 = [1, p4[1] + k * by, p4[2] - k * bx];
2498 
2499                 dataX = dataX.concat([p2[1], p3[1], p4[1]]);
2500                 dataY = dataY.concat([p2[2], p3[2], p4[2]]);
2501                 p1 = p4.slice(0);
2502             }
2503 
2504             if (withLegs) {
2505                 dataX = dataX.concat([ p4[1] + 0.333 * (x - p4[1]), p4[1] + 0.666 * (x - p4[1]), x]);
2506                 dataY = dataY.concat([ p4[2] + 0.333 * (y - p4[2]), p4[2] + 0.666 * (y - p4[2]), y]);
2507             }
2508 
2509             return [dataX, dataY];
2510         },
2511 
2512         /****************************************/
2513         /****           PROJECTIONS          ****/
2514         /****************************************/
2515 
2516         /**
2517          * Calculates the coordinates of the projection of a given point on a given circle. I.o.w. the
2518          * nearest one of the two intersection points of the line through the given point and the circles
2519          * center.
2520          * @param {JXG.Point,JXG.Coords} point Point to project or coords object to project.
2521          * @param {JXG.Circle} circle Circle on that the point is projected.
2522          * @param {JXG.Board} [board=point.board] Reference to the board
2523          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given circle.
2524          */
2525         projectPointToCircle: function (point, circle, board) {
2526             var dist, P, x, y, factor,
2527                 M = circle.center.coords.usrCoords;
2528 
2529             if (!Type.exists(board)) {
2530                 board = point.board;
2531             }
2532 
2533             // gave us a point
2534             if (Type.isPoint(point)) {
2535                 dist = point.coords.distance(Const.COORDS_BY_USER, circle.center.coords);
2536                 P = point.coords.usrCoords;
2537             // gave us coords
2538             } else {
2539                 dist = point.distance(Const.COORDS_BY_USER, circle.center.coords);
2540                 P = point.usrCoords;
2541             }
2542 
2543             if (Math.abs(dist) < Mat.eps) {
2544                 dist = Mat.eps;
2545             }
2546 
2547             factor = circle.Radius() / dist;
2548             x = M[1] + factor * (P[1] - M[1]);
2549             y = M[2] + factor * (P[2] - M[2]);
2550 
2551             return new Coords(Const.COORDS_BY_USER, [x, y], board);
2552         },
2553 
2554         /**
2555          * Calculates the coordinates of the orthogonal projection of a given point on a given line. I.o.w. the
2556          * intersection point of the given line and its perpendicular through the given point.
2557          * @param {JXG.Point|JXG.Coords} point Point to project.
2558          * @param {JXG.Line} line Line on that the point is projected.
2559          * @param {JXG.Board} [board=point.board|board=line.board] Reference to a board.
2560          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given line.
2561          */
2562         projectPointToLine: function (point, line, board) {
2563             var v = [0, line.stdform[1], line.stdform[2]],
2564                 coords;
2565 
2566             if (!Type.exists(board)) {
2567                 if (Type.exists(point.coords)) {
2568                     board = point.board;
2569                 } else {
2570                     board = line.board;
2571                 }
2572             }
2573 
2574             if (Type.exists(point.coords)) {
2575                 coords = point.coords.usrCoords;
2576             } else {
2577                 coords = point.usrCoords;
2578             }
2579 
2580             v = Mat.crossProduct(v, coords);
2581             return new Coords(Const.COORDS_BY_USER, Mat.crossProduct(v, line.stdform), board);
2582         },
2583 
2584         /**
2585          * Calculates the coordinates of the orthogonal projection of a given coordinate array on a given line
2586          * segment defined by two coordinate arrays.
2587          * @param {Array} p Point to project.
2588          * @param {Array} q1 Start point of the line segment on that the point is projected.
2589          * @param {Array} q2 End point of the line segment on that the point is projected.
2590          * @returns {Array} The coordinates of the projection of the given point on the given segment
2591          * and the factor that determines the projected point as a convex combination of the
2592          * two endpoints q1 and q2 of the segment.
2593          */
2594         projectCoordsToSegment: function (p, q1, q2) {
2595             var t, denom,
2596                 s = [q2[1] - q1[1], q2[2] - q1[2]],
2597                 v = [p[1] - q1[1], p[2] - q1[2]];
2598 
2599             /**
2600              * If the segment has length 0, i.e. is a point,
2601              * the projection is equal to that point.
2602              */
2603             if (Math.abs(s[0]) < Mat.eps && Math.abs(s[1]) < Mat.eps) {
2604                 return [q1, 0];
2605             }
2606 
2607             t = Mat.innerProduct(v, s);
2608             denom = Mat.innerProduct(s, s);
2609             t /= denom;
2610 
2611             return [ [1, t * s[0] + q1[1], t * s[1] + q1[2]], t];
2612         },
2613 
2614         /**
2615          * Finds the coordinates of the closest point on a Bezier segment of a
2616          * {@link JXG.Curve} to a given coordinate array.
2617          * @param {Array} pos Point to project in homogeneous coordinates.
2618          * @param {JXG.Curve} curve Curve of type "plot" having Bezier degree 3.
2619          * @param {Number} start Number of the Bezier segment of the curve.
2620          * @returns {Array} The coordinates of the projection of the given point
2621          * on the given Bezier segment and the preimage of the curve which
2622          * determines the closest point.
2623          */
2624         projectCoordsToBeziersegment: function (pos, curve, start) {
2625             var t0,
2626                 /** @ignore */
2627                 minfunc = function (t) {
2628                     var z = [1, curve.X(start + t), curve.Y(start + t)];
2629 
2630                     z[1] -= pos[1];
2631                     z[2] -= pos[2];
2632 
2633                     return z[1] * z[1] + z[2] * z[2];
2634                 };
2635 
2636             t0 = JXG.Math.Numerics.fminbr(minfunc, [0.0, 1.0]);
2637 
2638             return [[1, curve.X(t0 + start), curve.Y(t0 + start)], t0];
2639         },
2640 
2641         /**
2642          * Calculates the coordinates of the projection of a given point on a given curve.
2643          * Uses {@link JXG.Math.Geometry.projectCoordsToCurve}.
2644          *
2645          * @param {JXG.Point} point Point to project.
2646          * @param {JXG.Curve} curve Curve on that the point is projected.
2647          * @param {JXG.Board} [board=point.board] Reference to a board.
2648          * @see #projectCoordsToCurve
2649          * @returns {Array} [JXG.Coords, position] The coordinates of the projection of the given
2650          * point on the given graph and the relative position on the curve (real number).
2651          */
2652         projectPointToCurve: function (point, curve, board) {
2653             if (!Type.exists(board)) {
2654                 board = point.board;
2655             }
2656 
2657             var x = point.X(),
2658                 y = point.Y(),
2659                 t = point.position || 0.0,
2660                 result = this.projectCoordsToCurve(x, y, t, curve, board);
2661 
2662             // point.position = result[1];
2663 
2664             return result;
2665         },
2666 
2667         /**
2668          * Calculates the coordinates of the projection of a coordinates pair on a given curve. In case of
2669          * function graphs this is the
2670          * intersection point of the curve and the parallel to y-axis through the given point.
2671          * @param {Number} x coordinate to project.
2672          * @param {Number} y coordinate to project.
2673          * @param {Number} t start value for newtons method
2674          * @param {JXG.Curve} curve Curve on that the point is projected.
2675          * @param {JXG.Board} [board=curve.board] Reference to a board.
2676          * @see #projectPointToCurve
2677          * @returns {JXG.Coords} Array containing the coordinates of the projection of the given point on the given curve and
2678          * the position on the curve.
2679          */
2680         projectCoordsToCurve: function (x, y, t, curve, board) {
2681             var newCoords, newCoordsObj, i, j,
2682                 mindist, dist, lbda, v, coords, d,
2683                 p1, p2, res,
2684                 minfunc, t_new, f_new, f_old, delta, steps,
2685                 minX, maxX,
2686                 infty = Number.POSITIVE_INFINITY;
2687 
2688             if (!Type.exists(board)) {
2689                 board = curve.board;
2690             }
2691 
2692 
2693             if (Type.evaluate(curve.visProp.curvetype) === 'plot') {
2694                 t = 0;
2695                 mindist = infty;
2696                 if (curve.numberPoints === 0) {
2697                     newCoords = [0, 1, 1];
2698                 } else {
2699                     newCoords = [curve.Z(0), curve.X(0), curve.Y(0)];
2700                 }
2701 
2702                 if (curve.numberPoints > 1) {
2703                     v = [1, x, y];
2704                     if (curve.bezierDegree === 3) {
2705                         j = 0;
2706                     } else {
2707                         p1 = [curve.Z(0), curve.X(0), curve.Y(0)];
2708                     }
2709                     for (i = 0; i < curve.numberPoints - 1; i++) {
2710                         if (curve.bezierDegree === 3) {
2711                             res = this.projectCoordsToBeziersegment(v, curve, j);
2712                         } else {
2713                             p2 = [curve.Z(i + 1), curve.X(i + 1), curve.Y(i + 1)];
2714                             res = this.projectCoordsToSegment(v, p1, p2);
2715                         }
2716                         lbda = res[1];
2717                         coords = res[0];
2718 
2719                         if (0.0 <= lbda && lbda <= 1.0) {
2720                             dist = this.distance(coords, v);
2721                             d = i + lbda;
2722                         } else if (lbda < 0.0) {
2723                             coords = p1;
2724                             dist = this.distance(p1, v);
2725                             d = i;
2726                         } else if (lbda > 1.0 && i === curve.numberPoints - 2) {
2727                             coords = p2;
2728                             dist = this.distance(coords, v);
2729                             d = curve.numberPoints - 1;
2730                         }
2731 
2732                         if (dist < mindist) {
2733                             mindist = dist;
2734                             t = d;
2735                             newCoords = coords;
2736                         }
2737 
2738                         if (curve.bezierDegree === 3) {
2739                             j++;
2740                             i += 2;
2741                         } else {
2742                             p1 = p2;
2743                         }
2744                     }
2745                 }
2746 
2747                 newCoordsObj = new Coords(Const.COORDS_BY_USER, newCoords, board);
2748             } else {   // 'parameter', 'polar', 'functiongraph'
2749                 /** @ignore */
2750                 minfunc = function (t) {
2751                     var dx, dy;
2752                     if (t < curve.minX() || t > curve.maxX()) {
2753                         return Infinity;
2754                     }
2755                     dx = x - curve.X(t);
2756                     dy = y - curve.Y(t);
2757                     return dx * dx + dy * dy;
2758                 };
2759 
2760                 f_old = minfunc(t);
2761                 steps = 50;
2762                 minX = curve.minX();
2763                 maxX = curve.maxX();
2764 
2765                 delta = (maxX - minX) / steps;
2766                 t_new = minX;
2767 
2768                 for (i = 0; i < steps; i++) {
2769                     f_new = minfunc(t_new);
2770 
2771                     if (f_new < f_old || f_old === Infinity || isNaN(f_old)) {
2772                         t = t_new;
2773                         f_old = f_new;
2774                     }
2775 
2776                     t_new += delta;
2777                 }
2778 
2779                 //t = Numerics.root(Numerics.D(minfunc), t);
2780                 t = Numerics.fminbr(minfunc, [Math.max(t - delta, minX), Math.min(t + delta, maxX)]);
2781 
2782                 // Distinction between closed and open curves is not necessary.
2783                 // If closed, the cyclic projection shift will work anyhow
2784                 // if (Math.abs(curve.X(minX) - curve.X(maxX)) < Mat.eps &&
2785                 //     Math.abs(curve.Y(minX) - curve.Y(maxX)) < Mat.eps) {
2786                 //     // Cyclically
2787                 //     if (t < minX) {
2788                 //         t = maxX + t - minX;
2789                 //     }
2790                 //     if (t > maxX) {
2791                 //         t = minX + t - maxX;
2792                 //     }
2793                 // } else {
2794                     t = (t < minX) ? minX : t;
2795                     t = (t > maxX) ? maxX : t;
2796                 // }
2797 
2798                 newCoordsObj = new Coords(Const.COORDS_BY_USER, [curve.X(t), curve.Y(t)], board);
2799             }
2800 
2801             return [curve.updateTransform(newCoordsObj), t];
2802         },
2803 
2804         /**
2805          * Calculates the coordinates of the closest orthogonal projection of a given coordinate array onto the
2806          * border of a polygon.
2807          * @param {Array} p Point to project.
2808          * @param {JXG.Polygon} pol Polygon element
2809          * @returns {Array} The coordinates of the closest projection of the given point to the border of the polygon.
2810          */
2811         projectCoordsToPolygon: function (p, pol) {
2812             var i,
2813                 len = pol.vertices.length,
2814                 d_best = Infinity,
2815                 d,
2816                 projection, proj,
2817                 bestprojection;
2818 
2819             for (i = 0; i < len - 1; i++) {
2820                 projection = JXG.Math.Geometry.projectCoordsToSegment(
2821                     p,
2822                     pol.vertices[i].coords.usrCoords,
2823                     pol.vertices[i + 1].coords.usrCoords
2824                 );
2825 
2826                 if (0 <= projection[1] && projection[1] <= 1) {
2827                     d = JXG.Math.Geometry.distance(projection[0], p, 3);
2828                     proj = projection[0];
2829                 } else if (projection[1] < 0) {
2830                     d = JXG.Math.Geometry.distance(pol.vertices[i].coords.usrCoords, p, 3);
2831                     proj = pol.vertices[i].coords.usrCoords;
2832                 } else {
2833                     d = JXG.Math.Geometry.distance(pol.vertices[i + 1].coords.usrCoords, p, 3);
2834                     proj = pol.vertices[i + 1].coords.usrCoords;
2835                 }
2836                 if (d < d_best) {
2837                     bestprojection = proj.slice(0);
2838                     d_best = d;
2839                 }
2840             }
2841             return bestprojection;
2842         },
2843 
2844         /**
2845          * Calculates the coordinates of the projection of a given point on a given turtle. A turtle consists of
2846          * one or more curves of curveType 'plot'. Uses {@link JXG.Math.Geometry.projectPointToCurve}.
2847          * @param {JXG.Point} point Point to project.
2848          * @param {JXG.Turtle} turtle on that the point is projected.
2849          * @param {JXG.Board} [board=point.board] Reference to a board.
2850          * @returns {Array} [JXG.Coords, position] Array containing the coordinates of the projection of the given point on the turtle and
2851          * the position on the turtle.
2852          */
2853         projectPointToTurtle: function (point, turtle, board) {
2854             var newCoords, t, x, y, i, dist, el, minEl,
2855                 res, newPos,
2856                 np = 0,
2857                 npmin = 0,
2858                 mindist = Number.POSITIVE_INFINITY,
2859                 len = turtle.objects.length;
2860 
2861             if (!Type.exists(board)) {
2862                 board = point.board;
2863             }
2864 
2865             // run through all curves of this turtle
2866             for (i = 0; i < len; i++) {
2867                 el = turtle.objects[i];
2868 
2869                 if (el.elementClass === Const.OBJECT_CLASS_CURVE) {
2870                     res = this.projectPointToCurve(point, el);
2871                     newCoords = res[0];
2872                     newPos = res[1];
2873                     dist = this.distance(newCoords.usrCoords, point.coords.usrCoords);
2874 
2875                     if (dist < mindist) {
2876                         x = newCoords.usrCoords[1];
2877                         y = newCoords.usrCoords[2];
2878                         t = newPos;
2879                         mindist = dist;
2880                         minEl = el;
2881                         npmin = np;
2882                     }
2883                     np += el.numberPoints;
2884                 }
2885             }
2886 
2887             newCoords = new Coords(Const.COORDS_BY_USER, [x, y], board);
2888             // point.position = t + npmin;
2889             // return minEl.updateTransform(newCoords);
2890             return [minEl.updateTransform(newCoords), t + npmin];
2891         },
2892 
2893         /**
2894          * Trivial projection of a point to another point.
2895          * @param {JXG.Point} point Point to project (not used).
2896          * @param {JXG.Point} dest Point on that the point is projected.
2897          * @returns {JXG.Coords} The coordinates of the projection of the given point on the given circle.
2898          */
2899         projectPointToPoint: function (point, dest) {
2900             return dest.coords;
2901         },
2902 
2903         /**
2904          *
2905          * @param {JXG.Point|JXG.Coords} point
2906          * @param {JXG.Board} [board]
2907          */
2908         projectPointToBoard: function (point, board) {
2909             var i, l, c,
2910                 brd = board || point.board,
2911                 // comparison factor, point coord idx, bbox idx, 1st bbox corner x & y idx, 2nd bbox corner x & y idx
2912                 config = [
2913                     // left
2914                     [1, 1, 0, 0, 3, 0, 1],
2915                     // top
2916                     [-1, 2, 1, 0, 1, 2, 1],
2917                     // right
2918                     [-1, 1, 2, 2, 1, 2, 3],
2919                     // bottom
2920                     [1, 2, 3, 0, 3, 2, 3]
2921                 ],
2922                 coords = point.coords || point,
2923                 bbox = brd.getBoundingBox();
2924 
2925             for (i = 0; i < 4; i++) {
2926                 c = config[i];
2927                 if (c[0] * coords.usrCoords[c[1]] < c[0] * bbox[c[2]]) {
2928                     // define border
2929                     l = Mat.crossProduct([1, bbox[c[3]], bbox[c[4]]], [1, bbox[c[5]], bbox[c[6]]]);
2930                     l[3] = 0;
2931                     l = Mat.normalize(l);
2932 
2933                     // project point
2934                     coords = this.projectPointToLine({coords: coords}, {stdform: l}, brd);
2935                 }
2936             }
2937 
2938             return coords;
2939         },
2940 
2941         /**
2942          * Calculates the distance of a point to a line. The point and the line are given by homogeneous
2943          * coordinates. For lines this can be line.stdform.
2944          * @param {Array} point Homogeneous coordinates of a point.
2945          * @param {Array} line Homogeneous coordinates of a line ([C,A,B] where A*x+B*y+C*z=0).
2946          * @returns {Number} Distance of the point to the line.
2947          */
2948         distPointLine: function (point, line) {
2949             var a = line[1],
2950                 b = line[2],
2951                 c = line[0],
2952                 nom;
2953 
2954             if (Math.abs(a) + Math.abs(b) < Mat.eps) {
2955                 return Number.POSITIVE_INFINITY;
2956             }
2957 
2958             nom = a * point[1] + b * point[2] + c;
2959             a *= a;
2960             b *= b;
2961 
2962             return Math.abs(nom) / Math.sqrt(a + b);
2963         },
2964 
2965         /**
2966          * Helper function to create curve which displays a Reuleaux polygons.
2967          * @param {Array} points Array of points which should be the vertices of the Reuleaux polygon. Typically,
2968          * these point list is the array vertices of a regular polygon.
2969          * @param {Number} nr Number of vertices
2970          * @returns {Array} An array containing the two functions defining the Reuleaux polygon and the two values
2971          * for the start and the end of the paramtric curve. array may be used as parent array of a
2972          * {@link JXG.Curve}.
2973          *
2974          * @example
2975          * var A = brd.create('point',[-2,-2]);
2976          * var B = brd.create('point',[0,1]);
2977          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0});
2978          * var reuleauxTriangle = brd.create('curve', JXG.Math.Geometry.reuleauxPolygon(pol.vertices, 3),
2979          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
2980          *
2981          * </pre><div class="jxgbox" id="JXG2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div>
2982          * <script type="text/javascript">
2983          * var brd = JXG.JSXGraph.initBoard('JXG2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false});
2984          * var A = brd.create('point',[-2,-2]);
2985          * var B = brd.create('point',[0,1]);
2986          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0});
2987          * var reuleauxTriangle = brd.create('curve', JXG.Math.Geometry.reuleauxPolygon(pol.vertices, 3),
2988          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
2989          * </script><pre>
2990          */
2991         reuleauxPolygon: function (points, nr) {
2992             var beta,
2993                 pi2 = Math.PI * 2,
2994                 pi2_n = pi2 / nr,
2995                 diag = (nr - 1) / 2,
2996                 d = 0,
2997                 makeFct = function (which, trig) {
2998                     return function (t, suspendUpdate) {
2999                         var t1 = (t % pi2 + pi2) % pi2,
3000                             j = Math.floor(t1 / pi2_n) % nr;
3001 
3002                         if (!suspendUpdate) {
3003                             d = points[0].Dist(points[diag]);
3004                             beta = Mat.Geometry.rad([points[0].X() + 1, points[0].Y()], points[0], points[diag % nr]);
3005                         }
3006 
3007                         if (isNaN(j)) {
3008                             return j;
3009                         }
3010 
3011                         t1 = t1 * 0.5 + j * pi2_n * 0.5 + beta;
3012 
3013                         return points[j][which]() + d * Math[trig](t1);
3014                     };
3015                 };
3016 
3017             return [makeFct('X', 'cos'), makeFct('Y', 'sin'), 0, pi2];
3018         },
3019 
3020 
3021         meet3Planes: function (n1, d1, n2, d2, n3, d3) {
3022             var p = [0, 0, 0],
3023                 n31, n12, n23, denom,
3024                 i;
3025 
3026             n31 = Mat.crossProduct(n3, n1);
3027             n12 = Mat.crossProduct(n1, n2);
3028             n23 = Mat.crossProduct(n2, n3);
3029             denom = Mat.innerProduct(n1, n23, 3);
3030             for (i = 0; i < 3; i++) {
3031                 p[i] = (d1 * n23[i] + d2 * n31[i] + d3 * n12[i]) / denom;
3032             }
3033             return p;
3034         },
3035 
3036         meetPlanePlane: function (v11, v12, v21, v22) {
3037             var i, no1, no2,
3038                 v = [0, 0, 0],
3039                 w = [0, 0, 0];
3040 
3041             for (i = 0; i < 3; i++) {
3042                 v[i] = Type.evaluate(v11[i]);
3043                 w[i] = Type.evaluate(v12[i]);
3044             }
3045             no1 = Mat.crossProduct(v, w);
3046 
3047             for (i = 0; i < 3; i++) {
3048                 v[i] = Type.evaluate(v21[i]);
3049                 w[i] = Type.evaluate(v22[i]);
3050             }
3051             no2 = Mat.crossProduct(v, w);
3052 
3053             return Mat.crossProduct(no1, no2);
3054         },
3055 
3056         project3DTo3DPlane: function (point, normal, foot) {
3057             // TODO: homogeneous 3D coordinates
3058             var sol = [0, 0, 0],
3059                 le, d1, d2, lbda;
3060 
3061             foot = foot || [0, 0, 0];
3062 
3063             le = Mat.norm(normal);
3064             d1 = Mat.innerProduct(point, normal, 3);
3065             d2 = Mat.innerProduct(foot, normal, 3);
3066             // (point - lbda * normal / le) * normal / le == foot * normal / le
3067             // => (point * normal - foot * normal) ==  lbda * le
3068             lbda = (d1 - d2) / le;
3069             sol = Mat.axpy(-lbda, normal, point);
3070 
3071             return sol;
3072         },
3073 
3074         getPlaneBounds: function (v1, v2, q, s, e) {
3075             var s1, s2, e1, e2, mat, rhs, sol;
3076 
3077             if (v1[2] + v2[0] !== 0) {
3078                 mat = [
3079                     [v1[0], v2[0]],
3080                     [v1[1], v2[1]]
3081                 ];
3082                 rhs = [s - q[0], s - q[1]];
3083 
3084                 sol = Numerics.Gauss(mat, rhs);
3085                 s1 = sol[0];
3086                 s2 = sol[1];
3087 
3088                 rhs = [e - q[0], e - q[1]];
3089                 sol = Numerics.Gauss(mat, rhs);
3090                 e1 = sol[0];
3091                 e2 = sol[1];
3092                 return [s1, e1, s2, e2];
3093             }
3094             return null;
3095         }
3096 
3097     });
3098 
3099     return Mat.Geometry;
3100 });
3101