qpint1d.pro 60.6 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677
;+
; NAME:
;   QPINT1D
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   One dimensional numerical integration of IDL function or expression
;
; MAJOR TOPICS:
;   Numerical Analysis.
;
; CALLING SEQUENCE:
;   value = QPINT1D(FUNCT, A, B, [ PRIVATE, /EXPRESSION, FUNCTARGS=, 
;                   ERROR=error, NFEV=nfev, STATUS=status, NSUBINTERVALS=nsub,
;                   EPSABS=, EPSREL=, LIMIT=, BREAKPOINTS=, NPOINTS=
;                   /SYMMETRIC, SYM_AXIS= ] )
;
; DESCRIPTION:
;
;  QPINT1D adaptively calculates an approximation result to a given
;  definite integral
;
;       result = Integral[ f(x) dx ] over [a,b]
;
;  hopefully satisfying a constraint on the accuracy of the solution.
;  QPINT1D is based on the QUADPACK fortran package originally by
;  Piessens, de Doncker, Ueberhuber and Kahaner (and implements
;  equivalents to the QAGSE, QAGPE, QAGIE, and DQKxx fortran routines).
;
;  The returned result is intended to satisfy the following claim for
;  accuracy: ABS(result-value) LE MAX([epsabs, epsrel*ABS(value)]),
;  where VALUE is the true value of the integral, and EPSABS and
;  EPSREL are the absolute and relative error tolerances defined
;  below.  An estimate of the error is returned in the ERROR keyword.
;  Either A or B may be finite or infinite (i.e., an improper
;  integral).
;
;  QPINT1D is "adaptive" in the sense that it locates regions of the
;  integration interval which contain the highest error, and focusses
;  its efforts on those regions.  The algorithm locates these regions
;  by successively bisecting the starting interval.  Each subinterval
;  is assigned an error estimate, and the region with the largest
;  error estimate is subdivided further, until each subinterval
;  carries approximately the same amount of error.  Convergence of the
;  procedure may be accelerated by the Epsilon algorithm due to Wynn.
;
;  The estimate of the integral and the estimate of the error in each
;  subinterval are computed using Gauss Kronrod quadrature.
;  Integrators based on the 15-, 21-, 31-, 41-, 51- and 61-point
;  Gauss-Kronrod rule are available, and selected using the NPOINTS
;  keyword.  Generally, the more points the greater the precision,
;  especially for rapidly varying functions.  However the default
;  value of 21 is often sufficient, especially because of the adaptive
;  nature of QPINT1D.
;
;  In the following sections the requirements for the form of the
;  integrand are established.  Also, a description of how QPINT1D
;  handles singularities and discontinuities is presented.
;
; INTEGRAND is a FUNCTION
;
;  The integrand can be specified in two forms, either as a standard
;  IDL function, or as an IDL expression.  If integrating a function,
;  then the FUNCT should be a string naming the function.  The
;  function must be declared as following:
;
;    FUNCTION MYFUNCT, X, P, KEYWORDS=...
;      RETURN, (compute function of X and P)
;    END
;
;  The function must accept at least one, but optionally two,
;  parameters.  The first, 'X', is a vector of abcissae where the
;  function is to be computed.  The function must return the same
;  number of function values as abcissae passed.  The second
;  positional parameter, 'P', is a purely optional PRIVATE parameter
;  as described below.  MYFUNCT may accept more positional parameters,
;  but QPINT1D will not use them.  The difference between X and P is
;  that X is the variable of integration, while P contains any other
;  information expected to remain essentially constant over the
;  integration.
;
; INTEGRAND is an EXPRESSION
;
;  The integrand can also be specfied as an IDL expression and setting
;  the EXPRESSION keyword.  Any expression that can accept a vector of
;  abcissae named 'X' and produce a corresponding vector of output is
;  a valid expression.  Here is an example:
;
;    RESULT = QPINT1D('X^2 * EXP(-X)', /EXPRESSION, 0D, 10D)
;
;  It is important to note that the variable of integration must
;  always be named 'X', and the expression must be vectorizable.  The
;  expression may also use the PRIVATE data, and as above, it would be
;  referred to according to the variable 'P'.  For example, if the
;  exponential decay constant is parameterized by PRIVATE(0), then the
;  expression would be:
;
;    RESULT = QPINT1D('X^2 * EXP(-X/P(0))', /EXPRESSION, 0D, 10D, PRIVATE)
;
;  The user is solely responsible for defining and using the PRIVATE
;  data.  QPINT1D does not access or modify PRIVATE / P; it only
;  passes it on to the user routine for convenience.
;
; IMPROPER INTEGRALS and DISCONTINUITIES
;
;  QPINT1D computes improper integrals, as well as integrands with
;  discontinuities or singularities.  
;
;  Improper integrals are integrals where one or both of the limits of
;  integration are "infinity."  (Formally, these integrals are defined
;  by taking the limit as the integration limit tends to infinity).
;  QPINT1D handles a small class of such integrals, generally for
;  integrands that are convergent and monotonic (i.e.,
;  non-oscillatory, and falling off as 1/ABS(X)^2 or steeper).  Such
;  integrals are handled by a transformation of the original interval
;  into the interval [0,1].
;
;  Integrals from negative infinity to positive infinity are done in
;  two subintervals.  By default the interval is split at X EQ 0,
;  however this can be controlled by using the SYM_AXIS keyword.
;  Users should note that if the first subinterval fails the second is
;  not attempted, and thus the return value VALUE should not be
;  trusted in those cases.
;
;  Infinite integration limits are specified by using the standard
;  values !VALUES.F_INFINITY or !VALUES.D_INFINITY.  No other special
;  invocation syntax is required.
;
;  The integration routine is able to handle integrands which have
;  integrable singularities at the endpoints.  For example, the
;  integral:
;
;    RESULT = QPINT1D('2*sqrt((1-x)/(1+x))/(1-x^2)', 0.0d, 1d, /expr)
;
;  has a singularity at a value of X EQ 1.  Still, the singularity is
;  integrable, and the value returned is a correct value of 2.
;
;  If known singularities are present within the interval of
;  integration, then users should pass the BREAKPOINTS keyword to list
;  the locations of these points.  QPINT1D will then integrate each
;  subinterval separately, while still maintaining an overall error
;  budget.
;
;  If known discontinuities exist in the integrand, then the user may
;  additionally list those points using the BREAKPOINTS keyword.
;
;  It should be noted that the algorithm used is different, depending
;  on whether the BREAKPOINTS keyword has been specified or not (this
;  is the difference between the QAGSE vs. QAGPE routines in the
;  original FORTRAN).  The algorithm *without* BREAKPOINTS is
;  generally thought to be more precise than *with*.  Thus, it may be
;  worth splitting the original integration interval manually and
;  invoking QPINT1D without BREAKPOINTS.
;
;
; INPUTS:
;
;  FUNCT - by default, a scalar string containing the name of an IDL
;          function to be integrated.  See above for the formal
;          definition of MYFUNCT.  (No default).
;
;          If the EXPRESSION keyword is set, then FUNCT is a scalar
;          string containing an IDL expression to be evaluated, as
;          described above.
;
;  A, B - a scalar number indicating the lower and upper limits of the
;         interval of integration (i.e., [A, B] is the interval of
;         integration).
;
;  PRIVATE - any optional variable to be passed on to the function to
;            be integrated.  For functions, PRIVATE is passed as the
;            second positional parameter; for expressions, PRIVATE can
;            be referenced by the variable 'P'.   QPINT1D does not
;            examine or alter PRIVATE.
;
; RETURNS:
;
;   The value of the integral.  If either A or B are double precision,
;   then the integral is computed in double precision; otherwise the
;   result is returned in single precision floating point.
;
; KEYWORD PARAMETERS:
;
;   BREAKPOINTS - an array of numbers specifying points within the
;                 integration interval where the integrand is
;                 discontinuous or singular.  Out of bounds points are
;                 ignored.
;                 Default: undefined, i.e., no such points
;
;   EPSABS - a scalar number, the absolute error tolerance requested
;            in the integral computation.  If less than or equal to
;            zero, then the value is ignored.
;            Default: 0
;
;   EPSREL - a scalar number, the relative (i.e., fractional) error
;            tolerance in the integral computation.  If less than or
;            equal to zero, then the value is ignored.
;            Default: 1e-4 for float, 1d-6 for double
;
;   EXPRESSION - if set, then FUNCT is an IDL expression.  Otherwise,
;                FUNCT is an IDL function.
;
;   ERROR - upon return, this keyword contains an estimate of the
;           error in the computation.
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by FUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.  By default, no extra parameters
;               are passed to the user-supplied function.
;
;   LIMIT - a scalar, the maximum number of subintervals to create
;           before terminating execution.  Upon return, a STATUS value
;           of 1 indicates such an overflow occurred.
;           Default: 100
;
;   NFEV - upon return, this keyword contains the number of function
;          calls executed (i.e., the number of abcissae).
;
;   NPOINTS - a scalar, the number of Gauss Kronrod points to use in
;             computing the integral over a subinterval.  A larger
;             number of points can in principle increase the precision
;             of the integral, but also makes the computation take
;             longer.  Possible values are 15, 21, 31, 41, 51, and 61.
;             NPOINTS is rounded up to the next nearest available set,
;             with a maximum of 61.
;             Default: 21
;
;   NSUBINTERVALS - upon return, this keyword contains the number of
;                   subintervals the integration interval was divided
;                   into.
;
;   STATUS - upon return, the status of the integration operation is
;            returned in this keyword as an integer value.  A value of
;            zero indicates success; otherwise an abnormal condition
;            has occurred and the returned value should be considered
;            erroneous or less reliable according to STATUS:
;
;                 any negative number - outright failure (reserved for
;                                       future use).
;
;                 -1 - the input parameters are invalid, because
;                      epsabs LE 0 and epsrel LT max([50*EPS,0.5d-28]),
;                      where EPS is the machine precision, or if LIMIT
;                      is smaller than the number of BREAKPOINTS.
;
;                  0 - success.
;
;                  1 - maximum number of subdivisions allowed has been
;                      achieved. One can allow more subdivisions by
;                      increasing the value of limit (and taking the
;                      according dimension adjustments into
;                      account). However, if this yields no
;                      improvement it is advised to analyze the
;                      integrand in order to determine the integration
;                      difficulties.  If the position of a local
;                      difficulty can be determined (i.e. singularity,
;                      discontinuity within the interval), it should
;                      be supplied to the routine as an element of the
;                      vector BREAKPOINTS.
;
;                  2 - The occurrence of roundoff error is detected,
;                      which prevents the requested tolerance from
;                      being achieved.  The error may be
;                      under-estimated.
;
;                  3 - Extremely "bad" integrand behaviour occurs at
;                      some points of the integration interval.
;
;                  4 - The algorithm does not converge.  Roundoff
;                      error is detected in the extrapolation table.
;                      It is presumed that the requested tolerance
;                      cannot be achieved, and that the returned
;                      result is the best which can be obtained.
;
;                  5 - The integral is probably divergent, or only
;                      slowly convergent.  It must be noted that
;                      divergence can occur with any other value of
;                      ier GT 0.
;
;   SYM_AXIS - a scalar number, the bisection point of the real line
;              for improper integrals from negative infinity to
;              positive infinity.  Otherwise ignored.
;              Default: 0.
;
;   
; EXAMPLES:
;
;  Shows how function and expression can be used for exponential
;  integrand:
;
;    IDL> print, qpint1d('EXP(X)', 0D, 10D, /expr)
;           22025.466
;    IDL> print, qpint1d('EXP', 0D, 10D)
;           22025.466
;
;  Normal definite integral, and then parameterized using a PRIVATE
;  value of 2.
;    IDL> print, qpint1d('X^2*EXP(-X)', 0D, 10D, /expr)
;           1.9944612
;    IDL> print, qpint1d('X^2*EXP(-X/P(0))', 0D, 10D, 2D, /expr)
;           14.005568
;
;  Improper integrals of the gaussian function
;    IDL> inf = !values.d_infinity
;    IDL> print, qpint1d('EXP(-X^2)', 0D, +inf, 2D, /expr)
;          0.88622693
;    IDL> print, qpint1d('EXP(-X^2)', -inf, +inf, 2D, /expr), sqrt(!dpi)
;           1.7724539       1.7724539
;  The second integral shows the comparison to the analytic value of
;  SQRT(!DPI).
;
; COMMON BLOCKS:
;
;   COMMON QPINT1D_MACHAR
;   COMMON QPINT1D_PROFILE_COMMON
;   COMMON QPINT1D_QKEVAL_COMMON
;
;   These common blocks are used internally only and should not be
;   accessed or modified.
;
; REFERENCES:
;
;    R. Piessens, E. deDoncker-Kapenga, C. Uberhuber, D. Kahaner
;    Quadpack: a Subroutine Package for Automatic Integration
;    Springer Verlag, 1983.    Series in Computational Mathematics v.1
;    515.43/Q1S  100394Z
;
;    Netlib repository: http://www.netlib.org/quadpack/
;
;    SLATEC Common Mathematical Library, Version 4.1, July 1993
;    a comprehensive software library containing over
;    1400 general purpose mathematical and statistical routines
;    written in Fortran 77.  (http://www.netlib.org/slatec/)
;
; MODIFICATION HISTORY:
;   Written, Feb-Jun, 2001, CM
;   Documented, 04 Jun, 2001, CM
;   Add usage message, error checking, 15 Mar 2002, CM
;   Correct usage message, 28 Apr 2002, CM
;   More error checking when user EXPRession fails, 10 Jun 2009, CM
;
;  $Id: qpint1d.pro,v 1.14 2009/06/10 22:03:34 craigm Exp $
;-
; Copyright (C) 2001, 2002, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

FORWARD_FUNCTION qpint1d_eval, qpint1d

; --------------------------------------------------------------
; Some functions used in testing the algorithm
;
; Test functions
; .COMP
; FUNCTION BESTEST, X, N
;  RETURN, BESELJ(ABS(X-10)+0.001, N)
; END
; print, qpint1d('BESTEST', .1d, 35d, 1L, nev=nev, nsub=nsub, status=ier)

; .comp
; function invsq, x
;   return, 1/(abs(x-10)+0.0001)^2
; end
; print, qpint1d('INVSQ', 0d, 50d, nev=nev, nsub=nsub, status=ier, error=err)

; .comp 
; function broken, x
;   f = x*0
;   wh = where(x LT 0, ct)
;   if ct GT 0 then f(wh) = exp(x(wh))
;   wh = where(x GE 0 AND x LT 10, ct)
;   if ct GT 0 then f(wh) = x(wh)^2
;   wh = where(x GE 10, ct)
;   if ct GT 0 then f(wh) = exp(-x(wh))
;   return, f
; end
; print, qpint1d('BROKEN', -10d, 50d, nev=nev, nsub=nsub, status=ier, error=err) & print, err, ier, nsub


;; Function used in testing speed
pro qpint1d_profile, clear=clear
  common qpint1d_profile_common, profvals
;  if n_elements(profvals) EQ 0 OR keyword_set(clear) then $
;    profvals = {functime: 0D, tottime: 0D, srttime: 0D, acctime: 0D, $
;                qktime: 0D}
end

;; Following are machine constants that can be loaded once.  I have
;; found that bizarre underflow messages can be produced in each call
;; to MACHAR(), so this structure minimizes the number of calls to
;; one.

pro qpint1d_setmachar, double=isdouble
  common qpint1d_machar, qpint1d_machar_vals

  ;; In earlier versions of IDL, MACHAR itself could produce a load of
  ;; error messages.  We try to mask some of that out here.
  if (!version.release) LT 5 then dummy = check_math(1, 1)

  mch = 0.
  mch = machar(double=keyword_set(isdouble))
  dmachep = mch.eps
  dmaxnum = mch.xmax
  dminnum = mch.xmin
  dmaxlog = alog(mch.xmax)
  dminlog = alog(mch.xmin)
  if keyword_set(isdouble) then $
    dmaxgam = 171.624376956302725D $
  else $
    dmaxgam = 171.624376956302725
  drdwarf = sqrt(dminnum*1.5) * 10
  drgiant = sqrt(dmaxnum) * 0.1

  qpint1d_machar_vals = {machep: dmachep, maxnum: dmaxnum, minnum: dminnum, $
                       maxlog: dmaxlog, minlog: dminlog, maxgam: dmaxgam, $
                       rdwarf: drdwarf, rgiant: drgiant}

  if (!version.release) LT 5 then dummy = check_math(0, 0)

  return
end


;; --------------------------------------------------------
;; Main workhorse routine
pro qpint1d_qagse, f, a0, b0, result, abserr, private, functargs=fa, $
                   epsabs=epsabs, epsrel=epsrel, npoints=npts0, $
                   status=ier, limit=limit, neval=neval, nsubintervals=last, $
                   breakpoints=bpoints0, isdouble=isdouble, $
                   alist=alist, blist=blist, rlist=rlist, $
                   elist=elist, iord=iord
                   

  ;; Derived from QUADPACK QAGSE
  ;;***PURPOSE  The routine calculates an approximation result to a given
  ;;            definite integral I = Integral of F over (A,B),
  ;;            hopefully satisfying following claim for accuracy
  ;;            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  ;;***LIBRARY   SLATEC (QUADPACK)

  common qpint1d_machar, machvals
  r1mach1 = machvals.minnum
  r1mach2 = machvals.maxnum
  r1mach4 = machvals.machep

  ;; Establish the precision we will be computing in, base on the
  ;; integration limits.
  a = a0
  b = b0
  zero = a*0.
  if n_elements(isdouble) EQ 0 then begin
      sz = size(zero)
      if sz(sz(0)+1) EQ 5 then isdouble = 1 else isdouble = 0
  endif
  isdouble = isdouble(0)

  ;; Establish default values
  if n_elements(neval) EQ 0 then neval = 0L
  if n_elements(limit) EQ 0 then limit = 100L
  if n_elements(epsabs) EQ 0 then epsabs = zero
  if n_elements(epsrel) EQ 0 then epsrel = zero + 1e-6
  if n_elements(npts0) EQ 0 then npoints = 21L $
  else                           npoints = floor(npts0(0))
  rlist2 = make_array(52, value=zero)
  res3la = rlist2(0:2)

  ;; These are the work arrays which manage the subintervals
  alist = make_array(limit, value=zero)
  blist = alist
  rlist = alist
  elist = alist
  iord = lonarr(limit)
  level = iord

  epmach = r1mach4

;            TEST ON VALIDITY OF PARAMETERS
  ier = 0L
  neval = 0
  last = 0
  result = zero
  abserr = zero

  ;; Extract the valid breakpoints
  nbp0 = n_elements(bpoints0)
  nbp = nbp0
  if nbp GT 0 then begin
      abmin = min([a,b], max=abmax)
      wh = where(bpoints0 GT abmin(0) AND bpoints0 LT abmax(0), nbp)
      if nbp GT 0 then begin
          ;; Valid breakpoints were found.  Sort them into place.
          bpoints = bpoints0(wh)
          if a LE b then bpoints = [a, bpoints(sort(bpoints)),  b] $
          else           bpoints = [a, bpoints(sort(-bpoints)), b]
          nbp = nbp + 2
          if nbp GE limit then goto, INPUT_ERROR
      endif
  endif
  if nbp EQ 0 then begin
      ;; By default, we have two breakpoints, namely the start and
      ;; stop of the interval.
      bpoints = [a, b]
      nbp = 2L
  endif

  ;; Search for infinite bounds
  dirsign = 1
  if finite(b) EQ 0 then begin
      bound = a
      dirsign = 2*(b GT 0) - 1
      a = zero
      b = zero + 1

      ;; Transform the breakpoints into the [0,1] interval
      bpoints = rotate(1/((bpoints-bound(0))*dirsign + 1),2)
  endif

  alist(0) = a
  blist(0) = b
  rlist(0) = zero
  elist(0) = zero
  xnum = (a(0)*b(0)*0.)
  sz = size(xnum)
  if sz(sz(0)+1) EQ 5 then xlim = 0.5d-28 else xlim = 0.5d-14

  if (epsabs LE 0 AND epsrel LT max([0.5e2*epmach,xlim])) then begin
      INPUT_ERROR:
      ier = 6L
      goto, LAB999
  endif

  nint = nbp-1
  npts = nbp-2
  npts2 = nbp

;           FIRST APPROXIMATION TO THE INTEGRAL
  uflow = r1mach1
  oflow = r1mach2
  resabs = zero
  a1 = bpoints(0)
  ndin = iord * 0
  for i = 0L,nint-1 do begin
      b1 = bpoints(i+1)

      qpint1d_qkeval, f,a1,b1,area1,error1,defabs,resa, private, $
        neval=neval, functargs=fa, inflow=bound, dirsign=dirsign, $
        isdouble=isdouble, npoints=npoints

      abserr = abserr+error1
      result = result+area1
      ndin(i) = (error1 EQ resa AND error1 NE zero)
      resabs = resabs+defabs
      level(i) = 0
      elist(i) = error1
      alist(i) = a1
      blist(i) = b1
      rlist(i) = area1
      iord(i) = i+1
      a1 = b1
  endfor

  if nbp0 GT 0 then begin
      wh = where(ndin EQ 1, ct)
      if ct GT 0 then elist(wh) = abserr
  endif
  errsum = total(elist(0:nint-1))

;           TEST ON ACCURACY.
  last = nint
  dres = abs(result)
  errbnd = max([epsabs,epsrel*dres])
  if (abserr LE 1.0e+02*epmach*resabs) AND (abserr GT errbnd) then $
    ier = 2L

  ;; Sort in descending order
  iord(0:nint-1) = sort(-elist(0:nint-1)) + 1

  if (limit LT npts2) then ier = 1L
  if (ier NE 0 OR abserr LE errbnd) then goto, LAB999

;           INITIALIZATION
  rlist2(1-1) = result
  maxerr = iord(0)
  errmax = elist(maxerr-1)
  area = result
  nrmax = 1L
  nres = 0L
  if nbp0 EQ 0 then begin
      numrl2 = 2L
      errsum = abserr
  endif else begin
      numrl2 = 1L
      erlarg = errsum
  endelse
  ktmin = 0
  extrap = 0
  noext = 0
  ertest = errbnd
  levmax = 1L
  iroff1 = 0
  iroff2 = 0
  iroff3 = 0
  ierro = 0
  abserr = oflow
  ksgn = -1
  if (dres GE (0.1E+01-0.5E+02*epmach)*resabs) then ksgn = 1
  
;           MAIN DO-LOOP
  for last = npts2, limit do begin

;           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST
;           ERROR ESTIMATE.
      levcur = level(maxerr)+1
      a1 = alist(maxerr-1)
      b1 = 0.5e+00*(alist(maxerr-1)+blist(maxerr-1))
      a2 = b1
      b2 = blist(maxerr-1)
      erlast = errmax
      qpint1d_qkeval, f,a1,b1,area1,error1,resa,defab1, private, $
        neval=neval, functargs=fa, inflow=bound, dirsign=dirsign, $
        isdouble=isdouble, npoints=npoints
      qpint1d_qkeval, f,a2,b2,area2,error2,resa,defab2, private, $
        neval=neval, functargs=fa, inflow=bound, dirsign=dirsign, $
        isdouble=isdouble, npoints=npoints

;           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
;           AND ERROR AND TEST FOR ACCURACY.
      area12 = area1+area2
      erro12 = error1+error2
      errsum = errsum+erro12-errmax
      area = area+(area12-rlist(maxerr-1))
      if NOT (defab1 EQ error1 OR defab2 EQ error2) then begin
          
          if NOT (abs(rlist(maxerr-1)-area12) GT 0.1e-04*abs(area12) $
                  OR erro12 LT 0.99e+00*errmax) then begin
              if (extrap) then iroff2 = iroff2+1
              if (NOT extrap) then iroff1 = iroff1+1
          endif
          if(last GT 10 AND erro12 GT errmax) then iroff3 = iroff3+1
      endif
      level(maxerr-1) = levcur
      level(last-1) = levcur
      rlist(maxerr-1) = area1
      rlist(last-1) = area2
      errbnd = max([epsabs,epsrel*abs(area)])

;           TEST FOR ROUNDOFF ERROR AND EVENTUALLY
;           SET ERROR FLAG.
      if (iroff1+iroff2 GE 10 OR iroff3 GE 20) then ier = 2L
      if (iroff2 GE 5) then ierro = 3

;           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
;           SUBINTERVALS EQUALS LIMIT.
      if last EQ limit then ier = 1L
      
;           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
;           AT A POINT OF THE INTEGRATION RANGE.
      if (max([abs(a1),abs(b2)]) LE (0.1e+01+0.1e+03*epmach)* $
          (abs(a2)+0.1e+04*uflow)) then ier = 4L

;           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
      if (error2 LE error1) then begin ;; IDL OK
          alist(last-1) = a2
          blist(maxerr-1) = b1
          blist(last-1) = b2
          elist(maxerr-1) = error1
          elist(last-1) = error2
      endif else begin
          alist(maxerr-1) = a2
          alist(last-1) = a1
          blist(last-1) = b1
          rlist(maxerr-1) = area2
          rlist(last-1) = area1
          elist(maxerr-1) = error2
          elist(last-1) = error1
      endelse

;           CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
;           IN THE LIST OF ERROR ESTIMATES AND SELECT THE
;           SUBINTERVAL WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE
;           BISECTED NEXT).

      qpint1d_qpsrt, limit,last,maxerr,errmax,elist,iord,nrmax

