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