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