; ***JUMP OUT OF DO-LOOP
      if (errsum LE errbnd) then goto, LAB115
; ***JUMP OUT OF DO-LOOP
      if (ier NE 0) then goto, LAB100
      if (nbp0 EQ 0 AND last EQ 2) then goto, LAB80
      if (noext) then goto, LAB90
      erlarg = erlarg-erlast
      if (nbp0 EQ 0) then begin
          if (abs(b1-a1) GT small) then erlarg = erlarg+erro12 
      endif else begin
          if (levcur+1 LE levmax) then erlarg = erlarg+erro12
      endelse
      if NOT (extrap) then begin
          bool = 0
          if nbp0 EQ 0 then begin
              if (abs(blist(maxerr-1)-alist(maxerr-1)) GT small) then $
                goto, LAB90
          endif else begin
              if (level(maxerr-1)+1 LE levmax) then goto, LAB90
          endelse
          extrap = 1
          nrmax = 2L
      endif
      
      if NOT (ierro EQ 3 OR erlarg LE ertest) then begin
;
;           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
;           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS
;           OVER THE LARGER INTERVALS (ERLARG) AND PERFORM
;           EXTRAPOLATION.

          id = nrmax
          jupbnd = last
          if (last GT (2+limit/2)) then jupbnd = limit+3-last
          for k = id, jupbnd do begin
              maxerr = iord(nrmax-1)
              errmax = elist(maxerr-1)
