Changeset View
Changeset View
Standalone View
Standalone View
kmplot/parser.cpp
Show All 38 Lines | |||||
39 | #include <QList> | 39 | #include <QList> | ||
40 | 40 | | |||
41 | // standard c(++) includes | 41 | // standard c(++) includes | ||
42 | #include <stdlib.h> | 42 | #include <stdlib.h> | ||
43 | #include <stdio.h> | 43 | #include <stdio.h> | ||
44 | #include <math.h> | 44 | #include <math.h> | ||
45 | #include <assert.h> | 45 | #include <assert.h> | ||
46 | #include <cmath> | 46 | #include <cmath> | ||
47 | #include <tr1/cmath> | ||||
47 | #include <limits> | 48 | #include <limits> | ||
48 | #include <locale.h> | 49 | #include <locale.h> | ||
49 | 50 | | |||
50 | double Parser::m_radiansPerAngleUnit = 0; | 51 | double Parser::m_radiansPerAngleUnit = 0; | ||
51 | 52 | | |||
52 | /** | 53 | /** | ||
53 | * List of predefined functions. | 54 | * List of predefined functions. | ||
54 | * \note Some function names include other function names (e.g. "sinh" has the | 55 | * \note Some function names include other function names (e.g. "sinh" has the | ||
▲ Show 20 Lines • Show All 55 Lines • ▼ Show 20 Line(s) | 60 | { | |||
110 | // legendre | 111 | // legendre | ||
111 | {"P_0", QString(), legendre0}, // lengedre polynomial (n=0) | 112 | {"P_0", QString(), legendre0}, // lengedre polynomial (n=0) | ||
112 | {"P_1", QString(), legendre1}, // lengedre polynomial (n=1) | 113 | {"P_1", QString(), legendre1}, // lengedre polynomial (n=1) | ||
113 | {"P_2", QString(), legendre2}, // lengedre polynomial (n=2) | 114 | {"P_2", QString(), legendre2}, // lengedre polynomial (n=2) | ||
114 | {"P_3", QString(), legendre3}, // lengedre polynomial (n=3) | 115 | {"P_3", QString(), legendre3}, // lengedre polynomial (n=3) | ||
115 | {"P_4", QString(), legendre4}, // lengedre polynomial (n=4) | 116 | {"P_4", QString(), legendre4}, // lengedre polynomial (n=4) | ||
116 | {"P_5", QString(), legendre5}, // lengedre polynomial (n=5) | 117 | {"P_5", QString(), legendre5}, // lengedre polynomial (n=5) | ||
117 | {"P_6", QString(), legendre6}, // lengedre polynomial (n=6) | 118 | {"P_6", QString(), legendre6}, // lengedre polynomial (n=6) | ||
119 | | ||||
120 | // bessel | ||||
121 | {"J0", QString(), besselj0}, // the cylindrical Bessel function of the first kind (n=0) | ||||
118 | }; | 122 | }; | ||
119 | 123 | | |||
120 | VectorFunction Parser::vectorFunctions[ VectorCount ]= | 124 | VectorFunction Parser::vectorFunctions[ VectorCount ]= | ||
121 | { | 125 | { | ||
122 | {"min", min}, // minimum of a set of reals | 126 | {"min", min}, // minimum of a set of reals | ||
123 | {"max", max}, // maximum of a set of reals | 127 | {"max", max}, // maximum of a set of reals | ||
124 | {"mod", mod}, // l2 modulus of a set of reals | 128 | {"mod", mod}, // l2 modulus of a set of reals | ||
125 | }; | 129 | }; | ||
▲ Show 20 Lines • Show All 1264 Lines • ▼ Show 20 Line(s) | |||||
1390 | } | 1394 | } | ||
1391 | double legendre5( double x ) { | 1395 | double legendre5( double x ) { | ||
1392 | return (63*x*x*x*x*x - 70*x*x*x + 15*x)/8; | 1396 | return (63*x*x*x*x*x - 70*x*x*x + 15*x)/8; | ||
1393 | } | 1397 | } | ||
1394 | double legendre6( double x ) { | 1398 | double legendre6( double x ) { | ||
1395 | return (231*x*x*x*x*x*x - 315*x*x*x*x + 105*x*x - 5)/16; | 1399 | return (231*x*x*x*x*x*x - 315*x*x*x*x + 105*x*x - 5)/16; | ||
1396 | } | 1400 | } | ||
1397 | 1401 | | |||
1402 | // The following function code is adopted from Boost. | ||||
1403 | // Use, modification and distribution are subject to the | ||||
1404 | // Boost Software License, Version 1.0. (See accompanying file | ||||
1405 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | ||||
1406 | | ||||
1407 | double evaluate_rational(const double* num, const double* denom, const double& x_, int count) | ||||
1408 | { | ||||
1409 | double x(x_); | ||||
1410 | double s1, s2; | ||||
1411 | if (x <= 1) | ||||
1412 | { | ||||
1413 | s1 = num[count-1]; | ||||
1414 | s2 = denom[count-1]; | ||||
1415 | for (int i = count - 2; i >= 0; --i) | ||||
1416 | { | ||||
1417 | s1 *= x; | ||||
1418 | s2 *= x; | ||||
1419 | s1 += num[i]; | ||||
1420 | s2 += denom[i]; | ||||
1421 | } | ||||
1422 | } | ||||
1423 | else | ||||
1424 | { | ||||
1425 | x = 1 / x; | ||||
1426 | s1 = num[0]; | ||||
1427 | s2 = denom[0]; | ||||
1428 | for(unsigned i = 1; i < (unsigned)count; ++i) | ||||
1429 | { | ||||
1430 | s1 *= x; | ||||
1431 | s2 *= x; | ||||
1432 | s1 += num[i]; | ||||
1433 | s2 += denom[i]; | ||||
1434 | } | ||||
1435 | } | ||||
1436 | return s1 / s2; | ||||
1437 | } | ||||
1438 | | ||||
1439 | // The following function code is adopted from Boost. | ||||
1440 | // Use, modification and distribution are subject to the | ||||
1441 | // Boost Software License, Version 1.0. (See accompanying file | ||||
1442 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | ||||
1443 | | ||||
1444 | double besselj0( double x ) { | ||||
1445 | static const double P1[] = { | ||||
1446 | -4.1298668500990866786e+11, | ||||
1447 | 2.7282507878605942706e+10, | ||||
1448 | -6.2140700423540120665e+08, | ||||
1449 | 6.6302997904833794242e+06, | ||||
1450 | -3.6629814655107086448e+04, | ||||
1451 | 1.0344222815443188943e+02, | ||||
1452 | -1.2117036164593528341e-01 | ||||
1453 | }; | ||||
1454 | static const double Q1[] = { | ||||
1455 | 2.3883787996332290397e+12, | ||||
1456 | 2.6328198300859648632e+10, | ||||
1457 | 1.3985097372263433271e+08, | ||||
1458 | 4.5612696224219938200e+05, | ||||
1459 | 9.3614022392337710626e+02, | ||||
1460 | 1.0, | ||||
1461 | 0.0 | ||||
1462 | }; | ||||
1463 | static const double P2[] = { | ||||
1464 | -1.8319397969392084011e+03, | ||||
1465 | -1.2254078161378989535e+04, | ||||
1466 | -7.2879702464464618998e+03, | ||||
1467 | 1.0341910641583726701e+04, | ||||
1468 | 1.1725046279757103576e+04, | ||||
1469 | 4.4176707025325087628e+03, | ||||
1470 | 7.4321196680624245801e+02, | ||||
1471 | 4.8591703355916499363e+01 | ||||
1472 | }; | ||||
1473 | static const double Q2[] = { | ||||
1474 | -3.5783478026152301072e+05, | ||||
1475 | 2.4599102262586308984e+05, | ||||
1476 | -8.4055062591169562211e+04, | ||||
1477 | 1.8680990008359188352e+04, | ||||
1478 | -2.9458766545509337327e+03, | ||||
1479 | 3.3307310774649071172e+02, | ||||
1480 | -2.5258076240801555057e+01, | ||||
1481 | 1.0 | ||||
1482 | }; | ||||
1483 | static const double PC[] = { | ||||
1484 | 2.2779090197304684302e+04, | ||||
1485 | 4.1345386639580765797e+04, | ||||
1486 | 2.1170523380864944322e+04, | ||||
1487 | 3.4806486443249270347e+03, | ||||
1488 | 1.5376201909008354296e+02, | ||||
1489 | 8.8961548424210455236e-01 | ||||
1490 | }; | ||||
1491 | static const double QC[] = { | ||||
1492 | 2.2779090197304684318e+04, | ||||
1493 | 4.1370412495510416640e+04, | ||||
1494 | 2.1215350561880115730e+04, | ||||
1495 | 3.5028735138235608207e+03, | ||||
1496 | 1.5711159858080893649e+02, | ||||
1497 | 1.0 | ||||
1498 | }; | ||||
1499 | static const double PS[] = { | ||||
1500 | -8.9226600200800094098e+01, | ||||
1501 | -1.8591953644342993800e+02, | ||||
1502 | -1.1183429920482737611e+02, | ||||
1503 | -2.2300261666214198472e+01, | ||||
1504 | -1.2441026745835638459e+00, | ||||
1505 | -8.8033303048680751817e-03 | ||||
1506 | }; | ||||
1507 | static const double QS[] = { | ||||
1508 | 5.7105024128512061905e+03, | ||||
1509 | 1.1951131543434613647e+04, | ||||
1510 | 7.2642780169211018836e+03, | ||||
1511 | 1.4887231232283756582e+03, | ||||
1512 | 9.0593769594993125859e+01, | ||||
1513 | 1.0 | ||||
1514 | }; | ||||
1515 | static const double x1 = 2.4048255576957727686e+00, | ||||
1516 | x2 = 5.5200781102863106496e+00, | ||||
1517 | x11 = 6.160e+02, | ||||
1518 | x12 = -1.42444230422723137837e-03, | ||||
1519 | x21 = 1.4130e+03, | ||||
1520 | x22 = 5.46860286310649596604e-04; | ||||
1521 | | ||||
1522 | double value, factor, r, rc, rs; | ||||
1523 | | ||||
1524 | if (x < 0) | ||||
1525 | { | ||||
1526 | x = -x; // even function | ||||
1527 | } | ||||
1528 | if (x == 0) | ||||
1529 | { | ||||
1530 | return 1; | ||||
1531 | } | ||||
1532 | if (x <= 4) // x in (0, 4] | ||||
1533 | { | ||||
1534 | double y = x * x; | ||||
1535 | r = evaluate_rational(P1, Q1, y, 7); | ||||
1536 | factor = (x + x1) * ((x - x11/256) - x12); | ||||
1537 | value = factor * r; | ||||
1538 | } | ||||
1539 | else if (x <= 8.0) // x in (4, 8] | ||||
1540 | { | ||||
1541 | double y = 1 - (x * x)/64; | ||||
1542 | r = evaluate_rational(P2, Q2, y, 8); | ||||
1543 | factor = (x + x2) * ((x - x21/256) - x22); | ||||
1544 | value = factor * r; | ||||
1545 | } | ||||
1546 | else // x in (8, \infty) | ||||
1547 | { | ||||
1548 | double y = 8 / x; | ||||
1549 | double y2 = y * y; | ||||
1550 | rc = evaluate_rational(PC, QC, y2, 6); | ||||
1551 | rs = evaluate_rational(PS, QS, y2, 6); | ||||
1552 | factor = 1/ sqrt(M_PI) / sqrt(x); | ||||
1553 | double sx = sin(x); | ||||
1554 | double cx = cos(x); | ||||
1555 | value = factor * (rc * (cx + sx) - y * rs * (sx - cx)); | ||||
1556 | } | ||||
1557 | | ||||
1558 | return value; | ||||
1559 | } | ||||
1560 | | ||||
1398 | double sign(double x) | 1561 | double sign(double x) | ||
1399 | { | 1562 | { | ||
1400 | if(x<0.) | 1563 | if(x<0.) | ||
1401 | return -1.; | 1564 | return -1.; | ||
1402 | else if(x>0.) | 1565 | else if(x>0.) | ||
1403 | return 1.; | 1566 | return 1.; | ||
1404 | return 0.; | 1567 | return 0.; | ||
1405 | } | 1568 | } | ||
▲ Show 20 Lines • Show All 385 Lines • Show Last 20 Lines |