1 /*
  2     Copyright 2008,2009
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Bianca Valentin,
  7         Alfred Wassermann,
  8         Peter Wilfahrt
  9 
 10     This file is part of JSXGraph.
 11 
 12     JSXGraph is free software: you can redistribute it and/or modify
 13     it under the terms of the GNU Lesser General Public License as published by
 14     the Free Software Foundation, either version 3 of the License, or
 15     (at your option) any later version.
 16 
 17     JSXGraph is distributed in the hope that it will be useful,
 18     but WITHOUT ANY WARRANTY; without even the implied warranty of
 19     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 20     GNU Lesser General Public License for more details.
 21 
 22     You should have received a copy of the GNU Lesser General Public License
 23     along with JSXGraph. If not, see <http://www.gnu.org/licenses/>.
 24 */
 25 
 26 /** 
 27  * @fileoverview This file contains the Math.Geometry namespace for calculating algebraic/geometric
 28  * stuff like intersection points, angles, midpoint, and so on.
 29  */
 30 
 31 
 32 /**
 33  * Math.Geometry namespace definition
 34  */
 35 JXG.Math.Geometry = {
 36 
 37     /****************************************/
 38     /**** GENERAL GEOMETRIC CALCULATIONS ****/
 39     /****************************************/
 40 
 41     /**
 42      * Calculates the angle defined by the points A, B, C.
 43      * @param {JXG.Point} A A point  or [x,y] array.
 44      * @param {JXG.Point} B Another point or [x,y] array.
 45      * @param {JXG.Point} C A circle - no, of course the third point or [x,y] array.
 46      * @deprecated Use {@link JXG.Math.Geometry#rad} instead.
 47      * @see #rad
 48      * @see #trueAngle
 49      * @returns {Number} The angle in radian measure.
 50      */
 51     angle: function(A, B, C) {
 52         var a = [],
 53             b = [],
 54             c = [],
 55             u, v, s, t;
 56 
 57         if (A.coords == null) {
 58             a[0] = A[0];
 59             a[1] = A[1];
 60         } else {
 61             a[0] = A.coords.usrCoords[1];
 62             a[1] = A.coords.usrCoords[2];
 63         }
 64         if (B.coords == null) {
 65             b[0] = B[0];
 66             b[1] = B[1];
 67         } else {
 68             b[0] = B.coords.usrCoords[1];
 69             b[1] = B.coords.usrCoords[2];
 70         }
 71         if (C.coords == null) {
 72             c[0] = C[0];
 73             c[1] = C[1];
 74         } else {
 75             c[0] = C.coords.usrCoords[1];
 76             c[1] = C.coords.usrCoords[2];
 77         }
 78         u = a[0] - b[0];
 79         v = a[1] - b[1];
 80         s = c[0] - b[0];
 81         t = c[1] - b[1];
 82         return Math.atan2(u * t - v * s, u * s + v * t);
 83     },
 84 
 85     /**
 86      * Calculates the angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
 87      * @param {JXG.Point} A Point or [x,y] array
 88      * @param {JXG.Point} B Point or [x,y] array
 89      * @param {JXG.Point} C Point or [x,y] array
 90      * @see #rad
 91      * @returns {Number} The angle in degrees.
 92      */
 93     trueAngle: function(A, B, C) {
 94         return this.rad(A, B, C) * 57.295779513082323; // *180.0/Math.PI;
 95     },
 96 
 97     /**
 98      * Calculates the internal angle defined by the three points A, B, C if you're going from A to C around B counterclockwise.
 99      * @param {JXG.Point} A Point or [x,y] array
100      * @param {JXG.Point} B Point or [x,y] array
101      * @param {JXG.Point} C Point or [x,y] array
102      * @see #trueAngle
103      * @returns {Number} Angle in radians.
104      */
105     rad: function(A, B, C) {
106         var ax, ay, bx, by, cx, cy,
107             phi;
108 
109         if (A.coords == null) {
110             ax = A[0];
111             ay = A[1];
112         } else {
113             ax = A.coords.usrCoords[1];
114             ay = A.coords.usrCoords[2];
115         }
116         if (B.coords == null) {
117             bx = B[0];
118             by = B[1];
119         } else {
120             bx = B.coords.usrCoords[1];
121             by = B.coords.usrCoords[2];
122         }
123         if (C.coords == null) {
124             cx = C[0];
125             cy = C[1];
126         } else {
127             cx = C.coords.usrCoords[1];
128             cy = C.coords.usrCoords[2];
129         }
130 
131         phi = Math.atan2(cy - by, cx - bx) - Math.atan2(ay - by, ax - bx);
132         if (phi < 0) phi += 6.2831853071795862;
133         return phi;
134     },
135 
136     /**
137      * Calculates the bisection between the three points A, B, C. The bisection is defined by two points:
138      * Parameter B and a point with the coordinates calculated in this function.
139      * @param {JXG.Point} A Point
140      * @param {JXG.Point} B Point
141      * @param {JXG.Point} C Point
142      * @param [board=A.board] Reference to the board
143      * @returns {JXG.Coords} Coordinates of the second point defining the bisection.
144      */
145     angleBisector: function(A, B, C, board) {
146         /* First point */
147         var Ac = A.coords.usrCoords,
148             Bc = B.coords.usrCoords,
149             Cc = C.coords.usrCoords,
150             x = Ac[1] - Bc[1],
151             y = Ac[2] - Bc[2],
152             d = Math.sqrt(x * x + y * y),
153             phiA, phiC, phi;
154 
155         if (!JXG.exists(board))
156             board = A.board;
157 
158         x /= d;
159         y /= d;
160 
161         phiA = Math.acos(x);
162         if (y < 0) {
163             phiA *= -1;
164         }
165         if (phiA < 0) {
166             phiA += 2 * Math.PI;
167         }
168 
169         /* Second point */
170         x = Cc[1] - Bc[1];
171         y = Cc[2] - Bc[2];
172         d = Math.sqrt(x * x + y * y);
173         x /= d;
174         y /= d;
175 
176         phiC = Math.acos(x);
177         if (y < 0) {
178             phiC *= -1;
179         }
180         if (phiC < 0) {
181             phiC += 2 * Math.PI;
182         }
183 
184         phi = (phiA + phiC) * 0.5;
185         if (phiA > phiC) {
186             phi += Math.PI;
187         }
188 
189         x = Math.cos(phi) + Bc[1];
190         y = Math.sin(phi) + Bc[2];
191 
192         return new JXG.Coords(JXG.COORDS_BY_USER, [x,y], board);
193     },
194 
195     /**
196      * Reflects the point along the line.
197      * @param {JXG.Line} line Axis of reflection.
198      * @param {JXG.Point} point Point to reflect.
199      * @param [board=point.board] Reference to the board
200      * @returns {JXG.Coords} Coordinates of the reflected point.
201      */
202     reflection: function(line, point, board) {
203         /* (v,w) defines the slope of the line */
204         var pc = point.coords.usrCoords,
205             p1c = line.point1.coords.usrCoords,
206             p2c = line.point2.coords.usrCoords,
207             x0, y0, x1, y1, v, w, mu;
208 
209         if (!JXG.exists(board))
210             board = point.board;
211 
212         v = p2c[1] - p1c[1];
213         w = p2c[2] - p1c[2];
214 
215         x0 = pc[1] - p1c[1];
216         y0 = pc[2] - p1c[2];
217 
218         mu = (v * y0 - w * x0) / (v * v + w * w);
219 
220         /* point + mu*(-y,x) waere Lotpunkt */
221         x1 = pc[1] + 2 * mu * w;
222         y1 = pc[2] - 2 * mu * v;
223 
224         return new JXG.Coords(JXG.COORDS_BY_USER, [x1,y1], board);
225     },
226 
227     /**
228      * Computes the new position of a point which is rotated
229      * around a second point (called rotpoint) by the angle phi.
230      * @param {JXG.Point} rotpoint Center of the rotation
231      * @param {JXG.Point} point point to be rotated
232      * @param {number} phi rotation angle in arc length
233      * @param {JXG.Board} [board=point.board] Reference to the board
234      * @returns {JXG.Coords} Coordinates of the new position.
235      */
236     rotation: function(rotpoint, point, phi, board) {
237         var pc = point.coords.usrCoords,
238             rotpc = rotpoint.coords.usrCoords,
239             x0, y0, c, s, x1, y1;
240 
241         if (!JXG.exists(board))
242             board = point.board;
243 
244         x0 = pc[1] - rotpc[1];
245         y0 = pc[2] - rotpc[2];
246 
247         c = Math.cos(phi);
248         s = Math.sin(phi);
249 
250         x1 = x0 * c - y0 * s + rotpc[1];
251         y1 = x0 * s + y0 * c + rotpc[2];
252 
253         return new JXG.Coords(JXG.COORDS_BY_USER, [x1,y1], board);
254     },
255 
256     /**
257      * Calculates the coordinates of a point on the perpendicular to the given line through
258      * the given point.
259      * @param {JXG.Line} line A line.
260      * @param {JXG.Point} point Intersection point of line to perpendicular.
261      * @param {JXG.Board} [board=point.board] Reference to the board
262      * @returns {JXG.Coords} Coordinates of a point on the perpendicular to the given line through the given point.
263      */
264     perpendicular: function(line, point, board) {
265         var A = line.point1.coords.usrCoords,
266             B = line.point2.coords.usrCoords,
267             C = point.coords.usrCoords,
268             x, y, change,
269             fmd, emc, d0, d1, den;
270 
271         if (!JXG.exists(board))
272             board = point.board;
273 
274         if (point == line.point1) { // Punkt ist erster Punkt der Linie
275             x = A[1] + B[2] - A[2];
276             y = A[2] - B[1] + A[1];
277             change = true;
278         } else if (point == line.point2) {  // Punkt ist zweiter Punkt der Linie
279             x = B[1] + A[2] - B[2];
280             y = B[2] - A[1] + B[1];
281             change = false;
282         } else if (((Math.abs(A[1] - B[1]) > JXG.Math.eps) &&
283             (Math.abs(C[2] - (A[2] - B[2]) * (C[1] - A[1]) / (A[1] - B[1]) - A[2]) < JXG.Math.eps)) ||
284             ((Math.abs(A[1] - B[1]) <= JXG.Math.eps) && (Math.abs(A[1] - C[1]) < JXG.Math.eps))) { // Punkt liegt auf der Linie
285             x = C[1] + B[2] - C[2];
286             y = C[2] - B[1] + C[1];
287             change = true;
288             if (Math.abs(x - C[1]) < JXG.Math.eps && Math.abs(y - C[2]) < JXG.Math.eps) {
289                 x = C[1] + A[2] - C[2];
290                 y = C[2] - A[1] + C[1];
291                 change = false;
292             }
293         } else { // Punkt liegt nicht auf der Linie -> als zweiter Punkt wird der Lotfusspunkt gewaehlt
294             fmd = A[2] - B[2];
295             emc = A[1] - B[1];
296             d0 = B[1] * fmd - B[2] * emc;
297             d1 = C[1] * emc + C[2] * fmd;
298             den = fmd * fmd + emc * emc;
299             if (Math.abs(den) < JXG.Math.eps) {
300                 den = JXG.Math.eps;
301             }
302             x = (d0 * fmd + d1 * emc) / den;
303             y = (d1 * fmd - d0 * emc) / den;
304             change = true;
305         }
306         return [new JXG.Coords(JXG.COORDS_BY_USER, [x, y], board),change];
307     },
308 
309     /**
310      * Calculates the midpoint of the circumcircle of the three given points.
311      * @param {JXG.Point} point1 Point
312      * @param {JXG.Point} point2 Point
313      * @param {JXG.Point} point3 Point
314      * @param {JXG.Board} [board=point1.board] Reference to the board
315      * @returns {JXG.Coords} Coordinates of the midpoint of the circumcircle of the given points.
316      */
317     circumcenterMidpoint: function(point1, point2, point3, board) {
318         var A = point1.coords.usrCoords,
319             B = point2.coords.usrCoords,
320             C = point3.coords.usrCoords,
321             u, v, den, x, y;
322 
323         if (!JXG.exists(board))
324             board = point1.board;
325 
326         u = ((A[1] - B[1]) * (A[1] + B[1]) + (A[2] - B[2]) * (A[2] + B[2])) * 0.5;
327         v = ((B[1] - C[1]) * (B[1] + C[1]) + (B[2] - C[2]) * (B[2] + C[2])) * 0.5;
328         den = (A[1] - B[1]) * (B[2] - C[2]) - (B[1] - C[1]) * (A[2] - B[2]);
329 
330         if (Math.abs(den) < JXG.Math.eps) {
331             den = JXG.Math.eps;
332         }
333 
334         x = (u * (B[2] - C[2]) - v * (A[2] - B[2])) / den;
335         y = (v * (A[1] - B[1]) - u * (B[1] - C[1])) / den;
336 
337         return new JXG.Coords(JXG.COORDS_BY_USER, [x, y], board);
338     },
339 
340     /**
341      * Calculates euclidean norm for two given arrays of the same length.
342      * @param {Array} array1 Array of float or integer.
343      * @param {Array} array2 Array of float or integer.
344      * @returns {Number} Euclidean distance of the given vectors.
345      */
346     distance: function(array1, array2) {
347         var sum = 0,
348             i, len;
349 
350         if (array1.length != array2.length) {
351             return NaN;
352         }
353         len = array1.length;
354         for (i = 0; i < len; i++) {
355             sum += (array1[i] - array2[i]) * (array1[i] - array2[i]);
356         }
357         return Math.sqrt(sum);
358     },
359 
360     /**
361      * Calculates euclidean distance for two given arrays of the same length.
362      * If one of the arrays contains a zero in the first coordinate, and the euclidean distance
363      * is different from zero it is a point at infinity and we return Infinity.
364      * @param {Array} array1 Array containing elements of number.
365      * @param {Array} array2 Array containing elements of type number.
366      * @returns {Number} Euclidean (affine) distance of the given vectors.
367      */
368     affineDistance: function(array1, array2) {
369         var d;
370         if (array1.length != array2.length) {
371             return NaN;
372         }
373         d = this.distance(array1, array2);
374         if (d > JXG.Math.eps && (Math.abs(array1[0]) < JXG.Math.eps || Math.abs(array2[0]) < JXG.Math.eps)) {
375             return Infinity;
376         } else {
377             return d;
378         }
379     },
380 
381 
382     /****************************************/
383     /****          INTERSECTIONS         ****/
384     /****************************************/
385 
386     /**
387      * Calculates the coordinates of the intersection of the given lines.
388      * @param {JXG.Line} line1 Line.
389      * @param {JXG.Line} line2 Line.
390      * @param {JXG.Board} [board=line1.board] Reference to the board
391      * @returns {JXG.Coords} Coordinates of the intersection point of the given lines.
392      */
393     intersectLineLine: function(line1, line2, board) {
394         var A = line1.point1.coords.usrCoords,
395             B = line1.point2.coords.usrCoords,
396             C = line2.point1.coords.usrCoords,
397             D = line2.point2.coords.usrCoords,
398             d0, d1, den, x, y;
399 
400         if (!JXG.exists(board))
401             board = line1.board;
402 
403         d0 = A[1] * B[2] - A[2] * B[1];
404         d1 = C[1] * D[2] - C[2] * D[1];
405         den = (B[2] - A[2]) * (C[1] - D[1]) - (A[1] - B[1]) * (D[2] - C[2]);
406 
407         if (Math.abs(den) < JXG.Math.eps) {
408             den = JXG.Math.eps;
409         }
410         x = (d0 * (C[1] - D[1]) - d1 * (A[1] - B[1])) / den;
411         y = (d1 * (B[2] - A[2]) - d0 * (D[2] - C[2])) / den;
412 
413         return new JXG.Coords(JXG.COORDS_BY_USER, [x, y], board);
414     },
415 
416     /**
417      * Calculates the coordinates of the intersection of the given line and circle.
418      * @param {JXG.Circle} circle Circle.
419      * @param {JXG.Line} line Line.
420      * @param {JXG.Board} [board=line.board] Reference to the board
421      * @returns {Array} The coordinates of the intersection points of the given circle with the given line and
422      * the amount of intersection points in the first component of the array.
423      */
424     intersectCircleLine: function(circle, line, board) {
425         var eA = line.point1.coords.usrCoords,
426             eB = line.point2.coords.usrCoords,
427             fM = circle.midpoint.coords.usrCoords,
428             s, d0, d1, b, w, h, r, n1, dx, dy, firstPointX, firstPointY, l, x, y, n1s, firstPoint, secondPoint, d;
429 
430         if (!JXG.exists(board))
431             board = line.board;
432 
433         s = line.point1.Dist(line.point2);
434         if (s > 0) {
435             d0 = circle.midpoint.Dist(line.point1);
436             d1 = circle.midpoint.Dist(line.point2);
437             b = ((d0 * d0) + (s * s) - (d1 * d1)) / (2 * s);
438             w = (d0 * d0) - (b * b);
439             w = (w < 0) ? 0 : w;
440             h = Math.sqrt(w);
441 
442             r = circle.Radius();
443             n1 = Math.sqrt((r * r) - h * h);
444             dx = eB[1] - eA[1];
445             dy = eB[2] - eA[2];
446             firstPointX = fM[1] + (h / s) * dy;
447             firstPointY = fM[2] - (h / s) * dx;
448             d0 = (eB[1] * dy) - (eB[2] * dx);
449             d1 = (firstPointX * dx) + (firstPointY * dy);
450             l = (dy * dy) + (dx * dx);
451             if (Math.abs(l) < JXG.Math.eps) {
452                 l = JXG.Math.eps;
453             }
454             x = ((d0 * dy) + (d1 * dx)) / l;
455             y = ((d1 * dy) - (d0 * dx)) / l;
456             n1s = n1 / s;
457             firstPoint = new JXG.Coords(JXG.COORDS_BY_USER, [x + n1s * dx, y + n1s * dy], board);
458             secondPoint = new JXG.Coords(JXG.COORDS_BY_USER, [x - n1s * dx, y - n1s * dy], board);
459             d = circle.midpoint.coords.distance(JXG.COORDS_BY_USER, firstPoint);
460 
461             if ((r < (d - 1)) || isNaN(d)) {
462                 return [0];
463             } else {
464                 return [2,firstPoint,secondPoint];
465             }
466         }
467         return [0];
468     },
469 
470     /**
471      * Calculates the coordinates of the intersection of the given circles.
472      * @param {JXG.Circle} circle1 Circle.
473      * @param {JXG.Circle} circle2 Circle.
474      * @param {JXG.Board} [board=circle1.board] Reference to the board
475      * @returns {Array} Coordinates of the intersection points of the given circles and the
476      * amount of intersection points in the first component of the array.
477      */
478     intersectCircleCircle: function(circle1, circle2, board) {
479         var intersection = {},
480             r1 = circle1.Radius(),
481             r2 = circle2.Radius(),
482             M1 = circle1.midpoint.coords.usrCoords,
483             M2 = circle2.midpoint.coords.usrCoords,
484             rSum, rDiff, s,
485             dx, dy, a, h;
486 
487         if (!JXG.exists(board))
488             board = circle1.board;
489 
490         rSum = r1 + r2;
491         rDiff = Math.abs(r1 - r2);
492         // Abstand der Mittelpunkte der beiden Kreise
493         s = circle1.midpoint.coords.distance(JXG.COORDS_BY_USER, circle2.midpoint.coords);
494         if (s > rSum) {
495             return [0]; // Kreise schneiden sich nicht, liegen nebeneinander
496         } else if (s < rDiff) {
497             return [0]; // Kreise schneiden sich nicht, liegen ineinander
498         } else {
499             if (s != 0) {
500                 intersection[0] = 1; // es gibt einen Schnitt
501                 dx = M2[1] - M1[1];
502                 dy = M2[2] - M1[2];
503                 a = (s * s - r2 * r2 + r1 * r1) / (2 * s);
504                 h = Math.sqrt(r1 * r1 - a * a);
505                 intersection[1] = new JXG.Coords(JXG.COORDS_BY_USER,
506                     [M1[1] + (a / s) * dx + (h / s) * dy,
507                         M1[2] + (a / s) * dy - (h / s) * dx],
508                     board);
509                 intersection[2] = new JXG.Coords(JXG.COORDS_BY_USER,
510                     [M1[1] + (a / s) * dx - (h / s) * dy,
511                         M1[2] + (a / s) * dy + (h / s) * dx],
512                     board);
513             } else {
514                 return [0]; // vorsichtshalber...
515             }
516             return intersection;
517         }
518     },
519 
520     /**
521      * Computes the intersection of a pair of lines, circles or both.
522      * It uses the internal data array stdform of these elements.
523      * @param {Array} el1 stdform of the first element (line or circle)
524      * @param {Array} el2 stdform of the second element (line or circle)
525      * @param {Number} i Index of the intersection point that should be returned.
526      * @param board Reference to the board.
527      * @returns {JXG.Coords} Coordinates of one of the possible two or more intersection points.
528      * Which point will be returned is determined by i.
529      */
530     meet: function(el1, el2, i, board) {
531         var eps = JXG.Math.eps; //    var eps = 0.000001;
532 
533         if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) < eps) { // line line
534             return this.meetLineLine(el1, el2, i, board);
535         } else if (Math.abs(el1[3]) >= eps && Math.abs(el2[3]) < eps) { // circle line
536             return this.meetLineCircle(el2, el1, i, board);
537         } else if (Math.abs(el1[3]) < eps && Math.abs(el2[3]) >= eps) { // line circle
538             return this.meetLineCircle(el1, el2, i, board);
539         } else {  // circle circle
540             return this.meetCircleCircle(el1, el2, i, board);
541         }
542     },
543 
544     /**
545      * Intersection of two lines using the stdform.
546      * @param {Array} l1 stdform of the first line
547      * @param {Array} l2 stdform of the second line
548      * @param {number} i unused
549      * @param {JXG.Board} board Reference to the board.
550      * @returns {JXG.Coords} Coordinates of the intersection point.
551      */
552     meetLineLine: function(l1, l2, i, board) {
553         var s = JXG.Math.crossProduct(l1, l2);
554         if (Math.abs(s[0]) > JXG.Math.eps) {
555             s[1] /= s[0];
556             s[2] /= s[0];
557             s[0] = 1.0;
558         }
559         return new JXG.Coords(JXG.COORDS_BY_USER, s, board);
560     },
561 
562     /**
563      * Intersection of line and circle using the stdform.
564      * @param {Array} lin stdform of the line
565      * @param {Array} circ stdform of the circle
566      * @param {number} i number of the returned intersection point.
567      *   i==0: use the positive square root,
568      *   i==1: use the negative square root.
569      * @param {JXG.Board} board Reference to a board.
570      * @returns {JXG.Coords} Coordinates of the intersection point
571      */
572     meetLineCircle: function(lin, circ, i, board) {
573         var a,b,c,d,n, A,B,C, k,t;
574 
575         if (circ[4] < JXG.Math.eps) { // Radius is zero, return center of circle
576             return new JXG.Coords(JXG.COORDS_BY_USER, circ.slice(1, 3), board);
577         }
578         c = circ[0];
579         b = circ.slice(1, 3);
580         a = circ[3];
581         d = lin[0];
582         n = lin.slice(1, 3);
583 
584         // Line is normalized, therefore nn==1 and we can skip some operations:
585         /*
586          var nn = n[0]*n[0]+n[1]*n[1];
587          A = a*nn;
588          B = (b[0]*n[1]-b[1]*n[0])*nn;
589          C = a*d*d - (b[0]*n[0]+b[1]*n[1])*d + c*nn;
590          */
591         A = a;
592         B = (b[0] * n[1] - b[1] * n[0]);
593         C = a * d * d - (b[0] * n[0] + b[1] * n[1]) * d + c;
594 
595         k = B * B - 4 * A * C;
596         if (k >= 0) {
597             k = Math.sqrt(k);
598             t = [(-B + k) / (2 * A),(-B - k) / (2 * A)];
599             return ((i == 0)
600                 ? new JXG.Coords(JXG.COORDS_BY_USER, [-t[0] * (-n[1]) - d * n[0],-t[0] * n[0] - d * n[1]], board)
601                 : new JXG.Coords(JXG.COORDS_BY_USER, [-t[1] * (-n[1]) - d * n[0],-t[1] * n[0] - d * n[1]], board)
602                 );
603             /*
604              new JXG.Coords(JXG.COORDS_BY_USER, [-t[0]*(-n[1])-d*n[0]/nn,-t[0]*n[0]-d*n[1]/nn], this.board),
605              new JXG.Coords(JXG.COORDS_BY_USER, [-t[1]*(-n[1])-d*n[0]/nn,-t[1]*n[0]-d*n[1]/nn], this.board)
606              */
607         } else {
608             return new JXG.Coords(JXG.COORDS_BY_USER, [NaN,NaN], board);
609         }
610         // Returns do not work with homogeneous coordinates, yet
611     },
612 
613     /**
614      * Intersection of two circles using the stdform.
615      * @param {Array} circ1 stdform of the first circle
616      * @param {Array} circ2 stdform of the second circle
617      * @param {number} i number of the returned intersection point.
618      *   i==0: use the positive square root,
619      *   i==1: use the negative square root.
620      * @param {JXG.Board} board Reference to the board.
621      * @returns {JXG.Coords} Coordinates of the intersection point
622      */
623     meetCircleCircle: function(circ1, circ2, i, board) {
624         var radicalAxis;
625         if (circ1[4] < JXG.Math.eps) { // Radius are zero, return center of circle, if on other circle
626             if (this.distance(circ1.slice(1, 3), circ2.slice(1, 3)) == circ2[4]) {
627                 return new JXG.Coords(JXG.COORDS_BY_USER, circ1.slice(1, 3), board);
628             } else {
629                 return new JXG.Coords(JXG.COORDS_BY_USER, [NaN,NaN], board);
630             }
631         }
632         if (circ2[4] < JXG.Math.eps) { // Radius are zero, return center of circle, if on other circle
633             if (this.distance(circ2.slice(1, 3), circ1.slice(1, 3)) == circ1[4]) {
634                 return new JXG.Coords(JXG.COORDS_BY_USER, circ2.slice(1, 3), board);
635             } else {
636                 return new JXG.Coords(JXG.COORDS_BY_USER, [NaN,NaN], board);
637             }
638         }
639         radicalAxis = [circ2[3] * circ1[0] - circ1[3] * circ2[0],
640             circ2[3] * circ1[1] - circ1[3] * circ2[1],
641             circ2[3] * circ1[2] - circ1[3] * circ2[2],
642             0,1,Infinity, Infinity, Infinity];
643         radicalAxis = JXG.Math.normalize(radicalAxis);
644         return this.meetLineCircle(radicalAxis, circ1, i, board);
645         // Returns do not work with homogeneous coordinates, yet
646     },
647 
648     /**
649      * Compute an intersection of the curves c1 and c2
650      * with a generalized Newton method.
651      * We want to find values t1, t2 such that
652      * c1(t1) = c2(t2), i.e.
653      * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
654      * We set
655      * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2))
656      *
657      * The Jacobian J is defined by
658      * J = (a, b)
659      *     (c, d)
660      * where
661      * a = c1_x'(t1)
662      * b = -c2_x'(t2)
663      * c = c1_y'(t1)
664      * d = -c2_y'(t2)
665      *
666      * The inverse J^(-1) of J is equal to
667      *  (d, -b)/
668      *  (-c, a) / (ad-bc)
669      *
670      * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
671      * If the function meetCurveCurve possesses the properties
672      * t1memo and t2memo then these are taken as start values
673      * for the Newton algorithm.
674      * After stopping of the Newton algorithm the values of t1 and t2 are stored in
675      * t1memo and t2memo.
676      *
677      * @param {JXG.Curve} c1 Curve, Line or Circle
678      * @param {JXG.Curve} c2 Curve, Line or Circle
679      * @param {Number} t1ini start value for t1
680      * @param {Number} t2ini start value for t2
681      * @param {JXG.Board} [board=c1.board] Reference to a board object.
682      * @returns {JXG.Coords} intersection point
683      */
684     meetCurveCurve: function(c1, c2, t1ini, t2ini, board) {
685         var count = 0,
686             t1, t2,
687             a, b, c, d, disc,
688             e, f, F,
689             D00, D01,
690             D10, D11;
691 
692         if (!JXG.exists(board))
693             board = c1.board;
694 
695         if (arguments.callee.t1memo) {
696             t1 = arguments.callee.t1memo;
697             t2 = arguments.callee.t2memo;
698         } else {
699             t1 = t1ini;
700             t2 = t2ini;
701         }
702         /*
703          if (t1>c1.maxX()) { t1 = c1.maxX(); }
704          if (t1<c1.minX()) { t1 = c1.minX(); }
705          if (t2>c2.maxX()) { t2 = c2.maxX(); }
706          if (t2<c2.minX()) { t2 = c2.minX(); }
707          */
708         e = c1.X(t1) - c2.X(t2);
709         f = c1.Y(t1) - c2.Y(t2);
710         F = e * e + f * f;
711 
712         D00 = c1.board.D(c1.X, c1);
713         D01 = c2.board.D(c2.X, c2);
714         D10 = c1.board.D(c1.Y, c1);
715         D11 = c2.board.D(c2.Y, c2);
716 
717         while (F > JXG.Math.eps && count < 10) {
718             a = D00(t1);
719             b = -D01(t2);
720             c = D10(t1);
721             d = -D11(t2);
722             disc = a * d - b * c;
723             t1 -= (d * e - b * f) / disc;
724             t2 -= (a * f - c * e) / disc;
725             e = c1.X(t1) - c2.X(t2);
726             f = c1.Y(t1) - c2.Y(t2);
727             F = e * e + f * f;
728             count++;
729         }
730         //console.log(t1+' '+t2);
731 
732         arguments.callee.t1memo = t1;
733         arguments.callee.t2memo = t2;
734 
735         //return (new JXG.Coords(JXG.COORDS_BY_USER, [2,2], this.board));
736         if (Math.abs(t1) < Math.abs(t2)) {
737             return (new JXG.Coords(JXG.COORDS_BY_USER, [c1.X(t1),c1.Y(t1)], board));
738         } else {
739             return (new JXG.Coords(JXG.COORDS_BY_USER, [c2.X(t2),c2.Y(t2)], board));
740         }
741     },
742 
743     /**
744      * order of input does not matter for el1 and el2.
745      * @param {JXG.Curve,JXG.Line} el1 Curve or Line
746      * @param {JXG.Curve,JXG.Line} el2 Curve or Line
747      * @param {?} nr
748      * @param {JXG.Board} [board=el1.board] Reference to a board object.
749      * @returns {JXG.Coords} Intersection point
750      */
751     meetCurveLine: function(el1, el2, nr, board) {
752         var t, t2, i, cu, li, func, z,
753             tnew, steps, delta, tstart, tend, cux, cuy;
754 
755         if (!JXG.exists(board))
756             board = el1.board;
757 
758 
759         //for (i=0;i<arguments.length-1;i++) {
760         for (i = 0; i <= 1; i++) {
761             if (arguments[i].elementClass == JXG.OBJECT_CLASS_CURVE) {
762                 cu = arguments[i];
763             } else if (arguments[i].elementClass == JXG.OBJECT_CLASS_LINE) {
764                 li = arguments[i];
765             } else
766                 throw new Error("JSXGraph: Can't call meetCurveLine with parent class " + (arguments[i].elementClass) + ".");
767         }
768 
769         func = function(t) {
770             return li.stdform[0] + li.stdform[1] * cu.X(t) + li.stdform[2] * cu.Y(t);
771         };
772 
773         if (arguments.callee.t1memo) {
774             tstart = arguments.callee.t1memo;
775             t = JXG.Math.Numerics.root(func, tstart);
776         } else {
777             tstart = cu.minX();
778             tend = cu.maxX();
779             t = JXG.Math.Numerics.root(func, [tstart,tend]);
780         }
781         arguments.callee.t1memo = t;
782         cux = cu.X(t);
783         cuy = cu.Y(t);
784 
785         if (nr == 1) {
786             if (arguments.callee.t2memo) {
787                 tstart = arguments.callee.t2memo;
788                 t2 = JXG.Math.Numerics.root(func, tstart);
789             }
790             if (!(Math.abs(t2 - t) > 0.1 && Math.abs(cux - cu.X(t2)) > 0.1 && Math.abs(cuy - cu.Y(t2)) > 0.1)) {
791                 steps = 20;
792                 delta = (cu.maxX() - cu.minX()) / steps;
793                 tnew = cu.minX();
794                 for (i = 0; i < steps; i++) {
795                     t2 = JXG.Math.Numerics.root(func, [tnew,tnew + delta]);
796                     if (Math.abs(t2 - t) > 0.1 && Math.abs(cux - cu.X(t2)) > 0.1 && Math.abs(cuy - cu.Y(t2)) > 0.1) {
797                         break;
798                     }
799                     tnew += delta;
800                 }
801             }
802             t = t2;
803             arguments.callee.t2memo = t;
804         }
805 
806         if (Math.abs(func(t)) > JXG.Math.eps) {
807             z = 0.0;
808         } else {
809             z = 1.0;
810         }
811 
812         return (new JXG.Coords(JXG.COORDS_BY_USER, [z, cu.X(t),cu.Y(t)], board));
813     },
814 
815 
816 
817     /****************************************/
818     /****           PROJECTIONS          ****/
819     /****************************************/
820 
821     /**
822      * Calculates the coordinates of the projection of a given point on a given circle. I.o.w. the
823      * nearest one of the two intersection points of the line through the given point and the circles
824      * midpoint.
825      * @param {JXG.Point} point Point to project.
826      * @param {JXG.Circle} circle Circle on that the point is projected.
827      * @param {JXG.Board} [board=point.board] Reference to the board
828      * @returns {JXG.Coords} The coordinates of the projection of the given point on the given circle.
829      */
830     projectPointToCircle: function(point, circle, board) {
831         var dist = point.coords.distance(JXG.COORDS_BY_USER, circle.midpoint.coords),
832             P = point.coords.usrCoords,
833             M = circle.midpoint.coords.usrCoords,
834             x, y, factor;
835 
836         if (!JXG.exists(board))
837             board = point.board;
838 
839         if (Math.abs(dist) < JXG.Math.eps) {
840             dist = JXG.Math.eps;
841         }
842         factor = circle.Radius() / dist;
843         x = M[1] + factor * (P[1] - M[1]);
844         y = M[2] + factor * (P[2] - M[2]);
845 
846         return new JXG.Coords(JXG.COORDS_BY_USER, [x, y], board);
847     },
848 
849     /**
850      * Calculates the coordinates of the projection of a given point on a given line. I.o.w. the
851      * intersection point of the given line and its perpendicular through the given point.
852      * @param {JXG.Point} point Point to project.
853      * @param {JXG.Line} line Line on that the point is projected.
854      * @param {JXG.Board} [board=point.board] Reference to a board.
855      * @returns {JXG.Coords} The coordinates of the projection of the given point on the given line.
856      */
857     projectPointToLine: function(point, line, board) {
858         // Homogeneous version
859         var v = [0,line.stdform[1],line.stdform[2]];
860 
861         if (!JXG.exists(board))
862             board = point.board;
863 
864         v = JXG.Math.crossProduct(v, point.coords.usrCoords);
865 
866         return this.meetLineLine(v, line.stdform, 0, board);
867     },
868 
869     /**
870      * Calculates the coordinates of the projection of a given point on a given curve.
871      * Uses {@link #projectCoordsToCurve}.
872      * @param {JXG.Point} point Point to project.
873      * @param {JXG.Curve} curve Curve on that the point is projected.
874      * @param {JXG.Board} [board=point.board] Reference to a board.
875      * @see #projectCoordsToCurve
876      * @returns {JXG.Coords} The coordinates of the projection of the given point on the given graph.
877      */
878     projectPointToCurve: function(point, curve, board) {
879         if (!JXG.exists(board))
880             board = point.board;
881 
882         var x = point.X(),
883             y = point.Y(),
884             t = point.position || 0.0, //(curve.minX()+curve.maxX())*0.5,
885             result = this.projectCoordsToCurve(x, y, t, curve, board);
886 
887         point.position = result[1];      // side effect !
888         return result[0];
889     },
890 
891     /**
892      * Calculates the coordinates of the projection of a coordinates pair on a given curve. In case of
893      * function graphs this is the
894      * intersection point of the curve and the parallel to y-axis through the given point.
895      * @param {Number} x coordinate to project.
896      * @param {Number} y coordinate to project.
897      * @param {Number} t start value for newtons method
898      * @param {JXG.Curve} curve Curve on that the point is projected.
899      * @param {JXG.Board} [board=curve.board] Reference to a board.
900      * @see #projectPointToCurve
901      * @returns {JXG.Coords} Array containing the coordinates of the projection of the given point on the given graph and
902      * the position on the curve.
903      */
904     projectCoordsToCurve: function(x, y, t, curve, board) {
905         var newCoords, x0, y0, x1, y1, den, i, mindist, dist, lbda, j,
906             infty = 1000000.0, minfunc, tnew, fnew, fold, delta, steps;
907 
908         if (!JXG.exists(board))
909             board = curve.board;
910 
911         if (curve.curveType == 'parameter' || curve.curveType == 'polar') {
912             // Function to minimize
913             minfunc = function(t) {
914                 var dx = x - curve.X(t),
915                     dy = y - curve.Y(t);
916                 return dx * dx + dy * dy;
917             };
918 
919             fold = minfunc(t);
920             steps = 20;
921             delta = (curve.maxX() - curve.minX()) / steps;
922             tnew = curve.minX();
923             for (j = 0; j < steps; j++) {
924                 fnew = minfunc(tnew);
925                 if (fnew < fold) {
926                     t = tnew;
927                     fold = fnew;
928                 }
929                 tnew += delta;
930             }
931             t = JXG.Math.Numerics.root(JXG.Math.Numerics.D(minfunc), t);
932 
933             if (t < curve.minX()) {
934                 t = curve.maxX() + t - curve.minX();
935             } // Cyclically
936             if (t > curve.maxX()) {
937                 t = curve.minX() + t - curve.maxX();
938             }
939             newCoords = new JXG.Coords(JXG.COORDS_BY_USER, [curve.X(t),curve.Y(t)], board);
940         } else if (curve.curveType == 'plot') {
941             mindist = infty;
942             for (i = 0; i < curve.numberPoints; i++) {
943                 x0 = x - curve.X(i);
944                 y0 = y - curve.Y(i);
945                 dist = Math.sqrt(x0 * x0 + y0 * y0);
946                 if (dist < mindist) {
947                     mindist = dist;
948                     t = i;
949                 }
950                 if (i == curve.numberPoints - 1) {
951                     continue;
952                 }
953 
954                 x1 = curve.X(i + 1) - curve.X(i);
955                 y1 = curve.Y(i + 1) - curve.Y(i);
956                 den = x1 * x1 + y1 * y1;
957                 if (den >= JXG.Math.eps) {
958                     lbda = (x0 * x1 + y0 * y1) / den;
959                     dist = Math.sqrt(x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1));
960                 } else {
961                     lbda = 0.0;
962                     dist = Math.sqrt(x0 * x0 + y0 * y0);
963                 }
964                 if (lbda >= 0.0 && lbda <= 1.0 && dist < mindist) {
965                     t = i + lbda;
966                     mindist = dist;
967                 }
968             }
969             i = Math.floor(t);
970             lbda = t - i;
971             if (i < curve.numberPoints - 1) {
972                 x = lbda * curve.X(i + 1) + (1.0 - lbda) * curve.X(i);
973                 y = lbda * curve.Y(i + 1) + (1.0 - lbda) * curve.Y(i);
974             } else {
975                 x = curve.X(i);
976                 y = curve.Y(i);
977             }
978             newCoords = new JXG.Coords(JXG.COORDS_BY_USER, [x,y], board);
979         } else {             // functiongraph
980             t = x;
981             x = t; //curve.X(t);
982             y = curve.Y(t);
983             newCoords = new JXG.Coords(JXG.COORDS_BY_USER, [x,y], board);
984         }
985         return [curve.updateTransform(newCoords),t];
986     },
987 
988     /**
989      * Calculates the coordinates of the projection of a given point on a given turtle. A turtle consists of
990      * one or more curves of curveType 'plot'. Uses {@link #projectPointToCurve}.
991      * @param {JXG.Point} point Point to project.
992      * @param {JXG.Turtle} turtle on that the point is projected.
993      * @param {JXG.Board} [board=point.board] Reference to a board.
994      * @returns {JXG.Coords} The coordinates of the projection of the given point on the given turtle.
995      */
996     projectPointToTurtle: function(point, turtle, board) {
997         var newCoords, t, x, y, i,
998             np = 0,
999             npmin = 0,
1000             mindist = 1000000.0,
1001             dist, el, minEl,
1002             len = turtle.objects.length;
1003 
1004         if (!JXG.exists(board))
1005             board = point.board;
1006 
1007         for (i = 0; i < len; i++) {  // run through all curves of this turtle
1008             el = turtle.objects[i];
1009             if (el.elementClass == JXG.OBJECT_CLASS_CURVE) {
1010                 newCoords = this.projectPointToCurve(point, el);
1011                 dist = this.distance(newCoords.usrCoords, point.coords.usrCoords);
1012                 if (dist < mindist) {
1013                     x = newCoords.usrCoords[1];
1014                     y = newCoords.usrCoords[2];
1015                     t = point.position;
1016                     mindist = dist;
1017                     minEl = el;
1018                     npmin = np;
1019                 }
1020                 np += el.numberPoints;
1021             }
1022         }
1023         newCoords = new JXG.Coords(JXG.COORDS_BY_USER, [x,y], board);
1024         point.position = t + npmin;
1025         return minEl.updateTransform(newCoords);
1026     }
1027 };
1028