; ***JUMP OUT OF DO-LOOP
              bool = 0
              if nbp0 EQ 0 then begin
                  if (abs(blist(maxerr-1)-alist(maxerr-1)) GT small) then $
                    goto, LAB90
              endif else begin
                  if (level(maxerr-1)+1 LE levmax) then goto, LAB90
              endelse
              nrmax = nrmax+1
          endfor
      endif

;           PERFORM EXTRAPOLATION.
      numrl2 = numrl2+1
      rlist2(numrl2-1) = area
      if (nbp0 GT 0 AND numrl2 LE 2) then goto, LAB72
      qpint1d_qelg, numrl2,rlist2,reseps,abseps,res3la,nres
      ktmin = ktmin+1
      if (ktmin GT 5 AND abserr LT 0.1e-02*errsum) then ier = 5L
      if NOT (abseps GE abserr) then begin
          ktmin = 0
          abserr = abseps
          result = reseps
          correc = erlarg
          ertest = max([epsabs,epsrel*abs(reseps)])
; ***JUMP OUT OF DO-LOOP
          if (abserr LE ertest) then goto, LAB100
      endif

;           PREPARE BISECTION OF THE SMALLEST INTERVAL.
      if (numrl2 EQ 1) then noext = 1
      if (ier EQ 5) then goto, LAB100

      LAB72:
      maxerr = iord(0)
      errmax = elist(maxerr-1)
      nrmax = 1L
      extrap = 0
      if nbp0 EQ 0 then small = small*0.5e+00
      levmax = levmax + 1
      erlarg = errsum
      goto, LAB90

      LAB80:
      small = abs(b-a)*0.375e+00
      erlarg = errsum
      ertest = errbnd
      rlist2(2-1) = area
      
      LAB90:
  endfor
      
  LAB100:
