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