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