;           SET FINAL RESULT AND ERROR ESTIMATE.
  if abserr EQ oflow then goto, LAB115
  if ier+ierro EQ 0 then goto, LAB110
  if ierro EQ 3 then abserr = abserr+correc
  if ier EQ 0 then ier = 3L
  if NOT (result NE zero AND area NE zero) then begin
      if abserr GT errsum then goto, LAB115
      if area EQ 0 then goto, LAB130
      goto, LAB110
  endif
  if(abserr/abs(result) GT errsum/abs(area)) then goto, LAB115

;           TEST ON DIVERGENCE.
  LAB110:
  if (ksgn EQ (-1) AND max([abs(result),abs(area)]) LE $
      resabs*0.1e-01) then goto, LAB130
  if (0.1e-01 GT (result/area) OR (result/area) GT 0.1e+03 $
      OR errsum GT abs(area)) then ier = 6L
  goto, LAB130

  LAB115:
  result = total(rlist(0:last-1))
  abserr = errsum

  LAB130:
  if (ier GT 2) then ier = ier-1

  LAB999:
  result = result * dirsign
  return
end

;; Originally DQK21.F, this routine performs Gauss-Kronrod quadrature
;; using 15, 21, 31, 41, 51 or 61 points.
pro qpint1d_qkeval, f, a, b, result, abserr, resabs, resasc, priv, $  
                  functargs=fa, neval=neval, reset=reset, isdouble=isdouble, $
                  inflow=bound, dirsign=dirsign, npoints=npts0

  ;; Derived from QUADPACK QK21-QK61, and QK15I
  ;;***PURPOSE  To compute I = Integral of F over (A,B), with error
  ;;                           estimate
  ;;                       J = Integral of ABS(F) over (A,B)
  ;;***LIBRARY   SLATEC (QUADPACK)

  common qpint1d_qkeval_common, wg, wgk, xgk, ig, nptsreq, nptsact, prec

  if n_elements(npts0) EQ 0 then npts0 = 21L
  if n_elements(isdouble) EQ 0 then isdouble = 1

  ;; Determine the number of points "requested", versus the number
  ;; "actually" granted.
  if n_elements(nptsreq) EQ 0 OR keyword_set(reset) then begin
      nptsreq = 0L
      prec = -1L
  endif
  if npts0(0) NE nptsreq OR isdouble(0) NE prec then begin
      ;; If this is the first time around, or if we need to use a
      ;; different set of GK points, then request the new set of
      ;; points from GKWEIGHTS.
      nptsreq = npts0(0)
      prec = isdouble(0) EQ 1
      qpint1d_gkweights, wg, wgk, xgk, ig, nptsreq, nptsact, prec
      if keyword_set(reset) then return
  endif

  common qpint1d_machar, machvals
  r1mach1 = machvals.minnum
  r1mach2 = machvals.maxnum
  r1mach4 = machvals.machep
  common qpint1d_profile_common, profvals
