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