1 /* 2 Copyright 2008-2025 3 Matthias Ehmann, 4 Michael Gerhaeuser, 5 Carsten Miller, 6 Alfred Wassermann 7 8 This file is part of JSXGraph. 9 10 JSXGraph is free software dual licensed under the GNU LGPL or MIT License. 11 12 You can redistribute it and/or modify it under the terms of the 13 14 * GNU Lesser General Public License as published by 15 the Free Software Foundation, either version 3 of the License, or 16 (at your option) any later version 17 OR 18 * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT 19 20 JSXGraph is distributed in the hope that it will be useful, 21 but WITHOUT ANY WARRANTY; without even the implied warranty of 22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 23 GNU Lesser General Public License for more details. 24 25 You should have received a copy of the GNU Lesser General Public License and 26 the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/> 27 and <https://opensource.org/licenses/MIT/>. 28 */ 29 30 /*global JXG: true, define: true*/ 31 /*jslint nomen: true, plusplus: true*/ 32 33 import JXG from "../jxg.js"; 34 import Const from "../base/constants.js"; 35 import Coords from "../base/coords.js"; 36 import Mat from "./math.js"; 37 import Extrapolate from "./extrapolate.js"; 38 import Numerics from "./numerics.js"; 39 import Statistics from "./statistics.js"; 40 import Geometry from "./geometry.js"; 41 import IntervalArithmetic from "./ia.js"; 42 import Type from "../utils/type.js"; 43 44 /** 45 * Functions for plotting of curves. 46 * @name JXG.Math.Plot 47 * @exports Mat.Plot as JXG.Math.Plot 48 * @namespace 49 */ 50 Mat.Plot = { 51 /** 52 * Check if at least one point on the curve is finite and real. 53 **/ 54 checkReal: function (points) { 55 var b = false, 56 i, 57 p, 58 len = points.length; 59 60 for (i = 0; i < len; i++) { 61 p = points[i].usrCoords; 62 if (!isNaN(p[1]) && !isNaN(p[2]) && Math.abs(p[0]) > Mat.eps) { 63 b = true; 64 break; 65 } 66 } 67 return b; 68 }, 69 70 //---------------------------------------------------------------------- 71 // Plot algorithm v0 72 //---------------------------------------------------------------------- 73 /** 74 * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#doadvancedplot} is <tt>false</tt>. 75 * @param {JXG.Curve} curve JSXGraph curve element 76 * @param {Number} mi Left bound of curve 77 * @param {Number} ma Right bound of curve 78 * @param {Number} len Number of data points 79 * @returns {JXG.Curve} Reference to the curve object. 80 */ 81 updateParametricCurveNaive: function (curve, mi, ma, len) { 82 var i, 83 t, 84 suspendUpdate = false, 85 stepSize = (ma - mi) / len; 86 87 for (i = 0; i < len; i++) { 88 t = mi + i * stepSize; 89 // The last parameter prevents rounding in usr2screen(). 90 curve.points[i].setCoordinates( 91 Const.COORDS_BY_USER, 92 [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], 93 false 94 ); 95 curve.points[i]._t = t; 96 suspendUpdate = true; 97 } 98 return curve; 99 }, 100 101 //---------------------------------------------------------------------- 102 // Plot algorithm v1 103 //---------------------------------------------------------------------- 104 /** 105 * Crude and cheap test if the segment defined by the two points <tt>(x0, y0)</tt> and <tt>(x1, y1)</tt> is 106 * outside the viewport of the board. All parameters have to be given in screen coordinates. 107 * 108 * @private 109 * @deprecated 110 * @param {Number} x0 111 * @param {Number} y0 112 * @param {Number} x1 113 * @param {Number} y1 114 * @param {JXG.Board} board 115 * @returns {Boolean} <tt>true</tt> if the given segment is outside the visible area. 116 */ 117 isSegmentOutside: function (x0, y0, x1, y1, board) { 118 return ( 119 (y0 < 0 && y1 < 0) || 120 (y0 > board.canvasHeight && y1 > board.canvasHeight) || 121 (x0 < 0 && x1 < 0) || 122 (x0 > board.canvasWidth && x1 > board.canvasWidth) 123 ); 124 }, 125 126 /** 127 * Compares the absolute value of <tt>dx</tt> with <tt>MAXX</tt> and the absolute value of <tt>dy</tt> 128 * with <tt>MAXY</tt>. 129 * 130 * @private 131 * @deprecated 132 * @param {Number} dx 133 * @param {Number} dy 134 * @param {Number} MAXX 135 * @param {Number} MAXY 136 * @returns {Boolean} <tt>true</tt>, if <tt>|dx| < MAXX</tt> and <tt>|dy| < MAXY</tt>. 137 */ 138 isDistOK: function (dx, dy, MAXX, MAXY) { 139 return Math.abs(dx) < MAXX && Math.abs(dy) < MAXY && !isNaN(dx + dy); 140 }, 141 142 /** 143 * @private 144 * @deprecated 145 */ 146 isSegmentDefined: function (x0, y0, x1, y1) { 147 return !(isNaN(x0 + y0) && isNaN(x1 + y1)); 148 }, 149 150 /** 151 * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#doadvancedplot} is <tt>true</tt>. 152 * Since 0.99 this algorithm is deprecated. It still can be used if {@link JXG.Curve#doadvancedplotold} is <tt>true</tt>. 153 * 154 * @deprecated 155 * @param {JXG.Curve} curve JSXGraph curve element 156 * @param {Number} mi Left bound of curve 157 * @param {Number} ma Right bound of curve 158 * @returns {JXG.Curve} Reference to the curve object. 159 */ 160 updateParametricCurveOld: function (curve, mi, ma) { 161 var i, t, d, x, y, 162 x0, y0,// t0, 163 top, 164 depth, 165 MAX_DEPTH, 166 MAX_XDIST, 167 MAX_YDIST, 168 suspendUpdate = false, 169 po = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 170 dyadicStack = [], 171 depthStack = [], 172 pointStack = [], 173 divisors = [], 174 distOK = false, 175 j = 0, 176 distFromLine = function (p1, p2, p0) { 177 var lbda, 178 x0 = p0[1] - p1[1], 179 y0 = p0[2] - p1[2], 180 x1 = p2[0] - p1[1], 181 y1 = p2[1] - p1[2], 182 den = x1 * x1 + y1 * y1; 183 184 if (den >= Mat.eps) { 185 lbda = (x0 * x1 + y0 * y1) / den; 186 if (lbda > 0) { 187 if (lbda <= 1) { 188 x0 -= lbda * x1; 189 y0 -= lbda * y1; 190 // lbda = 1.0; 191 } else { 192 x0 -= x1; 193 y0 -= y1; 194 } 195 } 196 } 197 return Mat.hypot(x0, y0); 198 }; 199 200 JXG.deprecated("Curve.updateParametricCurveOld()"); 201 202 if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) { 203 MAX_DEPTH = 15; 204 MAX_XDIST = 10; // 10 205 MAX_YDIST = 10; // 10 206 } else { 207 MAX_DEPTH = 21; 208 MAX_XDIST = 0.7; // 0.7 209 MAX_YDIST = 0.7; // 0.7 210 } 211 212 divisors[0] = ma - mi; 213 for (i = 1; i < MAX_DEPTH; i++) { 214 divisors[i] = divisors[i - 1] * 0.5; 215 } 216 217 i = 1; 218 dyadicStack[0] = 1; 219 depthStack[0] = 0; 220 221 t = mi; 222 po.setCoordinates( 223 Const.COORDS_BY_USER, 224 [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], 225 false 226 ); 227 228 // Now, there was a first call to the functions defining the curve. 229 // Defining elements like sliders have been evaluated. 230 // Therefore, we can set suspendUpdate to false, so that these defining elements 231 // need not be evaluated anymore for the rest of the plotting. 232 suspendUpdate = true; 233 x0 = po.scrCoords[1]; 234 y0 = po.scrCoords[2]; 235 // t0 = t; 236 237 t = ma; 238 po.setCoordinates( 239 Const.COORDS_BY_USER, 240 [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], 241 false 242 ); 243 x = po.scrCoords[1]; 244 y = po.scrCoords[2]; 245 246 pointStack[0] = [x, y]; 247 248 top = 1; 249 depth = 0; 250 251 curve.points = []; 252 curve.points[j++] = new Coords(Const.COORDS_BY_SCREEN, [x0, y0], curve.board, false); 253 254 do { 255 distOK = 256 this.isDistOK(x - x0, y - y0, MAX_XDIST, MAX_YDIST) || 257 this.isSegmentOutside(x0, y0, x, y, curve.board); 258 while ( 259 depth < MAX_DEPTH && 260 (!distOK || depth < 6) && 261 (depth <= 7 || this.isSegmentDefined(x0, y0, x, y)) 262 ) { 263 // We jump out of the loop if 264 // * depth>=MAX_DEPTH or 265 // * (depth>=6 and distOK) or 266 // * (depth>7 and segment is not defined) 267 268 dyadicStack[top] = i; 269 depthStack[top] = depth; 270 pointStack[top] = [x, y]; 271 top += 1; 272 273 i = 2 * i - 1; 274 // Here, depth is increased and may reach MAX_DEPTH 275 depth++; 276 // In that case, t is undefined and we will see a jump in the curve. 277 t = mi + i * divisors[depth]; 278 279 po.setCoordinates( 280 Const.COORDS_BY_USER, 281 [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], 282 false, 283 true 284 ); 285 x = po.scrCoords[1]; 286 y = po.scrCoords[2]; 287 distOK = 288 this.isDistOK(x - x0, y - y0, MAX_XDIST, MAX_YDIST) || 289 this.isSegmentOutside(x0, y0, x, y, curve.board); 290 } 291 292 if (j > 1) { 293 d = distFromLine( 294 curve.points[j - 2].scrCoords, 295 [x, y], 296 curve.points[j - 1].scrCoords 297 ); 298 if (d < 0.015) { 299 j -= 1; 300 } 301 } 302 303 curve.points[j] = new Coords(Const.COORDS_BY_SCREEN, [x, y], curve.board, false); 304 curve.points[j]._t = t; 305 j += 1; 306 307 x0 = x; 308 y0 = y; 309 // t0 = t; 310 311 top -= 1; 312 x = pointStack[top][0]; 313 y = pointStack[top][1]; 314 depth = depthStack[top] + 1; 315 i = dyadicStack[top] * 2; 316 } while (top > 0 && j < 500000); 317 318 curve.numberPoints = curve.points.length; 319 320 return curve; 321 }, 322 323 //---------------------------------------------------------------------- 324 // Plot algorithm v2 325 //---------------------------------------------------------------------- 326 327 /** 328 * Add a point to the curve plot. If the new point is too close to the previously inserted point, 329 * it is skipped. 330 * Used in {@link JXG.Curve._plotRecursive}. 331 * 332 * @private 333 * @param {JXG.Coords} pnt Coords to add to the list of points 334 */ 335 _insertPoint_v2: function (curve, pnt, t) { 336 var lastReal = !isNaN(this._lastCrds[1] + this._lastCrds[2]), // The last point was real 337 newReal = !isNaN(pnt.scrCoords[1] + pnt.scrCoords[2]), // New point is real point 338 cw = curve.board.canvasWidth, 339 ch = curve.board.canvasHeight, 340 off = 500; 341 342 newReal = 343 newReal && 344 pnt.scrCoords[1] > -off && 345 pnt.scrCoords[2] > -off && 346 pnt.scrCoords[1] < cw + off && 347 pnt.scrCoords[2] < ch + off; 348 349 /* 350 * Prevents two consecutive NaNs or points wich are too close 351 */ 352 if ( 353 (!newReal && lastReal) || 354 (newReal && 355 (!lastReal || 356 Math.abs(pnt.scrCoords[1] - this._lastCrds[1]) > 0.7 || 357 Math.abs(pnt.scrCoords[2] - this._lastCrds[2]) > 0.7)) 358 ) { 359 pnt._t = t; 360 curve.points.push(pnt); 361 this._lastCrds = pnt.copy("scrCoords"); 362 } 363 }, 364 365 /** 366 * Check if there is a single NaN function value at t0. 367 * @param {*} curve 368 * @param {*} t0 369 * @returns {Boolean} true if there is a second NaN point close by, false otherwise 370 */ 371 neighborhood_isNaN_v2: function (curve, t0) { 372 var is_undef, 373 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 374 t, 375 p; 376 377 t = t0 + Mat.eps; 378 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false); 379 p = pnt.usrCoords; 380 is_undef = isNaN(p[1] + p[2]); 381 if (!is_undef) { 382 t = t0 - Mat.eps; 383 pnt.setCoordinates( 384 Const.COORDS_BY_USER, 385 [curve.X(t, true), curve.Y(t, true)], 386 false 387 ); 388 p = pnt.usrCoords; 389 is_undef = isNaN(p[1] + p[2]); 390 if (!is_undef) { 391 return false; 392 } 393 } 394 return true; 395 }, 396 397 /** 398 * Investigate a function term at the bounds of intervals where 399 * the function is not defined, e.g. log(x) at x = 0. 400 * 401 * c is between a and b 402 * @private 403 * @param {JXG.Curve} curve JSXGraph curve element 404 * @param {Array} a Screen coordinates of the left interval bound 405 * @param {Array} b Screen coordinates of the right interval bound 406 * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2 407 * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates 408 * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates 409 * @param {Number} tc (ta + tb) / 2 = tc. Parameter which evaluates to b, i.e. [1, X(tc), Y(tc)] = c in screen coordinates 410 * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0. 411 * @returns {JXG.Boolean} true if the point is inserted and the recursion should stop, false otherwise. 412 */ 413 _borderCase: function (curve, a, b, c, ta, tb, tc, depth) { 414 var t, pnt, p, 415 p_good = null, 416 j, 417 max_it = 30, 418 is_undef = false, 419 t_nan, t_real;// t_real2; 420 // dx, dy, 421 // vx, vy, vx2, vy2; 422 // asymptote; 423 424 if (depth <= 1) { 425 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false); 426 // Test if there is a single undefined point. 427 // If yes, we ignore it. 428 if ( 429 isNaN(a[1] + a[2]) && 430 !isNaN(c[1] + c[2]) && 431 !this.neighborhood_isNaN_v2(curve, ta) 432 ) { 433 return false; 434 } 435 if ( 436 isNaN(b[1] + b[2]) && 437 !isNaN(c[1] + c[2]) && 438 !this.neighborhood_isNaN_v2(curve, tb) 439 ) { 440 return false; 441 } 442 if ( 443 isNaN(c[1] + c[2]) && 444 (!isNaN(a[1] + a[2]) || !isNaN(b[1] + b[2])) && 445 !this.neighborhood_isNaN_v2(curve, tc) 446 ) { 447 return false; 448 } 449 450 j = 0; 451 // Bisect a, b and c until the point t_real is inside of the definition interval 452 // and as close as possible at the boundary. 453 // t_real2 is the second closest point. 454 do { 455 // There are four cases: 456 // a | c | b 457 // --------------- 458 // inf | R | R 459 // R | R | inf 460 // inf | inf | R 461 // R | inf | inf 462 // 463 if (isNaN(a[1] + a[2]) && !isNaN(c[1] + c[2])) { 464 t_nan = ta; 465 t_real = tc; 466 // t_real2 = tb; 467 } else if (isNaN(b[1] + b[2]) && !isNaN(c[1] + c[2])) { 468 t_nan = tb; 469 t_real = tc; 470 // t_real2 = ta; 471 } else if (isNaN(c[1] + c[2]) && !isNaN(b[1] + b[2])) { 472 t_nan = tc; 473 t_real = tb; 474 // t_real2 = tb + (tb - tc); 475 } else if (isNaN(c[1] + c[2]) && !isNaN(a[1] + a[2])) { 476 t_nan = tc; 477 t_real = ta; 478 // t_real2 = ta - (tc - ta); 479 } else { 480 return false; 481 } 482 t = 0.5 * (t_nan + t_real); 483 pnt.setCoordinates( 484 Const.COORDS_BY_USER, 485 [curve.X(t, true), curve.Y(t, true)], 486 false 487 ); 488 p = pnt.usrCoords; 489 490 is_undef = isNaN(p[1] + p[2]); 491 if (is_undef) { 492 t_nan = t; 493 } else { 494 // t_real2 = t_real; 495 t_real = t; 496 } 497 ++j; 498 } while (is_undef && j < max_it); 499 500 // If bisection was successful, take this point. 501 // Useful only for general curves, for function graph 502 // the code below overwrite p_good from here. 503 if (j < max_it) { 504 p_good = p.slice(); 505 c = p.slice(); 506 t_real = t; 507 } 508 509 // OK, bisection has been done now. 510 // t_real contains the closest inner point to the border of the interval we could find. 511 // t_real2 is the second nearest point to this boundary. 512 // Now we approximate the derivative by computing the slope of the line through these two points 513 // and test if it is "infinite", i.e larger than 400 in absolute values. 514 // 515 // vx = curve.X(t_real, true); 516 // vx2 = curve.X(t_real2, true); 517 // vy = curve.Y(t_real, true); 518 // vy2 = curve.Y(t_real2, true); 519 // dx = (vx - vx2) / (t_real - t_real2); 520 // dy = (vy - vy2) / (t_real - t_real2); 521 522 if (p_good !== null) { 523 this._insertPoint_v2( 524 curve, 525 new Coords(Const.COORDS_BY_USER, p_good, curve.board, false) 526 ); 527 return true; 528 } 529 } 530 return false; 531 }, 532 533 /** 534 * Recursive interval bisection algorithm for curve plotting. 535 * Used in {@link JXG.Curve.updateParametricCurve}. 536 * @private 537 * @deprecated 538 * @param {JXG.Curve} curve JSXGraph curve element 539 * @param {Array} a Screen coordinates of the left interval bound 540 * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates 541 * @param {Array} b Screen coordinates of the right interval bound 542 * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates 543 * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0. 544 * @param {Number} delta If the distance of the bisection point at (ta + tb) / 2 from the point (a + b) / 2 is less then delta, 545 * the segment [a,b] is regarded as straight line. 546 * @returns {JXG.Curve} Reference to the curve object. 547 */ 548 _plotRecursive_v2: function (curve, a, ta, b, tb, depth, delta) { 549 var tc, 550 c, 551 ds, 552 mindepth = 0, 553 isSmooth, 554 isJump, 555 isCusp, 556 cusp_threshold = 0.5, 557 jump_threshold = 0.99, 558 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false); 559 560 if (curve.numberPoints > 65536) { 561 return; 562 } 563 564 // Test if the function is undefined in an interval 565 if (depth < this.nanLevel && this._isUndefined(curve, a, ta, b, tb)) { 566 return this; 567 } 568 569 if (depth < this.nanLevel && this._isOutside(a, ta, b, tb, curve.board)) { 570 return this; 571 } 572 573 tc = (ta + tb) * 0.5; 574 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(tc, true), curve.Y(tc, true)], false); 575 c = pnt.scrCoords; 576 577 if (this._borderCase(curve, a, b, c, ta, tb, tc, depth)) { 578 return this; 579 } 580 581 ds = this._triangleDists(a, b, c); // returns [d_ab, d_ac, d_cb, d_cd] 582 583 isSmooth = depth < this.smoothLevel && ds[3] < delta; 584 585 isJump = 586 ( 587 depth <= this.jumpLevel && (isNaN(ds[0]) || isNaN(ds[1]) || isNaN(ds[2])) 588 ) || ( 589 depth < this.jumpLevel && 590 ( 591 ds[2] > jump_threshold * ds[0] || 592 ds[1] > jump_threshold * ds[0] || 593 ds[0] === Infinity || 594 ds[1] === Infinity || 595 ds[2] === Infinity 596 ) 597 ); 598 599 isCusp = depth < this.smoothLevel + 2 && ds[0] < cusp_threshold * (ds[1] + ds[2]); 600 601 if (isCusp) { 602 mindepth = 0; 603 isSmooth = false; 604 } 605 606 --depth; 607 608 if (isJump) { 609 this._insertPoint_v2( 610 curve, 611 new Coords(Const.COORDS_BY_SCREEN, [NaN, NaN], curve.board, false), 612 tc 613 ); 614 } else if (depth <= mindepth || isSmooth) { 615 this._insertPoint_v2(curve, pnt, tc); 616 //if (this._borderCase(a, b, c, ta, tb, tc, depth)) {} 617 } else { 618 this._plotRecursive_v2(curve, a, ta, c, tc, depth, delta); 619 620 if (!isNaN(pnt.scrCoords[1] + pnt.scrCoords[2])) { 621 this._insertPoint_v2(curve, pnt, tc); 622 } 623 624 this._plotRecursive_v2(curve, c, tc, b, tb, depth, delta); 625 } 626 627 return this; 628 }, 629 630 /** 631 * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#plotVersion} is <tt>3</tt>. 632 * 633 * @param {JXG.Curve} curve JSXGraph curve element 634 * @param {Number} mi Left bound of curve 635 * @param {Number} ma Right bound of curve 636 * @returns {JXG.Curve} Reference to the curve object. 637 */ 638 updateParametricCurve_v2: function (curve, mi, ma) { 639 var ta, tb, 640 a, b, 641 suspendUpdate = false, 642 pa = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 643 pb = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 644 depth, 645 delta, 646 w2, 647 // h2, 648 bbox, ret_arr; 649 650 //console.time("plot"); 651 if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) { 652 depth = curve.evalVisProp('recursiondepthlow') || 13; 653 delta = 2; 654 // this.smoothLevel = 5; //depth - 7; 655 this.smoothLevel = depth - 6; 656 this.jumpLevel = 3; 657 } else { 658 depth = curve.evalVisProp('recursiondepthhigh') || 17; 659 delta = 2; 660 // smoothLevel has to be small for graphs in a huge interval. 661 // this.smoothLevel = 3; //depth - 7; // 9 662 this.smoothLevel = depth - 9; // 9 663 this.jumpLevel = 2; 664 } 665 this.nanLevel = depth - 4; 666 667 curve.points = []; 668 669 if (this.xterm === "x") { 670 // For function graphs we can restrict the plot interval 671 // to the visible area + plus margin 672 bbox = curve.board.getBoundingBox(); 673 w2 = (bbox[2] - bbox[0]) * 0.3; 674 // h2 = (bbox[1] - bbox[3]) * 0.3; 675 ta = Math.max(mi, bbox[0] - w2); 676 tb = Math.min(ma, bbox[2] + w2); 677 } else { 678 ta = mi; 679 tb = ma; 680 } 681 pa.setCoordinates( 682 Const.COORDS_BY_USER, 683 [curve.X(ta, suspendUpdate), curve.Y(ta, suspendUpdate)], 684 false 685 ); 686 687 // The first function calls of X() and Y() are done. We can now 688 // switch `suspendUpdate` on. If supported by the functions, this 689 // avoids for the rest of the plotting algorithm, evaluation of any 690 // parent elements. 691 suspendUpdate = true; 692 693 pb.setCoordinates( 694 Const.COORDS_BY_USER, 695 [curve.X(tb, suspendUpdate), curve.Y(tb, suspendUpdate)], 696 false 697 ); 698 699 // Find start and end points of the visible area (plus a certain margin) 700 ret_arr = this._findStartPoint(curve, pa.scrCoords, ta, pb.scrCoords, tb); 701 pa.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false); 702 ta = ret_arr[1]; 703 ret_arr = this._findStartPoint(curve, pb.scrCoords, tb, pa.scrCoords, ta); 704 pb.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false); 705 tb = ret_arr[1]; 706 707 // Save the visible area. 708 // This can be used in Curve.hasPoint(). 709 this._visibleArea = [ta, tb]; 710 711 // Start recursive plotting algorithm 712 a = pa.copy("scrCoords"); 713 b = pb.copy("scrCoords"); 714 pa._t = ta; 715 curve.points.push(pa); 716 this._lastCrds = pa.copy("scrCoords"); // Used in _insertPoint 717 this._plotRecursive_v2(curve, a, ta, b, tb, depth, delta); 718 pb._t = tb; 719 curve.points.push(pb); 720 721 curve.numberPoints = curve.points.length; 722 //console.timeEnd("plot"); 723 724 return curve; 725 }, 726 727 //---------------------------------------------------------------------- 728 // Plot algorithm v3 729 //---------------------------------------------------------------------- 730 /** 731 * 732 * @param {JXG.Curve} curve JSXGraph curve element 733 * @param {*} pnt 734 * @param {*} t 735 * @param {*} depth 736 * @param {*} limes 737 * @private 738 */ 739 _insertLimesPoint: function (curve, pnt, t, depth, limes) { 740 var p0, p1, p2; 741 742 // Ignore jump point if it follows limes 743 if ( 744 (Math.abs(this._lastUsrCrds[1]) === Infinity && 745 Math.abs(limes.left_x) === Infinity) || 746 (Math.abs(this._lastUsrCrds[2]) === Infinity && Math.abs(limes.left_y) === Infinity) 747 ) { 748 // console.log("SKIP:", pnt.usrCoords, this._lastUsrCrds, limes); 749 return; 750 } 751 752 // // Ignore jump left from limes 753 // if (Math.abs(limes.left_x) > 100 * Math.abs(this._lastUsrCrds[1])) { 754 // x = Math.sign(limes.left_x) * Infinity; 755 // } else { 756 // x = limes.left_x; 757 // } 758 // if (Math.abs(limes.left_y) > 100 * Math.abs(this._lastUsrCrds[2])) { 759 // y = Math.sign(limes.left_y) * Infinity; 760 // } else { 761 // y = limes.left_y; 762 // } 763 // //pnt.setCoordinates(Const.COORDS_BY_USER, [x, y], false); 764 765 // Add points at a jump. pnt contains [NaN, NaN] 766 //console.log("Add", t, pnt.usrCoords, limes, depth) 767 p0 = new Coords(Const.COORDS_BY_USER, [limes.left_x, limes.left_y], curve.board); 768 p0._t = t; 769 curve.points.push(p0); 770 771 if ( 772 !isNaN(limes.left_x) && 773 !isNaN(limes.left_y) && 774 !isNaN(limes.right_x) && 775 !isNaN(limes.right_y) && 776 (Math.abs(limes.left_x - limes.right_x) > Mat.eps || 777 Math.abs(limes.left_y - limes.right_y) > Mat.eps) 778 ) { 779 p1 = new Coords(Const.COORDS_BY_SCREEN, pnt, curve.board); 780 p1._t = t; 781 curve.points.push(p1); 782 } 783 784 p2 = new Coords(Const.COORDS_BY_USER, [limes.right_x, limes.right_y], curve.board); 785 p2._t = t; 786 curve.points.push(p2); 787 this._lastScrCrds = p2.copy("scrCoords"); 788 this._lastUsrCrds = p2.copy("usrCoords"); 789 }, 790 791 /** 792 * Add a point to the curve plot. If the new point is too close to the previously inserted point, 793 * it is skipped. 794 * Used in {@link JXG.Curve._plotRecursive}. 795 * 796 * @private 797 * @param {JXG.Curve} curve JSXGraph curve element 798 * @param {JXG.Coords} pnt Coords to add to the list of points 799 */ 800 _insertPoint: function (curve, pnt, t, depth, limes) { 801 var last_is_real = !isNaN(this._lastScrCrds[1] + this._lastScrCrds[2]), // The last point was real 802 point_is_real = !isNaN(pnt[1] + pnt[2]), // New point is real point 803 cw = curve.board.canvasWidth, 804 ch = curve.board.canvasHeight, 805 p, 806 near = 0.8, 807 off = 500; 808 809 if (Type.exists(limes)) { 810 this._insertLimesPoint(curve, pnt, t, depth, limes); 811 return; 812 } 813 814 // Check if point has real coordinates and 815 // coordinates are not too far away from canvas. 816 point_is_real = 817 point_is_real && 818 pnt[1] > -off && 819 pnt[2] > -off && 820 pnt[1] < cw + off && 821 pnt[2] < ch + off; 822 823 // Prevent two consecutive NaNs 824 if (!last_is_real && !point_is_real) { 825 return; 826 } 827 828 // Prevent two consecutive points which are too close 829 if ( 830 point_is_real && 831 last_is_real && 832 Math.abs(pnt[1] - this._lastScrCrds[1]) < near && 833 Math.abs(pnt[2] - this._lastScrCrds[2]) < near 834 ) { 835 return; 836 } 837 838 // Prevent two consecutive points at infinity (either direction) 839 if ( 840 (Math.abs(pnt[1]) === Infinity && Math.abs(this._lastUsrCrds[1]) === Infinity) || 841 (Math.abs(pnt[2]) === Infinity && Math.abs(this._lastUsrCrds[2]) === Infinity) 842 ) { 843 return; 844 } 845 846 //console.log("add", t, pnt.usrCoords, depth) 847 // Add regular point 848 p = new Coords(Const.COORDS_BY_SCREEN, pnt, curve.board); 849 p._t = t; 850 curve.points.push(p); 851 this._lastScrCrds = p.copy("scrCoords"); 852 this._lastUsrCrds = p.copy("usrCoords"); 853 }, 854 855 /** 856 * Compute distances in screen coordinates between the points ab, 857 * ac, cb, and cd, where d = (a + b)/2. 858 * cd is used for the smoothness test, ab, ac, cb are used to detect jumps, cusps and poles. 859 * 860 * @private 861 * @param {Array} a Screen coordinates of the left interval bound 862 * @param {Array} b Screen coordinates of the right interval bound 863 * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2 864 * @returns {Array} array of distances in screen coordinates between: ab, ac, cb, and cd. 865 */ 866 _triangleDists: function (a, b, c) { 867 var d, d_ab, d_ac, d_cb, d_cd; 868 869 d = [a[0] * b[0], (a[1] + b[1]) * 0.5, (a[2] + b[2]) * 0.5]; 870 871 d_ab = Geometry.distance(a, b, 3); 872 d_ac = Geometry.distance(a, c, 3); 873 d_cb = Geometry.distance(c, b, 3); 874 d_cd = Geometry.distance(c, d, 3); 875 876 return [d_ab, d_ac, d_cb, d_cd]; 877 }, 878 879 /** 880 * Test if the function is undefined on an interval: 881 * If the interval borders a and b are undefined, 20 random values 882 * are tested if they are undefined, too. 883 * Only if all values are undefined, we declare the function to be undefined in this interval. 884 * 885 * @private 886 * @param {JXG.Curve} curve JSXGraph curve element 887 * @param {Array} a Screen coordinates of the left interval bound 888 * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates 889 * @param {Array} b Screen coordinates of the right interval bound 890 * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates 891 */ 892 _isUndefined: function (curve, a, ta, b, tb) { 893 var t, i, pnt; 894 895 if (!isNaN(a[1] + a[2]) || !isNaN(b[1] + b[2])) { 896 return false; 897 } 898 899 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false); 900 901 for (i = 0; i < 20; ++i) { 902 t = ta + Math.random() * (tb - ta); 903 pnt.setCoordinates( 904 Const.COORDS_BY_USER, 905 [curve.X(t, true), curve.Y(t, true)], 906 false 907 ); 908 if (!isNaN(pnt.scrCoords[0] + pnt.scrCoords[1] + pnt.scrCoords[2])) { 909 return false; 910 } 911 } 912 913 return true; 914 }, 915 916 /** 917 * Decide if a path segment is too far from the canvas that we do not need to draw it. 918 * @private 919 * @param {Array} a Screen coordinates of the start point of the segment 920 * @param {Array} ta Curve parameter of a (unused). 921 * @param {Array} b Screen coordinates of the end point of the segment 922 * @param {Array} tb Curve parameter of b (unused). 923 * @param {JXG.Board} board 924 * @returns {Boolean} True if the segment is too far away from the canvas, false otherwise. 925 */ 926 _isOutside: function (a, ta, b, tb, board) { 927 var off = 500, 928 cw = board.canvasWidth, 929 ch = board.canvasHeight; 930 931 return !!( 932 (a[1] < -off && b[1] < -off) || 933 (a[2] < -off && b[2] < -off) || 934 (a[1] > cw + off && b[1] > cw + off) || 935 (a[2] > ch + off && b[2] > ch + off) 936 ); 937 }, 938 939 /** 940 * Decide if a point of a curve is too far from the canvas that we do not need to draw it. 941 * @private 942 * @param {Array} a Screen coordinates of the point 943 * @param {JXG.Board} board 944 * @returns {Boolean} True if the point is too far away from the canvas, false otherwise. 945 */ 946 _isOutsidePoint: function (a, board) { 947 var off = 500, 948 cw = board.canvasWidth, 949 ch = board.canvasHeight; 950 951 return !!(a[1] < -off || a[2] < -off || a[1] > cw + off || a[2] > ch + off); 952 }, 953 954 /** 955 * For a curve c(t) defined on the interval [ta, tb] find the first point 956 * which is in the visible area of the board (plus some outside margin). 957 * <p> 958 * This method is necessary to restrict the recursive plotting algorithm 959 * {@link JXG.Curve._plotRecursive} to the visible area and not waste 960 * recursion to areas far outside of the visible area. 961 * <p> 962 * This method can also be used to find the last visible point 963 * by reversing the input parameters. 964 * 965 * @param {JXG.Curve} curve JSXGraph curve element 966 * @param {Array} ta Curve parameter of a. 967 * @param {Array} b Screen coordinates of the end point of the segment (unused) 968 * @param {Array} tb Curve parameter of b 969 * @return {Array} Array of length two containing the screen ccordinates of 970 * the starting point and the curve parameter at this point. 971 * @private 972 */ 973 _findStartPoint: function (curve, a, ta, b, tb) { 974 // The code below is too unstable. 975 // E.g. [function(t) { return Math.pow(t, 2) * (t + 5) * Math.pow(t - 5, 2); }, -8, 8] 976 // Therefore, we return here. 977 return [a, ta]; 978 979 // var i, 980 // delta, 981 // tc, 982 // td, 983 // z, 984 // isFound, 985 // w2, 986 // h2, 987 // pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 988 // steps = 40, 989 // eps = 0.01, 990 // fnX1, 991 // fnX2, 992 // fnY1, 993 // fnY2, 994 // bbox = curve.board.getBoundingBox(); 995 996 // if (true || !this._isOutsidePoint(a, curve.board)) { 997 // return [a, ta]; 998 // } 999 // w2 = (bbox[2] - bbox[0]) * 0.3; 1000 // h2 = (bbox[1] - bbox[3]) * 0.3; 1001 // bbox[0] -= w2; 1002 // bbox[1] += h2; 1003 // bbox[2] += w2; 1004 // bbox[3] -= h2; 1005 1006 // delta = (tb - ta) / steps; 1007 // tc = ta + delta; 1008 // isFound = false; 1009 1010 // fnX1 = function (t) { 1011 // return curve.X(t, true) - bbox[0]; 1012 // }; 1013 // fnY1 = function (t) { 1014 // return curve.Y(t, true) - bbox[1]; 1015 // }; 1016 // fnX2 = function (t) { 1017 // return curve.X(t, true) - bbox[2]; 1018 // }; 1019 // fnY2 = function (t) { 1020 // return curve.Y(t, true) - bbox[3]; 1021 // }; 1022 // for (i = 0; i < steps; ++i) { 1023 // // Left border 1024 // z = bbox[0]; 1025 // td = Numerics.root(fnX1, [tc - delta, tc], curve); 1026 // // td = Numerics.fzero(fnX1, [tc - delta, tc], this); 1027 // // console.log("A", tc - delta, tc, td, Math.abs(this.X(td, true) - z)); 1028 // if (Math.abs(curve.X(td, true) - z) < eps) { 1029 // //} * Math.abs(z)) { 1030 // isFound = true; 1031 // break; 1032 // } 1033 // // Top border 1034 // z = bbox[1]; 1035 // td = Numerics.root(fnY1, [tc - delta, tc], curve); 1036 // // td = Numerics.fzero(fnY1, [tc - delta, tc], this); 1037 // // console.log("B", tc - delta, tc, td, Math.abs(this.Y(td, true) - z)); 1038 // if (Math.abs(curve.Y(td, true) - z) < eps) { 1039 // // * Math.abs(z)) { 1040 // isFound = true; 1041 // break; 1042 // } 1043 // // Right border 1044 // z = bbox[2]; 1045 // td = Numerics.root(fnX2, [tc - delta, tc], curve); 1046 // // td = Numerics.fzero(fnX2, [tc - delta, tc], this); 1047 // // console.log("C", tc - delta, tc, td, Math.abs(this.X(td, true) - z)); 1048 // if (Math.abs(curve.X(td, true) - z) < eps) { 1049 // // * Math.abs(z)) { 1050 // isFound = true; 1051 // break; 1052 // } 1053 // // Bottom border 1054 // z = bbox[3]; 1055 // td = Numerics.root(fnY2, [tc - delta, tc], curve); 1056 // // td = Numerics.fzero(fnY2, [tc - delta, tc], this); 1057 // // console.log("D", tc - delta, tc, td, Math.abs(this.Y(td, true) - z)); 1058 // if (Math.abs(curve.Y(td, true) - z) < eps) { 1059 // // * Math.abs(z)) { 1060 // isFound = true; 1061 // break; 1062 // } 1063 // tc += delta; 1064 // } 1065 // if (isFound) { 1066 // pnt.setCoordinates( 1067 // Const.COORDS_BY_USER, 1068 // [curve.X(td, true), curve.Y(td, true)], 1069 // false 1070 // ); 1071 // return [pnt.scrCoords, td]; 1072 // } 1073 // console.log("TODO _findStartPoint", curve.Y.toString(), tc); 1074 // pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(ta, true), curve.Y(ta, true)], false); 1075 // return [pnt.scrCoords, ta]; 1076 }, 1077 1078 /** 1079 * Investigate a function term at the bounds of intervals where 1080 * the function is not defined, e.g. log(x) at x = 0. 1081 * 1082 * c is inbetween a and b 1083 * 1084 * @param {JXG.Curve} curve JSXGraph curve element 1085 * @param {Array} a Screen coordinates of the left interval bound 1086 * @param {Array} b Screen coordinates of the right interval bound 1087 * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2 1088 * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates 1089 * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates 1090 * @param {Number} tc (ta + tb) / 2 = tc. Parameter which evaluates to b, i.e. [1, X(tc), Y(tc)] = c in screen coordinates 1091 * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0. 1092 * @returns {JXG.Boolean} true if the point is inserted and the recursion should stop, false otherwise. 1093 * 1094 * @private 1095 */ 1096 _getBorderPos: function (curve, ta, a, tc, c, tb, b) { 1097 var t, pnt, p, j, 1098 max_it = 30, 1099 is_undef = false, 1100 t_good, t_bad; 1101 1102 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false); 1103 j = 0; 1104 // Bisect a, b and c until the point t_real is inside of the definition interval 1105 // and as close as possible at the boundary. 1106 // (t_real2 is/was the second closest point). 1107 // There are four cases: 1108 // a | c | b 1109 // --------------- 1110 // inf | R | R 1111 // R | R | inf 1112 // inf | inf | R 1113 // R | inf | inf 1114 // 1115 if (isNaN(a[1] + a[2]) && !isNaN(c[1] + c[2])) { 1116 t_bad = ta; 1117 t_good = tc; 1118 } else if (isNaN(b[1] + b[2]) && !isNaN(c[1] + c[2])) { 1119 t_bad = tb; 1120 t_good = tc; 1121 } else if (isNaN(c[1] + c[2]) && !isNaN(b[1] + b[2])) { 1122 t_bad = tc; 1123 t_good = tb; 1124 } else if (isNaN(c[1] + c[2]) && !isNaN(a[1] + a[2])) { 1125 t_bad = tc; 1126 t_good = ta; 1127 } else { 1128 return false; 1129 } 1130 do { 1131 t = 0.5 * (t_good + t_bad); 1132 pnt.setCoordinates( 1133 Const.COORDS_BY_USER, 1134 [curve.X(t, true), curve.Y(t, true)], 1135 false 1136 ); 1137 p = pnt.usrCoords; 1138 is_undef = isNaN(p[1] + p[2]); 1139 if (is_undef) { 1140 t_bad = t; 1141 } else { 1142 t_good = t; 1143 } 1144 ++j; 1145 } while (j < max_it && Math.abs(t_good - t_bad) > Mat.eps); 1146 return t; 1147 }, 1148 1149 /** 1150 * 1151 * @param {JXG.Curve} curve JSXGraph curve element 1152 * @param {Number} ta 1153 * @param {Number} tb 1154 */ 1155 _getCuspPos: function (curve, ta, tb) { 1156 var a = [curve.X(ta, true), curve.Y(ta, true)], 1157 b = [curve.X(tb, true), curve.Y(tb, true)], 1158 max_func = function (t) { 1159 var c = [curve.X(t, true), curve.Y(t, true)]; 1160 return -( 1161 Mat.hypot(a[0] - c[0], a[1] - c[1]) + 1162 Mat.hypot(b[0] - c[0], b[1] - c[1]) 1163 ); 1164 }; 1165 1166 return Numerics.fminbr(max_func, [ta, tb], curve); 1167 }, 1168 1169 /** 1170 * 1171 * @param {JXG.Curve} curve JSXGraph curve element 1172 * @param {Number} ta 1173 * @param {Number} tb 1174 */ 1175 _getJumpPos: function (curve, ta, tb) { 1176 var max_func = function (t) { 1177 var e = Mat.eps * Mat.eps, 1178 c1 = [curve.X(t, true), curve.Y(t, true)], 1179 c2 = [curve.X(t + e, true), curve.Y(t + e, true)]; 1180 return -Math.abs((c2[1] - c1[1]) / (c2[0] - c1[0])); 1181 }; 1182 1183 return Numerics.fminbr(max_func, [ta, tb], curve); 1184 }, 1185 1186 /** 1187 * 1188 * @param {JXG.Curve} curve JSXGraph curve element 1189 * @param {Number} t 1190 * @private 1191 */ 1192 _getLimits: function (curve, t) { 1193 var res, 1194 step = 2 / (curve.maxX() - curve.minX()), 1195 x_l, 1196 x_r, 1197 y_l, 1198 y_r; 1199 1200 // From left 1201 res = Extrapolate.limit(t, -step, curve.X); 1202 x_l = res[0]; 1203 if (res[1] === "infinite") { 1204 x_l = Math.sign(x_l) * Infinity; 1205 } 1206 1207 res = Extrapolate.limit(t, -step, curve.Y); 1208 y_l = res[0]; 1209 if (res[1] === "infinite") { 1210 y_l = Math.sign(y_l) * Infinity; 1211 } 1212 1213 // From right 1214 res = Extrapolate.limit(t, step, curve.X); 1215 x_r = res[0]; 1216 if (res[1] === "infinite") { 1217 x_r = Math.sign(x_r) * Infinity; 1218 } 1219 1220 res = Extrapolate.limit(t, step, curve.Y); 1221 y_r = res[0]; 1222 if (res[1] === "infinite") { 1223 y_r = Math.sign(y_r) * Infinity; 1224 } 1225 1226 return { 1227 left_x: x_l, 1228 left_y: y_l, 1229 right_x: x_r, 1230 right_y: y_r, 1231 t: t 1232 }; 1233 }, 1234 1235 /** 1236 * 1237 * @param {JXG.Curve} curve JSXGraph curve element 1238 * @param {Array} a 1239 * @param {Number} tc 1240 * @param {Array} c 1241 * @param {Number} tb 1242 * @param {Array} b 1243 * @param {String} may_be_special 1244 * @param {Number} depth 1245 * @private 1246 */ 1247 _getLimes: function (curve, ta, a, tc, c, tb, b, may_be_special, depth) { 1248 var t; 1249 1250 if (may_be_special === "border") { 1251 t = this._getBorderPos(curve, ta, a, tc, c, tb, b); 1252 } else if (may_be_special === "cusp") { 1253 t = this._getCuspPos(curve, ta, tb); 1254 } else if (may_be_special === "jump") { 1255 t = this._getJumpPos(curve, ta, tb); 1256 } 1257 return this._getLimits(curve, t); 1258 }, 1259 1260 /** 1261 * Recursive interval bisection algorithm for curve plotting. 1262 * Used in {@link JXG.Curve.updateParametricCurve}. 1263 * @private 1264 * @param {JXG.Curve} curve JSXGraph curve element 1265 * @param {Array} a Screen coordinates of the left interval bound 1266 * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates 1267 * @param {Array} b Screen coordinates of the right interval bound 1268 * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates 1269 * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0. 1270 * @param {Number} delta If the distance of the bisection point at (ta + tb) / 2 from the point (a + b) / 2 is less then delta, 1271 * the segment [a,b] is regarded as straight line. 1272 * @returns {JXG.Curve} Reference to the curve object. 1273 */ 1274 _plotNonRecursive: function (curve, a, ta, b, tb, d) { 1275 var tc, 1276 c, 1277 ds, 1278 mindepth = 0, 1279 limes = null, 1280 a_nan, 1281 b_nan, 1282 isSmooth = false, 1283 may_be_special = "", 1284 x, 1285 y, 1286 oc, 1287 depth, 1288 ds0, 1289 stack = [], 1290 stack_length = 0, 1291 item; 1292 1293 oc = curve.board.origin.scrCoords; 1294 stack[stack_length++] = [a, ta, b, tb, d, Infinity]; 1295 while (stack_length > 0) { 1296 // item = stack.pop(); 1297 item = stack[--stack_length]; 1298 a = item[0]; 1299 ta = item[1]; 1300 b = item[2]; 1301 tb = item[3]; 1302 depth = item[4]; 1303 ds0 = item[5]; 1304 1305 isSmooth = false; 1306 may_be_special = ""; 1307 limes = null; 1308 //console.log(stack.length, item) 1309 1310 if (curve.points.length > 65536) { 1311 return; 1312 } 1313 1314 if (depth < this.nanLevel) { 1315 // Test if the function is undefined in the whole interval [ta, tb] 1316 if (this._isUndefined(curve, a, ta, b, tb)) { 1317 continue; 1318 } 1319 // Test if the graph is far outside the visible are for the interval [ta, tb] 1320 if (this._isOutside(a, ta, b, tb, curve.board)) { 1321 continue; 1322 } 1323 } 1324 1325 tc = (ta + tb) * 0.5; 1326 1327 // Screen coordinates of point at tc 1328 x = curve.X(tc, true); 1329 y = curve.Y(tc, true); 1330 c = [1, oc[1] + x * curve.board.unitX, oc[2] - y * curve.board.unitY]; 1331 ds = this._triangleDists(a, b, c); // returns [d_ab, d_ac, d_cb, d_cd] 1332 1333 a_nan = isNaN(a[1] + a[2]); 1334 b_nan = isNaN(b[1] + b[2]); 1335 if ((a_nan && !b_nan) || (!a_nan && b_nan)) { 1336 may_be_special = "border"; 1337 } else if ( 1338 ds[0] > 0.66 * ds0 || 1339 ds[0] < this.cusp_threshold * (ds[1] + ds[2]) || 1340 ds[1] > 5 * ds[2] || 1341 ds[2] > 5 * ds[1] 1342 ) { 1343 may_be_special = "cusp"; 1344 } else if ( 1345 ds[2] > this.jump_threshold * ds[0] || 1346 ds[1] > this.jump_threshold * ds[0] || 1347 ds[0] === Infinity || 1348 ds[1] === Infinity || 1349 ds[2] === Infinity 1350 ) { 1351 may_be_special = "jump"; 1352 } 1353 isSmooth = 1354 may_be_special === "" && 1355 depth < this.smoothLevel && 1356 ds[3] < this.smooth_threshold; 1357 1358 if (depth < this.testLevel && !isSmooth) { 1359 if (may_be_special === "") { 1360 isSmooth = true; 1361 } else { 1362 limes = this._getLimes(curve, ta, a, tc, c, tb, b, may_be_special, depth); 1363 } 1364 } 1365 1366 if (limes !== null) { 1367 c = [1, NaN, NaN]; 1368 this._insertPoint(curve, c, tc, depth, limes); 1369 } else if (depth <= mindepth || isSmooth) { 1370 this._insertPoint(curve, c, tc, depth, null); 1371 } else { 1372 stack[stack_length++] = [c, tc, b, tb, depth - 1, ds[0]]; 1373 stack[stack_length++] = [a, ta, c, tc, depth - 1, ds[0]]; 1374 } 1375 } 1376 1377 return this; 1378 }, 1379 1380 /** 1381 * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#plotVersion} is <tt>3</tt>. 1382 * This is an experimental plot version, <b>not recommended</b> to be used. 1383 * @param {JXG.Curve} curve JSXGraph curve element 1384 * @param {Number} mi Left bound of curve 1385 * @param {Number} ma Right bound of curve 1386 * @returns {JXG.Curve} Reference to the curve object. 1387 */ 1388 updateParametricCurve_v3: function (curve, mi, ma) { 1389 var ta, 1390 tb, 1391 a, 1392 b, 1393 suspendUpdate = false, 1394 pa = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 1395 pb = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false), 1396 depth, 1397 w2, // h2, 1398 bbox, 1399 ret_arr; 1400 1401 // console.log("-----------------------------------------------------------"); 1402 // console.time("plot"); 1403 if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) { 1404 depth = curve.evalVisProp('recursiondepthlow') || 14; 1405 } else { 1406 depth = curve.evalVisProp('recursiondepthhigh') || 17; 1407 } 1408 1409 // smoothLevel has to be small for graphs in a huge interval. 1410 this.smoothLevel = 7; //depth - 10; 1411 this.nanLevel = depth - 4; 1412 this.testLevel = 4; 1413 this.cusp_threshold = 0.5; 1414 this.jump_threshold = 0.99; 1415 this.smooth_threshold = 2; 1416 1417 curve.points = []; 1418 1419 if (curve.xterm === "x") { 1420 // For function graphs we can restrict the plot interval 1421 // to the visible area +plus margin 1422 bbox = curve.board.getBoundingBox(); 1423 w2 = (bbox[2] - bbox[0]) * 0.3; 1424 //h2 = (bbox[1] - bbox[3]) * 0.3; 1425 ta = Math.max(mi, bbox[0] - w2); 1426 tb = Math.min(ma, bbox[2] + w2); 1427 } else { 1428 ta = mi; 1429 tb = ma; 1430 } 1431 pa.setCoordinates( 1432 Const.COORDS_BY_USER, 1433 [curve.X(ta, suspendUpdate), curve.Y(ta, suspendUpdate)], 1434 false 1435 ); 1436 1437 // The first function calls of X() and Y() are done. We can now 1438 // switch `suspendUpdate` on. If supported by the functions, this 1439 // avoids for the rest of the plotting algorithm, evaluation of any 1440 // parent elements. 1441 suspendUpdate = true; 1442 1443 pb.setCoordinates( 1444 Const.COORDS_BY_USER, 1445 [curve.X(tb, suspendUpdate), curve.Y(tb, suspendUpdate)], 1446 false 1447 ); 1448 1449 // Find start and end points of the visible area (plus a certain margin) 1450 ret_arr = this._findStartPoint(curve, pa.scrCoords, ta, pb.scrCoords, tb); 1451 pa.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false); 1452 ta = ret_arr[1]; 1453 ret_arr = this._findStartPoint(curve, pb.scrCoords, tb, pa.scrCoords, ta); 1454 pb.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false); 1455 tb = ret_arr[1]; 1456 1457 // Save the visible area. 1458 // This can be used in Curve.hasPoint(). 1459 this._visibleArea = [ta, tb]; 1460 1461 // Start recursive plotting algorithm 1462 a = pa.copy("scrCoords"); 1463 b = pb.copy("scrCoords"); 1464 pa._t = ta; 1465 curve.points.push(pa); 1466 this._lastScrCrds = pa.copy("scrCoords"); // Used in _insertPoint 1467 this._lastUsrCrds = pa.copy("usrCoords"); // Used in _insertPoint 1468 1469 this._plotNonRecursive(curve, a, ta, b, tb, depth); 1470 1471 pb._t = tb; 1472 curve.points.push(pb); 1473 1474 curve.numberPoints = curve.points.length; 1475 // console.timeEnd("plot"); 1476 // console.log("number of points:", this.numberPoints); 1477 1478 return curve; 1479 }, 1480 1481 //---------------------------------------------------------------------- 1482 // Plot algorithm v4 1483 //---------------------------------------------------------------------- 1484 1485 /** 1486 * TODO 1487 * @param {Array} vec 1488 * @param {Number} le 1489 * @param {Number} level 1490 * @returns Object 1491 * @private 1492 */ 1493 _criticalInterval: function (vec, le, level) { 1494 var i, 1495 j, 1496 le1, 1497 med, 1498 sgn, 1499 sgnChange, 1500 isGroup = false, 1501 abs_vec, 1502 last = -Infinity, 1503 very_small = false, 1504 smooth = false, 1505 group = 0, 1506 groups = [], 1507 types = [], 1508 positions = []; 1509 1510 abs_vec = Statistics.abs(vec); 1511 med = Statistics.median(abs_vec); 1512 1513 if (med < 1.0e-7) { 1514 med = 1.0e-7; 1515 very_small = true; 1516 } else { 1517 med *= this.criticalThreshold; 1518 } 1519 1520 //console.log("Median", med); 1521 for (i = 0; i < le; i++) { 1522 // Start a group if not yet done and 1523 // add position to group 1524 if (abs_vec[i] > med /*&& abs_vec[i] > 0.01*/) { 1525 positions.push({ i: i, v: vec[i], group: group }); 1526 last = i; 1527 if (!isGroup) { 1528 isGroup = true; 1529 } 1530 } else { 1531 if (isGroup && i > last + 4) { 1532 // End the group 1533 if (positions.length > 0) { 1534 groups.push(positions.slice(0)); 1535 } 1536 positions = []; 1537 isGroup = false; 1538 group++; 1539 } 1540 } 1541 } 1542 if (isGroup) { 1543 if (positions.length > 1) { 1544 groups.push(positions.slice(0)); 1545 } 1546 } 1547 1548 if (very_small && groups.length === 0) { 1549 smooth = true; 1550 } 1551 1552 // Decide if there is a singular critical point 1553 // or if a whole interval is problematic. 1554 // The latter is the case if the differences have many sign changes. 1555 for (j = 0; j < groups.length; j++) { 1556 types[j] = "point"; 1557 le1 = groups[j].length; 1558 if (le1 < 64) { 1559 continue; 1560 } 1561 sgnChange = 0; 1562 sgn = Math.sign(groups[j][0].v); 1563 for (i = 1; i < le1; i++) { 1564 if (Math.sign(groups[j][i].v) !== sgn) { 1565 sgnChange++; 1566 sgn = Math.sign(groups[j][i].v); 1567 } 1568 } 1569 if (sgnChange * 6 > le1) { 1570 types[j] = "interval"; 1571 } 1572 } 1573 1574 return { smooth: smooth, groups: groups, types: types }; 1575 }, 1576 1577 Component: function () { 1578 this.left_isNaN = false; 1579 this.right_isNaN = false; 1580 this.left_t = null; 1581 this.right_t = null; 1582 this.t_values = []; 1583 this.x_values = []; 1584 this.y_values = []; 1585 this.len = 0; 1586 }, 1587 1588 findComponents: function (curve, mi, ma, steps) { 1589 var i, t, h, 1590 x, y, 1591 components = [], 1592 comp, 1593 comp_nr = 0, 1594 cnt = 0, 1595 cntNaNs = 0, 1596 comp_started = false, 1597 suspended = false; 1598 1599 h = (ma - mi) / steps; 1600 components[comp_nr] = new this.Component(); 1601 comp = components[comp_nr]; 1602 1603 for (i = 0, t = mi; i <= steps; i++, t += h) { 1604 x = curve.X(t, suspended); 1605 y = curve.Y(t, suspended); 1606 1607 if (isNaN(x) || isNaN(y)) { 1608 cntNaNs++; 1609 // Wait for - at least - two consecutive NaNs 1610 // This avoids starting a new component if 1611 // the function value has infinity as intermediate value. 1612 if (cntNaNs > 1 && comp_started) { 1613 // Finalize a component 1614 comp.right_isNaN = true; 1615 comp.right_t = t - h; 1616 comp.len = cnt; 1617 1618 // Prepare a new component 1619 comp_started = false; 1620 comp_nr++; 1621 components[comp_nr] = new this.Component(); 1622 comp = components[comp_nr]; 1623 cntNaNs = 0; 1624 } 1625 } else { 1626 // Now there is a non-NaN entry. 1627 if (!comp_started) { 1628 // Start the component 1629 comp_started = true; 1630 cnt = 0; 1631 if (cntNaNs > 0) { 1632 comp.left_t = t - h; 1633 comp.left_isNaN = true; 1634 } 1635 } 1636 cntNaNs = 0; 1637 // Add the value to the component 1638 comp.t_values[cnt] = t; 1639 comp.x_values[cnt] = x; 1640 comp.y_values[cnt] = y; 1641 cnt++; 1642 } 1643 if (i === 0) { 1644 suspended = true; 1645 } 1646 } 1647 if (comp_started) { 1648 comp.len = cnt; 1649 } else { 1650 components.pop(); 1651 } 1652 1653 return components; 1654 }, 1655 1656 getPointType: function (curve, pos, t_approx, t_values, x_table, y_table, len) { 1657 var x_values = x_table[0], 1658 y_values = y_table[0], 1659 full_len = t_values.length, 1660 result = { 1661 idx: pos, 1662 t: t_approx, //t_values[pos], 1663 x: x_values[pos], 1664 y: y_values[pos], 1665 type: "other" 1666 }; 1667 1668 if (pos < 5) { 1669 result.type = "borderleft"; 1670 result.idx = 0; 1671 result.t = t_values[0]; 1672 result.x = x_values[0]; 1673 result.y = y_values[0]; 1674 1675 // console.log('Border left', result.t); 1676 return result; 1677 } 1678 if (pos > len - 6) { 1679 result.type = "borderright"; 1680 result.idx = full_len - 1; 1681 result.t = t_values[full_len - 1]; 1682 result.x = x_values[full_len - 1]; 1683 result.y = y_values[full_len - 1]; 1684 1685 // console.log('Border right', result.t, full_len - 1); 1686 return result; 1687 } 1688 1689 return result; 1690 }, 1691 1692 newtonApprox: function (idx, t, h, level, table) { 1693 var i, 1694 s = 0.0; 1695 for (i = level; i > 0; i--) { 1696 s = ((s + table[i][idx]) * (t - (i - 1) * h)) / i; 1697 } 1698 return s + table[0][idx]; 1699 }, 1700 1701 // Thiele's interpolation formula, 1702 // https://en.wikipedia.org/wiki/Thiele%27s_interpolation_formula 1703 // unused 1704 thiele: function (t, recip, t_values, idx, degree) { 1705 var i, 1706 v = 0.0; 1707 for (i = degree; i > 1; i--) { 1708 v = (t - t_values[idx + i]) / (recip[i][idx + 1] - recip[i - 2][idx + 1] + v); 1709 } 1710 return recip[0][idx + 1] + (t - t_values[idx + 1]) / (recip[1][idx + 1] + v); 1711 }, 1712 1713 differenceMethodExperiments: function (component, curve) { 1714 var i, 1715 level, 1716 le, 1717 up, 1718 t_values = component.t_values, 1719 x_values = component.x_values, 1720 y_values = component.y_values, 1721 x_diffs = [], 1722 y_diffs = [], 1723 x_slopes = [], 1724 y_slopes = [], 1725 x_table = [], 1726 y_table = [], 1727 x_recip = [], 1728 y_recip = [], 1729 h, 1730 numerator, 1731 // x_med, y_med, 1732 foundCriticalPoint = 0, 1733 pos, 1734 ma, 1735 j, 1736 v, 1737 groups, 1738 criticalPoints = []; 1739 1740 h = t_values[1] - t_values[0]; 1741 x_table.push([]); 1742 y_table.push([]); 1743 x_recip.push([]); 1744 y_recip.push([]); 1745 le = y_values.length; 1746 for (i = 0; i < le; i++) { 1747 x_table[0][i] = x_values[i]; 1748 y_table[0][i] = y_values[i]; 1749 x_recip[0][i] = x_values[i]; 1750 y_recip[0][i] = y_values[i]; 1751 } 1752 1753 x_table.push([]); 1754 y_table.push([]); 1755 x_recip.push([]); 1756 y_recip.push([]); 1757 numerator = h; 1758 le = y_values.length - 1; 1759 for (i = 0; i < le; i++) { 1760 x_diffs[i] = x_values[i + 1] - x_values[i]; 1761 y_diffs[i] = y_values[i + 1] - y_values[i]; 1762 x_slopes[i] = x_diffs[i]; 1763 y_slopes[i] = y_diffs[i]; 1764 x_table[1][i] = x_diffs[i]; 1765 y_table[1][i] = y_diffs[i]; 1766 x_recip[1][i] = numerator / x_diffs[i]; 1767 y_recip[1][i] = numerator / y_diffs[i]; 1768 } 1769 le--; 1770 1771 up = Math.min(8, y_values.length - 1); 1772 for (level = 1; level < up; level++) { 1773 x_table.push([]); 1774 y_table.push([]); 1775 x_recip.push([]); 1776 y_recip.push([]); 1777 numerator *= h; 1778 for (i = 0; i < le; i++) { 1779 x_diffs[i] = x_diffs[i + 1] - x_diffs[i]; 1780 y_diffs[i] = y_diffs[i + 1] - y_diffs[i]; 1781 x_table[level + 1][i] = x_diffs[i]; 1782 y_table[level + 1][i] = y_diffs[i]; 1783 x_recip[level + 1][i] = 1784 numerator / (x_recip[level][i + 1] - x_recip[level][i]) + 1785 x_recip[level - 1][i + 1]; 1786 y_recip[level + 1][i] = 1787 numerator / (y_recip[level][i + 1] - y_recip[level][i]) + 1788 y_recip[level - 1][i + 1]; 1789 } 1790 1791 // if (level == 1) { 1792 // console.log("bends level=", level, y_diffs.toString()); 1793 // } 1794 1795 // Store point location which may be centered around 1796 // critical points. 1797 // If the level is suitable, step out of the loop. 1798 groups = this._criticalPoints(y_diffs, le, level); 1799 if (groups === false) { 1800 // Its seems, the degree of the polynomial is equal to level 1801 console.log("Polynomial of degree", level); 1802 groups = []; 1803 break; 1804 } 1805 if (groups.length > 0) { 1806 foundCriticalPoint++; 1807 if (foundCriticalPoint > 1 && level % 2 === 0) { 1808 break; 1809 } 1810 } 1811 le--; 1812 } 1813 1814 // console.log("Last diffs", y_diffs, "level", level); 1815 1816 // Analyze the groups which have been found. 1817 for (i = 0; i < groups.length; i++) { 1818 // console.log("Group", i, groups[i]) 1819 // Identify the maximum difference, i.e. the center of the "problem" 1820 ma = -Infinity; 1821 for (j = 0; j < groups[i].length; j++) { 1822 v = Math.abs(groups[i][j].v); 1823 if (v > ma) { 1824 ma = v; 1825 pos = j; 1826 } 1827 } 1828 pos = Math.floor(groups[i][pos].i + level / 2); 1829 // Analyze the critical point 1830 criticalPoints.push( 1831 this.getPointType( 1832 curve, 1833 pos, 1834 t_values, 1835 x_values, 1836 y_values, 1837 x_slopes, 1838 y_slopes, 1839 le + 1 1840 ) 1841 ); 1842 } 1843 1844 return [criticalPoints, x_table, y_table, x_recip, y_recip]; 1845 }, 1846 1847 getCenterOfCriticalInterval: function (group, degree, t_values) { 1848 var ma, 1849 j, 1850 pos, 1851 v, 1852 num = 0.0, 1853 den = 0.0, 1854 h = t_values[1] - t_values[0], 1855 pos_mean, 1856 range = []; 1857 1858 // Identify the maximum difference, i.e. the center of the "problem" 1859 // If there are several equal maxima, store the positions 1860 // in the array range and determine the center of the array. 1861 1862 ma = -Infinity; 1863 range = []; 1864 for (j = 0; j < group.length; j++) { 1865 v = Math.abs(group[j].v); 1866 if (v > ma) { 1867 range = [j]; 1868 ma = v; 1869 pos = j; 1870 } else if (ma === v) { 1871 range.push(j); 1872 } 1873 } 1874 if (range.length > 0) { 1875 pos_mean = 1876 range.reduce(function (total, val) { 1877 return total + val; 1878 }, 0) / range.length; 1879 pos = Math.floor(pos_mean); 1880 pos_mean += group[0].i; 1881 } 1882 1883 if (ma < Infinity) { 1884 for (j = 0; j < group.length; j++) { 1885 num += Math.abs(group[j].v) * group[j].i; 1886 den += Math.abs(group[j].v); 1887 } 1888 pos_mean = num / den; 1889 } 1890 pos_mean += degree / 2; 1891 return [ 1892 group[pos].i + degree / 2, 1893 pos_mean, 1894 t_values[Math.floor(pos_mean)] + h * (pos_mean - Math.floor(pos_mean)) 1895 ]; 1896 }, 1897 1898 differenceMethod: function (component, curve) { 1899 var i, 1900 level, 1901 le, 1902 up, 1903 t_values = component.t_values, 1904 x_values = component.x_values, 1905 y_values = component.y_values, 1906 x_table = [], 1907 y_table = [], 1908 foundCriticalPoint = 0, 1909 degree_x = -1, 1910 degree_y = -1, 1911 pos, 1912 res, 1913 res_x, 1914 res_y, 1915 t_approx, 1916 groups = [], 1917 types, 1918 criticalPoints = []; 1919 1920 le = y_values.length; 1921 // x_table.push([]); 1922 // y_table.push([]); 1923 // for (i = 0; i < le; i++) { 1924 // x_table[0][i] = x_values[i]; 1925 // y_table[0][i] = y_values[i]; 1926 // } 1927 x_table.push(new Float64Array(x_values)); 1928 y_table.push(new Float64Array(y_values)); 1929 1930 le--; 1931 up = Math.min(12, le); 1932 for (level = 0; level < up; level++) { 1933 // Old style method: 1934 // x_table.push([]); 1935 // y_table.push([]); 1936 // for (i = 0; i < le; i++) { 1937 // x_table[level + 1][i] = x_table[level][i + 1] - x_table[level][i]; 1938 // y_table[level + 1][i] = y_table[level][i + 1] - y_table[level][i]; 1939 // } 1940 // New method: 1941 x_table.push(new Float64Array(le)); 1942 y_table.push(new Float64Array(le)); 1943 x_table[level + 1] = x_table[level].map(function (v, idx, arr) { 1944 return arr[idx + 1] - v; 1945 }); 1946 y_table[level + 1] = y_table[level].map(function (v, idx, arr) { 1947 return arr[idx + 1] - v; 1948 }); 1949 1950 // Store point location which may be centered around critical points. 1951 // If the level is suitable, step out of the loop. 1952 res_y = this._criticalInterval(y_table[level + 1], le, level); 1953 if (res_y.smooth === true) { 1954 // Its seems, the degree of the polynomial is equal to level 1955 // If the values in level + 1 are zero, it might be a polynomial of degree level. 1956 // Seems to work numerically stable until degree 6. 1957 degree_y = level; 1958 groups = []; 1959 } 1960 res_x = this._criticalInterval(x_table[level + 1], le, level); 1961 if (degree_x === -1 && res_x.smooth === true) { 1962 // Its seems, the degree of the polynomial is equal to level 1963 // If the values in level + 1 are zero, it might be a polynomial of degree level. 1964 // Seems to work numerically stable until degree 6. 1965 degree_x = level; 1966 } 1967 if (degree_y >= 0) { 1968 break; 1969 } 1970 1971 if (res_y.groups.length > 0) { 1972 foundCriticalPoint++; 1973 if (foundCriticalPoint > 2 && (level + 1) % 2 === 0) { 1974 groups = res_y.groups; 1975 types = res_y.types; 1976 break; 1977 } 1978 } 1979 le--; 1980 } 1981 1982 // console.log("Last diffs", y_table[Math.min(level + 1, up)], "level", level + 1); 1983 // Analyze the groups which have been found. 1984 for (i = 0; i < groups.length; i++) { 1985 if (types[i] === "interval") { 1986 continue; 1987 } 1988 // console.log("Group", i, groups[i], types[i], level + 1) 1989 res = this.getCenterOfCriticalInterval(groups[i], level + 1, t_values); 1990 pos = res_y[0]; 1991 pos = Math.floor(res[1]); 1992 t_approx = res[2]; 1993 // console.log("Critical points:", groups, res, pos) 1994 1995 // Analyze the type of the critical point 1996 // Result is of type 'borderleft', borderright', 'other' 1997 criticalPoints.push( 1998 this.getPointType(curve, pos, t_approx, t_values, x_table, y_table, le + 1) 1999 ); 2000 } 2001 2002 // if (level === up) { 2003 // console.log("No convergence!"); 2004 // } else { 2005 // console.log("Convergence level", level); 2006 // } 2007 return [criticalPoints, x_table, y_table, degree_x, degree_y]; 2008 }, 2009 2010 _insertPoint_v4: function (curve, crds, t, doLog) { 2011 var p, 2012 prev = null, 2013 x, 2014 y, 2015 near = 0.8; 2016 2017 if (curve.points.length > 0) { 2018 prev = curve.points[curve.points.length - 1].scrCoords; 2019 } 2020 2021 // Add regular point 2022 p = new Coords(Const.COORDS_BY_USER, crds, curve.board); 2023 2024 if (prev !== null) { 2025 x = p.scrCoords[1] - prev[1]; 2026 y = p.scrCoords[2] - prev[2]; 2027 if (x * x + y * y < near * near) { 2028 // Math.abs(p.scrCoords[1] - prev[1]) < near && 2029 // Math.abs(p.scrCoords[2] - prev[2]) < near) { 2030 return; 2031 } 2032 } 2033 2034 p._t = t; 2035 curve.points.push(p); 2036 }, 2037 2038 getInterval: function (curve, ta, tb) { 2039 var t_int, 2040 // x_int, 2041 y_int; 2042 2043 //console.log('critical point', ta, tb); 2044 IntervalArithmetic.disable(); 2045 2046 t_int = IntervalArithmetic.Interval(ta, tb); 2047 curve.board.mathLib = IntervalArithmetic; 2048 curve.board.mathLibJXG = IntervalArithmetic; 2049 // x_int = curve.X(t_int, true); 2050 y_int = curve.Y(t_int, true); 2051 curve.board.mathLib = Math; 2052 curve.board.mathLibJXG = JXG.Math; 2053 2054 //console.log(x_int, y_int); 2055 return y_int; 2056 }, 2057 2058 sign: function (v) { 2059 if (v < 0) { 2060 return -1; 2061 } 2062 if (v > 0) { 2063 return 1; 2064 } 2065 return 0; 2066 }, 2067 2068 handleBorder: function (curve, comp, group, x_table, y_table) { 2069 var idx = group.idx, 2070 t, 2071 t1, 2072 t2, 2073 size = 32, 2074 y_int, 2075 x, 2076 y, 2077 lo, 2078 hi, 2079 i, 2080 components2, 2081 le, 2082 h; 2083 2084 // console.log("HandleBorder at t =", t_approx); 2085 // console.log("component:", comp) 2086 // console.log("Group:", group); 2087 2088 h = comp.t_values[1] - comp.t_values[0]; 2089 if (group.type === "borderleft") { 2090 t = comp.left_isNaN ? comp.left_t : group.t - h; 2091 t1 = t; 2092 t2 = t1 + h; 2093 } else if (group.type === "borderright") { 2094 t = comp.right_isNaN ? comp.right_t : group.t + h; 2095 t2 = t; 2096 t1 = t2 - h; 2097 } else { 2098 console.log("No bordercase!!!"); 2099 } 2100 2101 components2 = this.findComponents(curve, t1, t2, size); 2102 if (components2.length === 0) { 2103 return; 2104 } 2105 if (group.type === "borderleft") { 2106 t1 = components2[0].left_t; 2107 t2 = components2[0].t_values[0]; 2108 h = components2[0].t_values[1] - components2[0].t_values[0]; 2109 t1 = t1 === null ? t2 - h : t1; 2110 t = t1; 2111 y_int = this.getInterval(curve, t1, t2); 2112 if (Type.isObject(y_int)) { 2113 lo = y_int.lo; 2114 hi = y_int.hi; 2115 2116 x = curve.X(t, true); 2117 y = y_table[1][idx] < 0 ? hi : lo; 2118 this._insertPoint_v4(curve, [1, x, y], t); 2119 } 2120 } 2121 2122 le = components2[0].t_values.length; 2123 for (i = 0; i < le; i++) { 2124 t = components2[0].t_values[i]; 2125 x = components2[0].x_values[i]; 2126 y = components2[0].y_values[i]; 2127 this._insertPoint_v4(curve, [1, x, y], t); 2128 } 2129 2130 if (group.type === "borderright") { 2131 t1 = components2[0].t_values[le - 1]; 2132 t2 = components2[0].right_t; 2133 h = components2[0].t_values[1] - components2[0].t_values[0]; 2134 t2 = t2 === null ? t1 + h : t2; 2135 2136 t = t2; 2137 y_int = this.getInterval(curve, t1, t2); 2138 if (Type.isObject(y_int)) { 2139 lo = y_int.lo; 2140 hi = y_int.hi; 2141 x = curve.X(t, true); 2142 y = y_table[1][idx] > 0 ? hi : lo; 2143 this._insertPoint_v4(curve, [1, x, y], t); 2144 } 2145 } 2146 }, 2147 2148 _seconditeration_v4: function (curve, comp, group, x_table, y_table) { 2149 var i, t1, t2, ret, components2, comp2, idx, groups2, g, x_table2, y_table2, start, le; 2150 2151 // Look at two points, hopefully left and right from the critical point 2152 t1 = comp.t_values[group.idx - 2]; 2153 t2 = comp.t_values[group.idx + 2]; 2154 components2 = this.findComponents(curve, t1, t2, 64); 2155 for (idx = 0; idx < components2.length; idx++) { 2156 comp2 = components2[idx]; 2157 ret = this.differenceMethod(comp2, curve); 2158 groups2 = ret[0]; 2159 x_table2 = ret[1]; 2160 y_table2 = ret[2]; 2161 start = 0; 2162 for (g = 0; g <= groups2.length; g++) { 2163 if (g === groups2.length) { 2164 le = comp2.len; 2165 } else { 2166 le = groups2[g].idx; 2167 } 2168 2169 // Insert all uncritical points until next critical point 2170 for (i = start; i < le; i++) { 2171 if (!isNaN(comp2.x_values[i]) && !isNaN(comp2.y_values[i])) { 2172 this._insertPoint_v4( 2173 curve, 2174 [1, comp2.x_values[i], comp2.y_values[i]], 2175 comp2.t_values[i] 2176 ); 2177 } 2178 } 2179 // Handle next critical point 2180 if (g < groups2.length) { 2181 this.handleSingularity(curve, comp2, groups2[g], x_table2, y_table2); 2182 start = groups2[g].idx + 1; 2183 } 2184 } 2185 le = comp2.len; 2186 if (idx < components2.length - 1) { 2187 this._insertPoint_v4(curve, [1, NaN, NaN], comp2.right_t); 2188 } 2189 } 2190 return this; 2191 }, 2192 2193 _recurse_v4: function (curve, t1, t2, x1, y1, x2, y2, level) { 2194 var tol = 2, 2195 t = (t1 + t2) * 0.5, 2196 x = curve.X(t, true), 2197 y = curve.Y(t, true), 2198 dx, 2199 dy; 2200 2201 //console.log("Level", level) 2202 if (level === 0) { 2203 this._insertPoint_v4(curve, [1, NaN, NaN], t); 2204 return; 2205 } 2206 // console.log("R", t1, t2) 2207 dx = (x - x1) * curve.board.unitX; 2208 dy = (y - y1) * curve.board.unitY; 2209 // console.log("D1", Math.sqrt(dx * dx + dy * dy)) 2210 if (Mat.hypot(dx, dy) > tol) { 2211 this._recurse_v4(curve, t1, t, x1, y1, x, y, level - 1); 2212 } else { 2213 this._insertPoint_v4(curve, [1, x, y], t); 2214 } 2215 dx = (x - x2) * curve.board.unitX; 2216 dy = (y - y2) * curve.board.unitY; 2217 // console.log("D2", Math.sqrt(dx * dx + dy * dy), x-x2, y-y2) 2218 if (Mat.hypot(dx, dy) > tol) { 2219 this._recurse_v4(curve, t, t2, x, y, x2, y2, level - 1); 2220 } else { 2221 this._insertPoint_v4(curve, [1, x, y], t); 2222 } 2223 }, 2224 2225 handleSingularity: function (curve, comp, group, x_table, y_table) { 2226 var idx = group.idx, 2227 t, 2228 t1, 2229 t2, 2230 y_int, 2231 i1, 2232 i2, 2233 x, 2234 // y, 2235 lo, 2236 hi, 2237 d_lft, 2238 d_rgt, 2239 d_thresh = 100, 2240 // d1, 2241 // d2, 2242 di1 = 5, 2243 di2 = 3; 2244 2245 t = group.t; 2246 console.log("HandleSingularity at t =", t); 2247 // console.log(comp.t_values[idx - 1], comp.y_values[idx - 1], comp.t_values[idx + 1], comp.y_values[idx + 1]); 2248 // console.log(group); 2249 2250 // Look at two points, hopefully left and right from the critical point 2251 t1 = comp.t_values[idx - di1]; 2252 t2 = comp.t_values[idx + di1]; 2253 2254 y_int = this.getInterval(curve, t1, t2); 2255 if (Type.isObject(y_int)) { 2256 lo = y_int.lo; 2257 hi = y_int.hi; 2258 } else { 2259 if (y_table[0][idx - 1] < y_table[0][idx + 1]) { 2260 lo = y_table[0][idx - 1]; 2261 hi = y_table[0][idx + 1]; 2262 } else { 2263 lo = y_table[0][idx + 1]; 2264 hi = y_table[0][idx - 1]; 2265 } 2266 } 2267 2268 x = curve.X(t, true); 2269 2270 d_lft = 2271 (y_table[0][idx - di2] - y_table[0][idx - di1]) / 2272 (comp.t_values[idx - di2] - comp.t_values[idx - di1]); 2273 d_rgt = 2274 (y_table[0][idx + di2] - y_table[0][idx + di1]) / 2275 (comp.t_values[idx + di2] - comp.t_values[idx + di1]); 2276 2277 console.log(":::", d_lft, d_rgt); 2278 2279 //this._insertPoint_v4(curve, [1, NaN, NaN], 0); 2280 2281 if (d_lft < -d_thresh) { 2282 // Left branch very steep downwards -> add the minimum 2283 this._insertPoint_v4(curve, [1, x, lo], t, true); 2284 if (d_rgt <= d_thresh) { 2285 // Right branch not very steep upwards -> interrupt the curve 2286 // I.e. it looks like -infty / (finite or infty) and not like -infty / -infty 2287 this._insertPoint_v4(curve, [1, NaN, NaN], t); 2288 } 2289 } else if (d_lft > d_thresh) { 2290 // Left branch very steep upwards -> add the maximum 2291 this._insertPoint_v4(curve, [1, x, hi], t); 2292 if (d_rgt >= -d_thresh) { 2293 // Right branch not very steep downwards -> interrupt the curve 2294 // I.e. it looks like infty / (finite or -infty) and not like infty / infty 2295 this._insertPoint_v4(curve, [1, NaN, NaN], t); 2296 } 2297 } else { 2298 if (lo === -Infinity) { 2299 this._insertPoint_v4(curve, [1, x, lo], t, true); 2300 this._insertPoint_v4(curve, [1, NaN, NaN], t); 2301 } 2302 if (hi === Infinity) { 2303 this._insertPoint_v4(curve, [1, NaN, NaN], t); 2304 this._insertPoint_v4(curve, [1, x, hi], t, true); 2305 } 2306 2307 if (group.t < comp.t_values[idx]) { 2308 i1 = idx - 1; 2309 i2 = idx; 2310 } else { 2311 i1 = idx; 2312 i2 = idx + 1; 2313 } 2314 t1 = comp.t_values[i1]; 2315 t2 = comp.t_values[i2]; 2316 this._recurse_v4( 2317 curve, 2318 t1, 2319 t2, 2320 x_table[0][i1], 2321 y_table[0][i1], 2322 x_table[0][i2], 2323 y_table[0][i2], 2324 10 2325 ); 2326 2327 // x = (x_table[0][idx] - x_table[0][idx - 1]) * curve.board.unitX; 2328 // y = (y_table[0][idx] - y_table[0][idx - 1]) * curve.board.unitY; 2329 // d1 = Math.sqrt(x * x + y * y); 2330 // x = (x_table[0][idx + 1] - x_table[0][idx]) * curve.board.unitX; 2331 // y = (y_table[0][idx + 1] - y_table[0][idx]) * curve.board.unitY; 2332 // d2 = Math.sqrt(x * x + y * y); 2333 2334 // console.log("end", t1, t2, t); 2335 // if (true || (d1 > 2 || d2 > 2)) { 2336 2337 // console.log(d1, d2, y_table[0][idx]) 2338 // // Finite jump 2339 // this._insertPoint_v4(curve, [1, NaN, NaN], t); 2340 // } else { 2341 // if (lo !== -Infinity && hi !== Infinity) { 2342 // // Critical point which can be ignored 2343 // this._insertPoint_v4(curve, [1, x_table[0][idx], y_table[0][idx]], comp.t_values[idx]); 2344 // } else { 2345 // if (lo === -Infinity) { 2346 // this._insertPoint_v4(curve, [1, x, lo], t, true); 2347 // this._insertPoint_v4(curve, [1, NaN, NaN], t); 2348 // } 2349 // if (hi === Infinity) { 2350 // this._insertPoint_v4(curve, [1, NaN, NaN], t); 2351 // this._insertPoint_v4(curve, [1, x, hi], t, true); 2352 // } 2353 // } 2354 // } 2355 } 2356 if (d_rgt < -d_thresh) { 2357 // Right branch very steep downwards -> add the maximum 2358 this._insertPoint_v4(curve, [1, x, hi], t); 2359 } else if (d_rgt > d_thresh) { 2360 // Right branch very steep upwards -> add the minimum 2361 this._insertPoint_v4(curve, [1, x, lo], t); 2362 } 2363 }, 2364 2365 /** 2366 * Number of equidistant points where the function is evaluated 2367 */ 2368 steps: 1021, //2053, // 1021, 2369 2370 /** 2371 * If the absolute maximum of the set of differences is larger than 2372 * criticalThreshold * median of these values, it is regarded as critical point. 2373 * @see JXG.Math.Plot._criticalInterval 2374 */ 2375 criticalThreshold: 1000, 2376 2377 plot_v4: function (curve, ta, tb, steps) { 2378 var i, 2379 // j, 2380 le, 2381 components, 2382 idx, 2383 comp, 2384 groups, 2385 g, 2386 start, 2387 ret, 2388 x_table, y_table, 2389 t, t1, t2, 2390 // good, 2391 // bad, 2392 // x_int, 2393 y_int, 2394 // degree_x, 2395 // degree_y, 2396 h = (tb - ta) / steps, 2397 Ypl = function (x) { 2398 return curve.Y(x, true); 2399 }, 2400 Ymi = function (x) { 2401 return -curve.Y(x, true); 2402 }, 2403 h2 = h * 0.5; 2404 2405 components = this.findComponents(curve, ta, tb, steps); 2406 for (idx = 0; idx < components.length; idx++) { 2407 comp = components[idx]; 2408 ret = this.differenceMethod(comp, curve); 2409 groups = ret[0]; 2410 x_table = ret[1]; 2411 y_table = ret[2]; 2412 2413 // degree_x = ret[3]; 2414 // degree_y = ret[4]; 2415 // if (degree_x >= 0) { 2416 // console.log("x polynomial of degree", degree_x); 2417 // } 2418 // if (degree_y >= 0) { 2419 // console.log("y polynomial of degree", degree_y); 2420 // } 2421 if (groups.length === 0 || groups[0].type !== "borderleft") { 2422 groups.unshift({ 2423 idx: 0, 2424 t: comp.t_values[0], 2425 x: comp.x_values[0], 2426 y: comp.y_values[0], 2427 type: "borderleft" 2428 }); 2429 } 2430 if (groups[groups.length - 1].type !== "borderright") { 2431 le = comp.t_values.length; 2432 groups.push({ 2433 idx: le - 1, 2434 t: comp.t_values[le - 1], 2435 x: comp.x_values[le - 1], 2436 y: comp.y_values[le - 1], 2437 type: "borderright" 2438 }); 2439 } 2440 2441 start = 0; 2442 for (g = 0; g <= groups.length; g++) { 2443 if (g === groups.length) { 2444 le = comp.len; 2445 } else { 2446 le = groups[g].idx - 1; 2447 } 2448 2449 // good = 0; 2450 // bad = 0; 2451 // Insert all uncritical points until next critical point 2452 for (i = start; i < le - 2; i++) { 2453 this._insertPoint_v4( 2454 curve, 2455 [1, comp.x_values[i], comp.y_values[i]], 2456 comp.t_values[i] 2457 ); 2458 // j = Math.max(0, i - 2); 2459 // Add more points in critical intervals 2460 if ( 2461 //degree_y === -1 && // No polynomial 2462 i >= start + 3 && 2463 i < le - 3 && // Do not do this if too close to a critical point 2464 y_table.length > 3 && 2465 Math.abs(y_table[2][i]) > 0.2 * Math.abs(y_table[0][i]) 2466 ) { 2467 t = comp.t_values[i]; 2468 h2 = h * 0.25; 2469 y_int = this.getInterval(curve, t, t + h); 2470 if (Type.isObject(y_int)) { 2471 if (y_table[2][i] > 0) { 2472 this._insertPoint_v4(curve, [1, t + h2, y_int.lo], t + h2); 2473 } else { 2474 this._insertPoint_v4( 2475 curve, 2476 [1, t + h - h2, y_int.hi], 2477 t + h - h2 2478 ); 2479 } 2480 } else { 2481 t1 = Numerics.fminbr(Ypl, [t, t + h]); 2482 t2 = Numerics.fminbr(Ymi, [t, t + h]); 2483 if (t1 < t2) { 2484 this._insertPoint_v4( 2485 curve, 2486 [1, curve.X(t1, true), curve.Y(t1, true)], 2487 t1 2488 ); 2489 this._insertPoint_v4( 2490 curve, 2491 [1, curve.X(t2, true), curve.Y(t2, true)], 2492 t2 2493 ); 2494 } else { 2495 this._insertPoint_v4( 2496 curve, 2497 [1, curve.X(t2, true), curve.Y(t2, true)], 2498 t2 2499 ); 2500 this._insertPoint_v4( 2501 curve, 2502 [1, curve.X(t1, true), curve.Y(t1, true)], 2503 t1 2504 ); 2505 } 2506 } 2507 // bad++; 2508 // } else { 2509 // good++; 2510 } 2511 } 2512 // console.log("GOOD", good, "BAD", bad); 2513 2514 // Handle next critical point 2515 if (g < groups.length) { 2516 //console.log("critical point / interval", groups[g]); 2517 2518 i = groups[g].idx; 2519 if (groups[g].type === "borderleft" || groups[g].type === "borderright") { 2520 this.handleBorder(curve, comp, groups[g], x_table, y_table); 2521 } else { 2522 this._seconditeration_v4(curve, comp, groups[g], x_table, y_table); 2523 } 2524 2525 start = groups[g].idx + 1 + 1; 2526 } 2527 } 2528 2529 le = comp.len; 2530 if (idx < components.length - 1) { 2531 this._insertPoint_v4(curve, [1, NaN, NaN], comp.right_t); 2532 } 2533 } 2534 }, 2535 2536 /** 2537 * Updates the data points of a parametric curve, plotVersion 4. This version is used if {@link JXG.Curve#plotVersion} is <tt>4</tt>. 2538 * @param {JXG.Curve} curve JSXGraph curve element 2539 * @param {Number} mi Left bound of curve 2540 * @param {Number} ma Right bound of curve 2541 * @returns {JXG.Curve} Reference to the curve object. 2542 */ 2543 updateParametricCurve_v4: function (curve, mi, ma) { 2544 var ta, tb, w2, bbox; 2545 2546 if (curve.xterm === "x") { 2547 // For function graphs we can restrict the plot interval 2548 // to the visible area +plus margin 2549 bbox = curve.board.getBoundingBox(); 2550 w2 = (bbox[2] - bbox[0]) * 0.3; 2551 // h2 = (bbox[1] - bbox[3]) * 0.3; 2552 ta = Math.max(mi, bbox[0] - w2); 2553 tb = Math.min(ma, bbox[2] + w2); 2554 } else { 2555 ta = mi; 2556 tb = ma; 2557 } 2558 2559 curve.points = []; 2560 2561 //console.log("--------------------"); 2562 this.plot_v4(curve, ta, tb, this.steps); 2563 2564 curve.numberPoints = curve.points.length; 2565 //console.log(curve.numberPoints); 2566 }, 2567 2568 //---------------------------------------------------------------------- 2569 // Plot algorithm alias 2570 //---------------------------------------------------------------------- 2571 2572 /** 2573 * Updates the data points of a parametric curve, alias for {@link JXG.Curve#updateParametricCurve_v2}. 2574 * This is needed for backwards compatibility, if this method has been 2575 * used directly in an application. 2576 * @param {JXG.Curve} curve JSXGraph curve element 2577 * @param {Number} mi Left bound of curve 2578 * @param {Number} ma Right bound of curve 2579 * @returns {JXG.Curve} Reference to the curve object. 2580 * 2581 * @see JXG.Curve#updateParametricCurve_v2 2582 */ 2583 updateParametricCurve: function (curve, mi, ma) { 2584 return this.updateParametricCurve_v2(curve, mi, ma); 2585 } 2586 }; 2587 2588 export default Mat.Plot; 2589