;  prof_start1 = systime(1)

  if n_elements(neval) EQ 0 then neval = 0L
  epmach = r1mach4
  uflow = r1mach1

  zero = a*b*0.
  centr = 0.5e+00*(a+b)
  hlgth = 0.5e+00*(b-a)
  dhlgth = abs(hlgth)

  x = centr+hlgth*xgk
  if n_elements(bound) GT 0 then begin
      ;; Transformation for infinite integrals
      u = temporary(x)
      x = bound(0) + dirsign*(1/u - 1)
  endif

;  prof_start2 = systime(1)
  ;; Call with or without PRIVATE and _EXTRA keywords
  if n_elements(priv) GT 0 then begin
      if n_elements(fa) GT 0 then fv = call_function(f, x, priv, _EXTRA=fa) $
      else                        fv = call_function(f, x, priv)
  endif else begin
      if n_elements(fa) GT 0 then fv = call_function(f, x, _EXTRA=fa) $
      else                        fv = call_function(f, x)
  endelse
;  profvals.functime = profvals.functime + (systime(1) - prof_start2)
  neval = neval + n_elements(x)
  if n_elements(fv) NE n_elements(x) then $
    message, 'ERROR: Integrand function '+strupcase(f)+$
    ' must return a vector of values'
  
  if n_elements(bound) GT 0 then begin
      ;; Complete the transformation for infinite integrals
      fv = fv / u^2
  endif

  resk = total(wgk*fv)
  resabs = total(wgk*abs(fv))
  resg = total(wg*fv(ig))

  resasc = total(wgk*abs(fv-resk*0.5e+00))

  result = resk*hlgth
  resabs = resabs*dhlgth
  resasc = resasc*dhlgth
  abserr = abs((resk-resg)*hlgth)
  if (resasc NE 0.0e+00 AND abserr NE 0.0e+00) then $
    abserr = resasc*min([0.1e+01,(0.2e+03*abserr/resasc)^1.5])
  if (resabs GT uflow/(0.5e+02*epmach)) then $
    abserr = max([(epmach*0.5e+02)*resabs,abserr])
;  profvals.qktime = profvals.qktime + (systime(1) - prof_start1)
  return
end



pro qpint1d_qelg, n, epstab, result, abserr, res3la, nres

  ;; Derived from QUADPACK QELG
  ;;***PURPOSE  The routine determines the limit of a given sequence of
  ;;            approximations, by means of the Epsilon algorithm of
  ;;            P. Wynn. An estimate of the absolute error is also given.
  ;;            The condensed Epsilon table is computed. Only those
  ;;            elements needed for the computation of the next diagonal
  ;;            are preserved.
  ;;***LIBRARY   SLATEC

  common qpint1d_machar, machvals
  r1mach1 = machvals.minnum
  r1mach2 = machvals.maxnum
  r1mach4 = machvals.machep
  common qpint1d_profile_common, profvals
;  prof_start = systime(1)

  epmach = r1mach4
  oflow = r1mach2
  nres = nres+1
  abserr = oflow
  result = epstab(n-1)
  if (n LT 3) then goto, LAB100
  limexp = 50
  epstab(n+2-1) = epstab(n-1)
  newelm = (n-1)/2
  epstab(n-1) = oflow
  num = n
  k1 = n-1  ;; OK IDL 
  for i = 1, newelm do begin
      res = epstab(k1+2) ;; OK IDL
      e0 = epstab(k1-2)  ;; OK IDL
      e1 = epstab(k1-1)  ;; OK IDL
      e2 = res
      e1abs = abs(e1)
      delta2 = e2-e1
      err2 = abs(delta2)
      tol2 = max([abs(e2),e1abs])*epmach
      delta3 = e1-e0
      err3 = abs(delta3)
      tol3 = max([e1abs,abs(e0)])*epmach

      if NOT (err2 GT tol2 OR err3 GT tol3) then begin
;           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
;           ACCURACY, CONVERGENCE IS ASSUMED.
          result = res
          abserr = err2+err3
; ***JUMP OUT OF DO-LOOP
          goto, LAB100
      endif
      e3 = epstab(k1)
      epstab(k1) = e1
      delta1 = e1-e3
      err1 = abs(delta1)
      tol1 = max([e1abs,abs(e3)])*epmach

;           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
;           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
      if NOT (err1 LE tol1 OR err2 LE tol2 OR err3 LE tol3) then begin

          ;; Avoid underflow errors
          if abs(delta1) GE 0.5*abs(oflow) then odelta1 = delta1*0 $
          else                                  odelta1 = 1/delta1
          ss = odelta1+0.1e+01/delta2-0.1e+01/delta3
          epsinf = abs(ss*e1)
          
;           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
;           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
;           OF N.
          if (epsinf GT 0.1e-03) then goto, LAB30
      endif

      
      n = i+i-1
; ***JUMP OUT OF DO-LOOP
      goto, LAB50
      
;           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
;           THE VALUE OF RESULT.
      LAB30:
      res = e1+0.1e+01/ss
      epstab(k1) = res
      k1 = k1-2
      error = err2+abs(res-e2)+err3
      if NOT (error GT abserr) then begin
          abserr = error
          result = res
      endif
  endfor
  
;           SHIFT THE TABLE.
  LAB50:
  if (n EQ limexp) then n = 2*(limexp/2)-1
  ib = 1
  if ((num/2)*2 EQ num) then ib = 2
  ie = newelm+1
  ii = ib-1+lindgen(ie)*2
  epstab(ii) = epstab(ii+2) ;; OK IDL
  if NOT (num EQ n) then begin
      indx = num-n+1
      ii = lindgen(n)
      epstab(ii) = epstab(indx-1+ii) ;; OK IDL
  endif
  
  if NOT (nres GE 4) then begin
      res3la(nres-1) = result
      abserr = oflow
  endif else begin
      
;           COMPUTE ERROR ESTIMATE - OK IDL
      abserr = (abs(result-res3la(2))+abs(result-res3la(1)) $
                +abs(result-res3la(0)))
      res3la(0) = res3la(1)
      res3la(1) = res3la(2)
      res3la(2) = result
  endelse
  
  LAB100:
  abserr = max([abserr,0.5e+01*epmach*abs(result)])
;  profvals.acctime = profvals.acctime + (systime(1) - prof_start)
  return
end
  
pro qpint1d_qpsrt, limit, last, maxerr, ermax, elist, iord, nrmax
;;                              OOOOOO  OOOOO         OOOO  OOOOO

  ;; Derived from QUADPACK QPSRT
  ;;***PURPOSE  Subsidiary to QAGE, QAGIE, QAGPE, QAGSE, QAWCE, QAWOE and
  ;;            QAWSE
  ;;***LIBRARY   SLATEC

  ;; ELIST - list of unsorted errors
  ;; IORD - (O) list of indices which sort ELIST in descending order
  ;; LAST - index of last element (new elt) in ELIST
  ;; LIMIT - maximum size of ELIST

  ;; NRMAX - (O) position in IORD of the maximum error
  ;; MAXERR - (O) position in ELIST of the maximum error
  ;; ERMAX - (O) amount of maximum error
  common qpint1d_profile_common, profvals
;  prof_start = systime(1)

  if (last LE 2) then begin
      iord(1-1) = 1
      iord(2-1) = 2
      goto, LAB90
  endif

;           THIS PART OF THE ROUTINE IS ONLY EXECUTED
;           IF, DUE TO A DIFFICULT INTEGRAND, SUBDIVISION
;           INCREASED THE ERROR ESTIMATE. IN THE NORMAL CASE
;           THE INSERT PROCEDURE SHOULD START AFTER THE
;           NRMAX-TH LARGEST ERROR ESTIMATE.

  errmax = elist(maxerr-1)
  if (nrmax NE 1) then begin
      ido = nrmax-1
      
      for i = 1, ido do begin
          isucc = iord(nrmax-1-1)
; ***JUMP OUT OF DO-LOOP
          if (errmax LE elist(isucc-1)) then goto, LAB30
          iord(nrmax-1) = isucc
          nrmax = nrmax-1
      endfor
  endif

  LAB30:
  jupbn = last
  if (last GT (limit/2+2)) then jupbn = limit+3-last
  errmin = elist(last-1)

;           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
;           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
  jbnd = jupbn-1
  ibeg = nrmax+1
  if (ibeg LE jbnd) then begin
      for i = ibeg, jbnd do begin
          isucc = iord(i-1)
; ***JUMP OUT OF DO-LOOP
          if (errmax GE elist(isucc-1)) then goto, LAB60
          iord(i-1-1) = isucc
      endfor
  endif
  iord(jbnd-1) = maxerr
  iord(jupbn-1) = last
  goto, LAB90
  
;           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
  LAB60:
  iord(i-1-1) = maxerr
  k = jbnd
  for j=i, jbnd do begin
      isucc = iord(k-1)
; ***JUMP OUT OF DO-LOOP
      if (errmin LT elist(isucc-1)) then goto, LAB80
      iord(k+1-1) = isucc
      k = k-1
  endfor
  iord(i-1) = last
  goto, LAB90
  LAB80:
  iord(k+1-1) = last

;           SET MAXERR AND ERMAX.
  LAB90:
  maxerr = iord(nrmax-1)
  ermax = elist(maxerr-1)
;  profvals.srttime = profvals.srttime + (systime(1) - prof_start)
  return
end

pro qpint1d_gkweights, wg, wgk, xgk, ig, nptsreq, nptsact, prec


  ;; The user has requested NPTSREQ points, however that number may
  ;; not be available.  Filter the number and record the *actual*
  ;; number of points used, NPTSACT


;           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
;           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
;           CORRESPONDING WEIGHTS ARE GIVEN.
;
;           XGK    - ABSCISSAE OF THE 21-POINT KRONROD RULE
;                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 10-POINT
;                    GAUSS RULE
;                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
;                    ADDED TO THE 10-POINT GAUSS RULE
;
;           WGK    - WEIGHTS OF THE 21-POINT KRONROD RULE
;
;           WG     - WEIGHTS OF THE 10-POINT GAUSS RULE

      ;; gauss quadrature weights and kronron quadrature abscissae and
      ;; weights as evaluated with 80 decimal digit arithmetic by
      ;; l. w. fullerton, bell labs, nov. 1981.
    
  if nptsreq LE 15 then begin
      nptsact = 15
      wg = [ $
             0.129484966168869693270611432679082d0, $
             0.279705391489276667901467771423780d0, $
             0.381830050505118944950369775488975d0, $
             0.417959183673469387755102040816327d0  ]

      xgk = [ $
              0.991455371120812639206854697526329d0, $
              0.949107912342758524526189684047851d0, $
              0.864864423359769072789712788640926d0, $
              0.741531185599394439863864773280788d0, $
              0.586087235467691130294144838258730d0, $
              0.405845151377397166906606412076961d0, $
              0.207784955007898467600689403773245d0, $
              0.000000000000000000000000000000000d0  ]

      wgk = [ $
              0.022935322010529224963732008058970d0, $
              0.063092092629978553290700663189204d0, $
              0.104790010322250183839876322541518d0, $
              0.140653259715525918745189590510238d0, $
              0.169004726639267902826583426598550d0, $
              0.190350578064785409913256402421014d0, $
              0.204432940075298892414161999234649d0, $
              0.209482141084727828012999174891714d0  ]
              
  endif else if nptsreq LE 21 then begin
      nptsact = 21
      wg = [  0.066671344308688137593568809893332d0, $
              0.149451349150580593145776339657697d0, $
              0.219086362515982043995534934228163d0, $
              0.269266719309996355091226921569469d0, $
              0.295524224714752870173892994651338d0  ]

      xgk = [ 0.995657163025808080735527280689003d0, $
              0.973906528517171720077964012084452d0, $
              0.930157491355708226001207180059508d0, $
              0.865063366688984510732096688423493d0, $
              0.780817726586416897063717578345042d0, $
              0.679409568299024406234327365114874d0, $
              0.562757134668604683339000099272694d0, $
              0.433395394129247190799265943165784d0, $
              0.294392862701460198131126603103866d0, $
              0.148874338981631210884826001129720d0, $
              0.000000000000000000000000000000000d0  ]

      wgk = [ 0.011694638867371874278064396062192d0, $
              0.032558162307964727478818972459390d0, $
              0.054755896574351996031381300244580d0, $
              0.075039674810919952767043140916190d0, $
              0.093125454583697605535065465083366d0, $
              0.109387158802297641899210590325805d0, $
              0.123491976262065851077958109831074d0, $
              0.134709217311473325928054001771707d0, $
              0.142775938577060080797094273138717d0, $
              0.147739104901338491374841515972068d0, $
              0.149445554002916905664936468389821d0  ]

  endif else if nptsreq LE 31 then begin
      nptsact = 31
      wg = [ $
             0.03075324199611726834628393577204d0, $
             0.07036604748810812479267416450667d0, $
             0.10715922046717193501869546685869d0, $
             0.13957067792615431447804794511028d0, $
             0.16626920581699393353200860481209d0, $
             0.18616100001556221106800561866423d0, $
             0.19843148532711157646118326443839d0, $
             0.20257824192556127280620199967519d0  ]

      xgk = [ $
              0.99800229869339706025172840152271d0, $
              0.98799251802048542849565718586613d0, $
              0.96773907567913913427347978784337d0, $
              0.93727339240070590437758947710209d0, $
              0.89726453234408190082509656454496d0, $
              0.84820658341042721620648320774217d0, $
              0.79041850144246593297649294817947d0, $
              0.72441773136017004746186054613938d0, $
              0.65099674129741697053735895313275d0, $
              0.57097217260853884757226737253911d0, $
              0.48508186364023968063655740232351d0, $
              0.39415134707756336987207370981045d0, $
              0.29918000715316881216780024266389d0, $
              0.20119409399743452230628303394596d0, $
              0.10114206691871749907074231447392d0, $
              0.00000000000000000000000000000000d0  ]

      wgk = [ $
              0.00537747987292334897792051430128d0, $
              0.01500794732931612258374763075807d0, $
              0.02546084732671532016874001019653d0, $
              0.03534636079137584622037948478360d0, $
              0.04458975132476487668227299373280d0, $
              0.05348152469092808725343147239430d0, $
              0.06200956780067064025139230960803d0, $
              0.06985412131872825879520077099147d0, $
              0.07684968075772037884432777482659d0, $
              0.08308050282313302108289247286104d0, $
              0.08856444305621177067275443693774d0, $
              0.09312659817082532125486872747346d0, $
              0.09664272698362367855179907627589d0, $
              0.09917359872179195932393173484603d0, $
              0.10076984552387559504946662617570d0, $
              0.10133000701479154907374792767493d0  ]
  endif else if nptsreq LE 41 then begin
      nptsact = 41
      wg = [ $
             0.017614007139152118311861962351853d0, $
             0.040601429800386941331039952274932d0, $
             0.062672048334109063569506535187042d0, $
             0.083276741576704748724758143222046d0, $
             0.101930119817240435036750135480350d0, $
             0.118194531961518417312377377711382d0, $
             0.131688638449176626898494499748163d0, $
             0.142096109318382051329298325067165d0, $
             0.149172986472603746787828737001969d0, $
             0.152753387130725850698084331955098d0  ]
      xgk = [ $
              0.998859031588277663838315576545863d0, $
              0.993128599185094924786122388471320d0, $
              0.981507877450250259193342994720217d0, $
              0.963971927277913791267666131197277d0, $
              0.940822633831754753519982722212443d0, $
              0.912234428251325905867752441203298d0, $
              0.878276811252281976077442995113078d0, $
              0.839116971822218823394529061701521d0, $
              0.795041428837551198350638833272788d0, $
              0.746331906460150792614305070355642d0, $
              0.693237656334751384805490711845932d0, $
              0.636053680726515025452836696226286d0, $
              0.575140446819710315342946036586425d0, $
              0.510867001950827098004364050955251d0, $
              0.443593175238725103199992213492640d0, $
              0.373706088715419560672548177024927d0, $
              0.301627868114913004320555356858592d0, $
              0.227785851141645078080496195368575d0, $
              0.152605465240922675505220241022678d0, $
              0.076526521133497333754640409398838d0, $
              0.000000000000000000000000000000000d0  ]
      wgk = [ $
              0.003073583718520531501218293246031d0, $
              0.008600269855642942198661787950102d0, $
              0.014626169256971252983787960308868d0, $
              0.020388373461266523598010231432755d0, $
              0.025882133604951158834505067096153d0, $
              0.031287306777032798958543119323801d0, $
              0.036600169758200798030557240707211d0, $
              0.041668873327973686263788305936895d0, $
              0.046434821867497674720231880926108d0, $
              0.050944573923728691932707670050345d0, $
              0.055195105348285994744832372419777d0, $
              0.059111400880639572374967220648594d0, $
              0.062653237554781168025870122174255d0, $
              0.065834597133618422111563556969398d0, $
              0.068648672928521619345623411885368d0, $
              0.071054423553444068305790361723210d0, $
              0.073030690332786667495189417658913d0, $
              0.074582875400499188986581418362488d0, $
              0.075704497684556674659542775376617d0, $
              0.076377867672080736705502835038061d0, $
              0.076600711917999656445049901530102d0  ]
  endif else if nptsreq LE 51 then begin
      nptsact = 51

      wg = [ $
             0.011393798501026287947902964113235d0, $
             0.026354986615032137261901815295299d0, $
             0.040939156701306312655623487711646d0, $
             0.054904695975835191925936891540473d0, $
             0.068038333812356917207187185656708d0, $
             0.080140700335001018013234959669111d0, $
             0.091028261982963649811497220702892d0, $
             0.100535949067050644202206890392686d0, $
             0.108519624474263653116093957050117d0, $
             0.114858259145711648339325545869556d0, $
             0.119455763535784772228178126512901d0, $
             0.122242442990310041688959518945852d0, $
             0.123176053726715451203902873079050d0  ]

      xgk = [ $
              0.999262104992609834193457486540341d0, $
              0.995556969790498097908784946893902d0, $
              0.988035794534077247637331014577406d0, $
              0.976663921459517511498315386479594d0, $
              0.961614986425842512418130033660167d0, $
              0.942974571228974339414011169658471d0, $
              0.920747115281701561746346084546331d0, $
              0.894991997878275368851042006782805d0, $
              0.865847065293275595448996969588340d0, $
              0.833442628760834001421021108693570d0, $
              0.797873797998500059410410904994307d0, $
              0.759259263037357630577282865204361d0, $
              0.717766406813084388186654079773298d0, $
              0.673566368473468364485120633247622d0, $
              0.626810099010317412788122681624518d0, $
              0.577662930241222967723689841612654d0, $
              0.526325284334719182599623778158010d0, $
              0.473002731445714960522182115009192d0, $
              0.417885382193037748851814394594572d0, $
              0.361172305809387837735821730127641d0, $
              0.303089538931107830167478909980339d0, $
              0.243866883720988432045190362797452d0, $
              0.183718939421048892015969888759528d0, $
              0.122864692610710396387359818808037d0, $
              0.061544483005685078886546392366797d0, $
              0.000000000000000000000000000000000d0  ]

      wgk = [ $
              0.001987383892330315926507851882843d0, $
              0.005561932135356713758040236901066d0, $
              0.009473973386174151607207710523655d0, $
              0.013236229195571674813656405846976d0, $
              0.016847817709128298231516667536336d0, $
              0.020435371145882835456568292235939d0, $
              0.024009945606953216220092489164881d0, $
              0.027475317587851737802948455517811d0, $
              0.030792300167387488891109020215229d0, $
              0.034002130274329337836748795229551d0, $
              0.037116271483415543560330625367620d0, $
              0.040083825504032382074839284467076d0, $
              0.042872845020170049476895792439495d0, $
              0.045502913049921788909870584752660d0, $
              0.047982537138836713906392255756915d0, $
              0.050277679080715671963325259433440d0, $
              0.052362885806407475864366712137873d0, $
              0.054251129888545490144543370459876d0, $
              0.055950811220412317308240686382747d0, $
              0.057437116361567832853582693939506d0, $
              0.058689680022394207961974175856788d0, $
              0.059720340324174059979099291932562d0, $
              0.060539455376045862945360267517565d0, $
              0.061128509717053048305859030416293d0, $
              0.061471189871425316661544131965264d0, $
              0.061580818067832935078759824240066d0  ]
;       note: wgk (26) was calculated from the values of wgk(1..25)

  endif else begin
      nptsact = 61

      wg = [ $
             0.007968192496166605615465883474674d0, $
             0.018466468311090959142302131912047d0, $
             0.028784707883323369349719179611292d0, $
             0.038799192569627049596801936446348d0, $
             0.048402672830594052902938140422808d0, $
             0.057493156217619066481721689402056d0, $
             0.065974229882180495128128515115962d0, $
             0.073755974737705206268243850022191d0, $
             0.080755895229420215354694938460530d0, $
             0.086899787201082979802387530715126d0, $
             0.092122522237786128717632707087619d0, $
             0.096368737174644259639468626351810d0, $
             0.099593420586795267062780282103569d0, $
             0.101762389748405504596428952168554d0, $
             0.102852652893558840341285636705415d0  ]

      xgk = [ $
              0.999484410050490637571325895705811d0, $
              0.996893484074649540271630050918695d0, $
              0.991630996870404594858628366109486d0, $
              0.983668123279747209970032581605663d0, $
              0.973116322501126268374693868423707d0, $
              0.960021864968307512216871025581798d0, $
              0.944374444748559979415831324037439d0, $
              0.926200047429274325879324277080474d0, $
              0.905573307699907798546522558925958d0, $
              0.882560535792052681543116462530226d0, $
              0.857205233546061098958658510658944d0, $
              0.829565762382768397442898119732502d0, $
              0.799727835821839083013668942322683d0, $
              0.767777432104826194917977340974503d0, $
              0.733790062453226804726171131369528d0, $
              0.697850494793315796932292388026640d0, $
              0.660061064126626961370053668149271d0, $
              0.620526182989242861140477556431189d0, $
              0.579345235826361691756024932172540d0, $
              0.536624148142019899264169793311073d0, $
              0.492480467861778574993693061207709d0, $
              0.447033769538089176780609900322854d0, $
              0.400401254830394392535476211542661d0, $
              0.352704725530878113471037207089374d0, $
              0.304073202273625077372677107199257d0, $
              0.254636926167889846439805129817805d0, $
              0.204525116682309891438957671002025d0, $
              0.153869913608583546963794672743256d0, $
              0.102806937966737030147096751318001d0, $
              0.051471842555317695833025213166723d0, $
              0.000000000000000000000000000000000d0  ]

      wgk = [ $
              0.001389013698677007624551591226760d0, $
              0.003890461127099884051267201844516d0, $
              0.006630703915931292173319826369750d0, $
              0.009273279659517763428441146892024d0, $
              0.011823015253496341742232898853251d0, $
              0.014369729507045804812451432443580d0, $
              0.016920889189053272627572289420322d0, $
              0.019414141193942381173408951050128d0, $
              0.021828035821609192297167485738339d0, $
              0.024191162078080601365686370725232d0, $
              0.026509954882333101610601709335075d0, $
              0.028754048765041292843978785354334d0, $
              0.030907257562387762472884252943092d0, $
              0.032981447057483726031814191016854d0, $
              0.034979338028060024137499670731468d0, $
              0.036882364651821229223911065617136d0, $
              0.038678945624727592950348651532281d0, $
              0.040374538951535959111995279752468d0, $
              0.041969810215164246147147541285970d0, $
              0.043452539701356069316831728117073d0, $
              0.044814800133162663192355551616723d0, $
              0.046059238271006988116271735559374d0, $
              0.047185546569299153945261478181099d0, $
              0.048185861757087129140779492298305d0, $
              0.049055434555029778887528165367238d0, $
              0.049795683427074206357811569379942d0, $
              0.050405921402782346840893085653585d0, $
              0.050881795898749606492297473049805d0, $
              0.051221547849258772170656282604944d0, $
              0.051426128537459025933862879215781d0, $
              0.051494729429451567558340433647099d0  ]
  endelse

  nhalf  = (nptsact-1)/2
  nquart = (nptsact-1)/4
  wgk = [  wgk(0:nhalf-1), rotate(wgk,2) ]
  xgk = [ -xgk(0:nhalf-1), rotate(xgk,2) ]
  wg  = [wg(0:nquart-1), rotate(wg,2)]
  ig  = lindgen(nhalf)*2 + 1

  ;; Convert to float from double if requested
  if prec EQ 0 then begin
      wgk = float(wgk)
      xgk = float(xgk)
      wg  = float(wg)
  endif
end

;; Evaluate a user-supplied expression
function qpint1d_eval, x, p, expression=expr
  y = 0
  cmd = 'Y = '+expr
  dummy = execute(cmd)
  if dummy NE 1 then begin
      err_string = ''+!error_state.msg
      message, 'ERROR: failed to execute expression "'+expr+'" ('+err_string+')'
  endif
  return, y
end

;; The outer routine which does most of the preparation and special
;; cases.
function qpint1d, f0, a0, b0, private, npoints=npoints, expression=expr, $
                  epsabs=epsabs, epsrel=epsrel, error=abserr, nfev=neval, $
                  status=ier, functargs=fa, limit=limit, nsubintervals=nsub, $
                  sym_axis=symaxis, symmetric=symmetric, breakpoints=bpoints0,$
                  alist=alist, blist=blist, rlist=rlist, $
                  elist=elist, iord=iord, nocatch=nocatch

  ;; Derived from QUADPACK QAGS
  ;;***PURPOSE  The routine calculates an approximation result to a given
  ;;            Definite integral  I = Integral of F over (A,B),
  ;;            Hopefully satisfying following claim for accuracy
  ;;           ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
  ;;***LIBRARY   SLATEC (QUADPACK)

  if keyword_set(query) then return, 1
  common qpint1d_profile_common, profvals
  qpint1d_profile

  if n_params() EQ 0 then begin
      USAGE:
      message, 'USAGE:', /info
      message, ' G = QPINT1D(FUNCNAME, A, B, $', /info
      message, '             [EPSABS=, EPSREL=, ERROR=, STATUS=])', /info
      message, '      (or)', /info
      message, ' G = QPINT1D(EXPR, A, B, /EXPRESSION, $', /info
      message, '             [EPSABS=, EPSREL=, ERROR=, STATUS=])', /info
      return, !values.d_nan
  endif
  if n_elements(f0) EQ 0 OR n_elements(a0) EQ 0 OR n_elements(b0) EQ 0 then $
    goto, USAGE

  sz = size(f0)
  if sz(sz(0)+1) NE 7 OR n_elements(f0) NE 1 then begin
      message, 'ERROR: FUNCT must be a scalar string', /info
      return, !values.d_nan
  endif

  ;; Handle error conditions gracefully
  if NOT keyword_set(nocatch) then begin
      catch, catcherror
      if catcherror NE 0 then begin
          catch, /cancel
          message, 'Error detected while integrating "'+f0(0)+'"', /info
          message, !err_string, /info
          errmsg = 0
          if NOT keyword_set(expr) then begin
              f1 = byte(strupcase(strtrim(f0(0),2)))
              ca = (byte('A'))(0) 
              cz = (byte('Z'))(0) 
              c0 = (byte('0'))(0) 
              c9 = (byte('9'))(0) 
              c_ = (byte('_'))(0) 
              wh = where((f1 GE ca AND f1 LE cz) EQ 0 AND f1 NE c_ $
                         AND (f1 GE c0 AND f1 LE c9) EQ 0, ct)
              if ct GT 0 OR (f1(0) GE c0 AND f1(0) LE c9) then begin
                  message, ('FUNCT appears to be an expression.  Did you '+$
                            'intend to pass the /EXPRESSION keyword?'), /info
                  errmsg = 1
              endif
          endif
          if errmsg EQ 0 then $
            message, ('Please verify that function works and conforms to '+$
                      'the documentation'), /info
          ier = -1L
          return, !values.d_nan
      endif
  endif
  
;  prof_start = systime(1)
  neval = 0L
  sza = size(a0)
  szb = size(b0)
  ;; Determine the data precision, i.e. single or double precision
  if sza(sza(0)+1) EQ 5 OR szb(szb(0)+1) EQ 5 then begin
      isdouble = 1
      zero = 0D
      reldef = 1e-6
      a = double(a0(0))
      b = double(b0(0))
  endif else begin
      isdouble = 0
      zero = 0.
      reldef = 1e-4
      a = float(a0(0))
      b = float(b0(0))
  endelse
  qpint1d_setmachar, double=isdouble

  ;; Default values
  if n_elements(limit) EQ 0 then limit = 100L
  if n_elements(epsrel) EQ 0 then epsrel = zero + reldef
  if n_elements(epsabs) EQ 0 then epsabs = zero

  ;; Establish the initial return values return
  abserr = zero
  result = zero
  if a EQ b then begin
      ier = 0L
      return, result
  endif

  ;; Now prepare for potentially one or two integrals, depending on
  ;; whether the integral is fully infinite, or only partially
  infa = finite(a) EQ 0
  infb = finite(b) EQ 0

  if infa AND infb then begin
      ;; Fully infinite: break it into two partial integrals
      if n_elements(symaxis0) EQ 0 then symaxis = zero $
      else                              symaxis = zero + symaxis0(0)

      a1 = symaxis & b1 = b
      a2 = symaxis & b2 = a
      ef = 2.  ;; Each integral contributes half of error budget
  endif else if infa then begin
      ;; Partially infinite, [a,inf]: do only second partial integral
      a1 = zero    & b1 = zero
      a2 = b       & b2 = a
      ef = 1.
  endif else begin
      ;; Partially infinite, [inf,b]: do only first partial integral
      a1 = a       & b1 = b
      a2 = zero    & b2 = zero
      ef = 1.
  endelse

  ;; Prepare for EXPRESSION if requested
  if keyword_set(expr) then begin
      f = 'QPINT1D_EVAL'
      fa = {expression: strtrim(f0(0),2)}
  endif else begin
      f = strtrim(f0(0),2)
  endelse

  ;; Call first partial integral if requested
  ier = 0L
  neval = 0L
  nsub = 0L
  if a1 NE b1 then begin
      qpint1d_qagse, f, a1, b1, result, abserr, private, functargs=fa, $
        epsabs=epsabs/ef, epsrel=epsrel/ef, breakpoints=bpoints0, $
        status=ier, limit=limit, neval=neval, nsubintervals=nsub, $
        alist=alist, blist=blist, rlist=rlist, elist=elist, iord=iord, $
        npoints=npoints
  endif

  ;; Return if an error condition was detected
  if ier EQ 6 then ier = -1L
  if ier EQ -1 OR ier EQ 3 then return, result

  ;; Now call the second partial integral if requested
  if a2 NE b2 then begin
      qpint1d_qagse, f, a2, b2, result2, abserr2, private, functargs=fa, $
        epsabs=epsabs/ef, epsrel=epsrel/ef, breakpoints=bpoints0, $
        status=ier, limit=limit, neval=neval2, nsubintervals=nsub2, $
        alist=alist, blist=blist, rlist=rlist, elist=elist, iord=iord, $
        npoints=npoints

      ;; Merge the two results together
      result = result - result2
      abserr = abserr + abserr2
      neval  = neval  + neval2
      nsub   = nsub   + nsub2
  endif

;  profvals.tottime = profvals.tottime + (systime(1) - prof_start)

  ;; Convert from QUADPACK to QPINT1D errors
  if ier EQ 6 then ier = -1L
  return, result
end