Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Namespaces | Functions
Garfield::Numerics Namespace Reference

Collection of numerical routines. More...

Namespaces

namespace  CERNLIB
 Linear algebra routines from CERNLIB.
 
namespace  QUADPACK
 

Functions

double Legendre (const unsigned int n, const double x)
 Legendre polynomials.
 
double BesselI0S (const double xx)
 
double BesselI1S (const double xx)
 
double BesselK0S (const double xx)
 
double BesselK0L (const double xx)
 
double BesselK1S (const double xx)
 
double BesselK1L (const double xx)
 
double Divdif (const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
 
bool Boxin2 (const std::vector< std::vector< double > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const int nx, const int ny, const double xx, const double yy, double &f, const int iOrder)
 
bool Boxin3 (const std::vector< std::vector< std::vector< double > > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const std::vector< double > &zAxis, const int nx, const int ny, const int nz, const double xx, const double yy, const double zz, double &f, const int iOrder)
 
bool LeastSquaresFit (std::function< double(double, const std::vector< double > &)> f, std::vector< double > &par, std::vector< double > &epar, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &ey, const unsigned int nMaxIter, const double diff, double &chi2, const double eps, const bool debug, const bool verbose)
 Least-squares minimisation.
 

Detailed Description

Collection of numerical routines.

Function Documentation

◆ BesselI0S()

double Garfield::Numerics::BesselI0S ( const double  xx)
inline

Modified Bessel functions. Series expansions from Abramowitz and Stegun.

Definition at line 135 of file Numerics.hh.

135 {
136 const double y = xx / 3.75;
137 const double y2 = y * y;
138 return 1. + 3.5156229 * y2 + 3.0899424 * y2 * y2 +
139 1.2067492 * pow(y2, 3) + 0.2659732 * pow(y2, 4) +
140 0.0360768 * pow(y2, 5) + 0.0045813 * pow(y2, 6);
141}

Referenced by BesselK0S().

◆ BesselI1S()

double Garfield::Numerics::BesselI1S ( const double  xx)
inline

Definition at line 143 of file Numerics.hh.

143 {
144 const double y = xx / 3.75;
145 const double y2 = y * y;
146 return xx *
147 (0.5 + 0.87890594 * y2 + 0.51498869 * y2 * y2 +
148 0.15084934 * pow(y2, 3) + 0.02658733 * pow(y2, 4) +
149 0.00301532 * pow(y2, 5) + 0.00032411 * pow(y2, 6));
150}

Referenced by BesselK1S().

◆ BesselK0L()

double Garfield::Numerics::BesselK0L ( const double  xx)
inline

Definition at line 161 of file Numerics.hh.

161 {
162 const double y = 2. / xx;
163 return (exp(-xx) / sqrt(xx)) *
164 (1.25331414 - 0.07832358 * y + 0.02189568 * y * y -
165 0.01062446 * pow(y, 3) + 0.00587872 * pow(y, 4) -
166 0.00251540 * pow(y, 5) + 0.00053208 * pow(y, 6));
167}

◆ BesselK0S()

double Garfield::Numerics::BesselK0S ( const double  xx)
inline

Definition at line 152 of file Numerics.hh.

152 {
153 const double y = xx / 2.;
154 const double y2 = y * y;
155 return -log(y) * BesselI0S(xx) - 0.57721566 +
156 0.42278420 * y2 + 0.23069756 * y2 * y2 +
157 0.03488590 * pow(y2, 3) + 0.00262698 * pow(y2, 4) +
158 0.00010750 * pow(y2, 5) + 0.00000740 * pow(y2, 6);
159}
Garfield::Numerics::BesselI0S
double BesselI0S(const double xx)
Definition: Numerics.hh:135

◆ BesselK1L()

double Garfield::Numerics::BesselK1L ( const double  xx)
inline

Definition at line 179 of file Numerics.hh.

179 {
180 const double y = 2. / xx;
181 return (exp(-xx) / sqrt(xx)) *
182 (1.25331414 + 0.23498619 * y - 0.03655620 * y * y +
183 0.01504268 * pow(y, 3) - 0.00780353 * pow(y, 4) +
184 0.00325614 * pow(y, 5) - 0.00068245 * pow(y, 6));
185}

◆ BesselK1S()

double Garfield::Numerics::BesselK1S ( const double  xx)
inline

Definition at line 169 of file Numerics.hh.

169 {
170 const double y = xx / 2.;
171 const double y2 = y * y;
172 return log(y) * BesselI1S(xx) +
173 (1. / xx) *
174 (1. + 0.15443144 * y2 - 0.67278579 * y2 * y2 -
175 0.18156897 * pow(y2, 3) - 0.01919402 * pow(y2, 4) -
176 0.00110404 * pow(y2, 5) - 0.00004686 * pow(y2, 6));
177}
Garfield::Numerics::BesselI1S
double BesselI1S(const double xx)
Definition: Numerics.hh:143

◆ Boxin2()

bool Garfield::Numerics::Boxin2 ( const std::vector< std::vector< double > > &  value,
const std::vector< double > &  xAxis,
const std::vector< double > &  yAxis,
const int  nx,
const int  ny,
const double  xx,
const double  yy,
double &  f,
const int  iOrder 
)

Interpolation of order 1 and 2 in an irregular rectangular two-dimensional grid.

Definition at line 1305 of file Numerics.cc.

1308 {
1309 //-----------------------------------------------------------------------
1310 // BOXIN2 - Interpolation of order 1 and 2 in an irregular rectangular
1311 // 2-dimensional grid.
1312 //-----------------------------------------------------------------------
1313 int iX0 = 0, iX1 = 0;
1314 int iY0 = 0, iY1 = 0;
1315 std::array<double, 3> fX;
1316 std::array<double, 3> fY;
1317 f = 0.;
1318 // Ensure we are in the grid.
1319 if ((xAxis[nx - 1] - x) * (x - xAxis[0]) < 0 ||
1320 (yAxis[ny - 1] - y) * (y - yAxis[0]) < 0) {
1321 std::cerr << "Boxin2: Point not in the grid; no interpolation.\n";
1322 return false;
1323 }
1324 // Make sure we have enough points.
1325 if (iOrder < 0 || iOrder > 2) {
1326 std::cerr << "Boxin2: Incorrect order; no interpolation.\n";
1327 return false;
1328 } else if (nx < 1 || ny < 1) {
1329 std::cerr << "Boxin2: Incorrect number of points; no interpolation.\n";
1330 return false;
1331 }
1332 if (iOrder == 0 || nx <= 1) {
1333 // Zeroth order interpolation in x.
1334 // Find the nearest node.
1335 double dist = fabs(x - xAxis[0]);
1336 int iNode = 0;
1337 for (int i = 1; i < nx; i++) {
1338 if (fabs(x - xAxis[i]) < dist) {
1339 dist = fabs(x - xAxis[i]);
1340 iNode = i;
1341 }
1342 }
1343 // Set the summing range.
1344 iX0 = iNode;
1345 iX1 = iNode;
1346 // Establish the shape functions.
1347 fX = {1., 0., 0.};
1348 } else if (iOrder == 1 || nx <= 2) {
1349 // First order interpolation in x.
1350 // Find the grid segment containing this point.
1351 int iGrid = 0;
1352 for (int i = 1; i < nx; i++) {
1353 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1354 iGrid = i;
1355 }
1356 }
1357 const double x0 = xAxis[iGrid - 1];
1358 const double x1 = xAxis[iGrid];
1359 // Ensure there won't be divisions by zero.
1360 if (x1 == x0) {
1361 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1362 return false;
1363 }
1364 // Compute local coordinates.
1365 const double xL = (x - x0) / (x1 - x0);
1366 // Set the summing range.
1367 iX0 = iGrid - 1;
1368 iX1 = iGrid;
1369 // Set the shape functions.
1370 fX = {1. - xL, xL, 0.};
1371 } else if (iOrder == 2) {
1372 // Second order interpolation in x.
1373 // Find the nearest node and the grid segment.
1374 double dist = fabs(x - xAxis[0]);
1375 int iNode = 0;
1376 for (int i = 1; i < nx; ++i) {
1377 if (fabs(x - xAxis[i]) < dist) {
1378 dist = fabs(x - xAxis[i]);
1379 iNode = i;
1380 }
1381 }
1382 // Find the nearest fitting 2x2 matrix.
1383 int iGrid = std::max(1, std::min(nx - 2, iNode));
1384 const double x0 = xAxis[iGrid - 1];
1385 const double x1 = xAxis[iGrid];
1386 const double x2 = xAxis[iGrid + 1];
1387 // Ensure there won't be divisions by zero.
1388 if (x2 == x0) {
1389 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1390 return false;
1391 }
1392 // Compute the alpha and local coordinate for this grid segment.
1393 const double xAlpha = (x1 - x0) / (x2 - x0);
1394 const double xL = (x - x0) / (x2 - x0);
1395 // Ensure there won't be divisions by zero.
1396 if (xAlpha <= 0 || xAlpha >= 1) {
1397 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1398 return false;
1399 }
1400 // Set the summing range.
1401 iX0 = iGrid - 1;
1402 iX1 = iGrid + 1;
1403 // Set the shape functions.
1404 const double xL2 = xL * xL;
1405 fX[0] = xL2 / xAlpha - xL * (1. + xAlpha) / xAlpha + 1.;
1406 fX[1] = (xL2 - xL) / (xAlpha * xAlpha - xAlpha);
1407 fX[2] = (xL2 - xL * xAlpha) / (1. - xAlpha);
1408 }
1409 if (iOrder == 0 || ny <= 1) {
1410 // Zeroth order interpolation in y.
1411 // Find the nearest node.
1412 double dist = fabs(y - yAxis[0]);
1413 int iNode = 0;
1414 for (int i = 1; i < ny; i++) {
1415 if (fabs(y - yAxis[i]) < dist) {
1416 dist = fabs(y - yAxis[i]);
1417 iNode = i;
1418 }
1419 }
1420 // Set the summing range.
1421 iY0 = iNode;
1422 iY1 = iNode;
1423 // Establish the shape functions.
1424 fY = {1., 0., 0.};
1425 } else if (iOrder == 1 || ny <= 2) {
1426 // First order interpolation in y.
1427 // Find the grid segment containing this point.
1428 int iGrid = 0;
1429 for (int i = 1; i < ny; ++i) {
1430 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0) {
1431 iGrid = i;
1432 }
1433 }
1434 const double y0 = yAxis[iGrid - 1];
1435 const double y1 = yAxis[iGrid];
1436 // Ensure there won't be divisions by zero.
1437 if (y1 == y0) {
1438 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1439 return false;
1440 }
1441 // Compute local coordinates.
1442 const double yL = (y - y0) / (y1 - y0);
1443 // Set the summing range.
1444 iY0 = iGrid - 1;
1445 iY1 = iGrid;
1446 // Set the shape functions.
1447 fY = {1. - yL, yL, 0.};
1448 } else if (iOrder == 2) {
1449 // Second order interpolation in y.
1450 // Find the nearest node and the grid segment.
1451 double dist = fabs(y - yAxis[0]);
1452 int iNode = 0;
1453 for (int i = 1; i < ny; ++i) {
1454 if (fabs(y - yAxis[i]) < dist) {
1455 dist = fabs(y - yAxis[i]);
1456 iNode = i;
1457 }
1458 }
1459 // Find the nearest fitting 2x2 matrix.
1460 int iGrid = std::max(1, std::min(ny - 2, iNode));
1461 const double y0 = yAxis[iGrid - 1];
1462 const double y1 = yAxis[iGrid];
1463 const double y2 = yAxis[iGrid + 1];
1464 // Ensure there won't be divisions by zero.
1465 if (y2 == y0) {
1466 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1467 return false;
1468 }
1469 // Compute the alpha and local coordinate for this grid segment.
1470 const double yAlpha = (y1 - y0) / (y2 - y0);
1471 const double yL = (y - y0) / (y2 - y0);
1472 // Ensure there won't be divisions by zero.
1473 if (yAlpha <= 0 || yAlpha >= 1) {
1474 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1475 return false;
1476 }
1477 // Set the summing range.
1478 iY0 = iGrid - 1;
1479 iY1 = iGrid + 1;
1480 // Set the shape functions.
1481 const double yL2 = yL * yL;
1482 fY[0] = yL2 / yAlpha - yL * (1. + yAlpha) / yAlpha + 1.;
1483 fY[1] = (yL2 - yL) / (yAlpha * yAlpha - yAlpha);
1484 fY[2] = (yL2 - yL * yAlpha) / (1. - yAlpha);
1485 }
1486
1487 // Sum the shape functions.
1488 for (int i = iX0; i <= iX1; ++i) {
1489 for (int j = iY0; j <= iY1; ++j) {
1490 f += value[i][j] * fX[i - iX0] * fY[j - iY0];
1491 }
1492 }
1493 return true;
1494}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ Boxin3()

bool Garfield::Numerics::Boxin3 ( const std::vector< std::vector< std::vector< double > > > &  value,
const std::vector< double > &  xAxis,
const std::vector< double > &  yAxis,
const std::vector< double > &  zAxis,
const int  nx,
const int  ny,
const int  nz,
const double  xx,
const double  yy,
const double  zz,
double &  f,
const int  iOrder 
)

Interpolation of order 1 and 2 in an irregular rectangular three-dimensional grid.

Definition at line 1496 of file Numerics.cc.

1500 {
1501 //-----------------------------------------------------------------------
1502 // BOXIN3 - interpolation of order 1 and 2 in an irregular rectangular
1503 // 3-dimensional grid.
1504 //-----------------------------------------------------------------------
1505 int iX0 = 0, iX1 = 0;
1506 int iY0 = 0, iY1 = 0;
1507 int iZ0 = 0, iZ1 = 0;
1508 std::array<double, 4> fX;
1509 std::array<double, 4> fY;
1510 std::array<double, 4> fZ;
1511
1512 f = 0.;
1513 // Ensure we are in the grid.
1514 const double x = std::min(std::max(xx, std::min(xAxis[0], xAxis[nx - 1])),
1515 std::max(xAxis[0], xAxis[nx - 1]));
1516 const double y = std::min(std::max(yy, std::min(yAxis[0], yAxis[ny - 1])),
1517 std::max(yAxis[0], yAxis[ny - 1]));
1518 const double z = std::min(std::max(zz, std::min(zAxis[0], zAxis[nz - 1])),
1519 std::max(zAxis[0], zAxis[nz - 1]));
1520
1521 // Make sure we have enough points.
1522 if (iOrder < 0 || iOrder > 2) {
1523 std::cerr << "Boxin3: Incorrect order; no interpolation.\n";
1524 return false;
1525 } else if (nx < 1 || ny < 1 || nz < 1) {
1526 std::cerr << "Boxin3: Incorrect number of points; no interpolation.\n";
1527 return false;
1528 }
1529 if (iOrder == 0 || nx == 1) {
1530 // Zeroth order interpolation in x.
1531 // Find the nearest node.
1532 double dist = fabs(x - xAxis[0]);
1533 int iNode = 0;
1534 for (int i = 1; i < nx; i++) {
1535 if (fabs(x - xAxis[i]) < dist) {
1536 dist = fabs(x - xAxis[i]);
1537 iNode = i;
1538 }
1539 }
1540 // Set the summing range.
1541 iX0 = iNode;
1542 iX1 = iNode;
1543 // Establish the shape functions.
1544 fX = {1., 0., 0., 0.};
1545 } else if (iOrder == 1 || nx == 2) {
1546 // First order interpolation in x.
1547 // Find the grid segment containing this point.
1548 int iGrid = 0;
1549 for (int i = 1; i < nx; i++) {
1550 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1551 iGrid = i;
1552 }
1553 }
1554 const double x0 = xAxis[iGrid - 1];
1555 const double x1 = xAxis[iGrid];
1556 // Ensure there won't be divisions by zero.
1557 if (x1 == x0) {
1558 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1559 return false;
1560 }
1561 // Compute local coordinates.
1562 const double xL = (x - x0) / (x1 - x0);
1563 // Set the summing range.
1564 iX0 = iGrid - 1;
1565 iX1 = iGrid;
1566 // Set the shape functions.
1567 fX = {1. - xL, xL, 0., 0.};
1568 } else if (iOrder == 2) {
1569 // Second order interpolation in x.
1570 // Find the grid segment containing this point.
1571 int iGrid = 0;
1572 for (int i = 1; i < nx; i++) {
1573 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1574 iGrid = i;
1575 }
1576 }
1577 if (iGrid == 1) {
1578 iX0 = iGrid - 1;
1579 iX1 = iGrid + 1;
1580 const double x0 = xAxis[iX0];
1581 const double x1 = xAxis[iX0 + 1];
1582 const double x2 = xAxis[iX0 + 2];
1583 if (x0 == x1 || x0 == x2 || x1 == x2) {
1584 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1585 << " No interpolation.\n";
1586 return false;
1587 }
1588 fX[0] = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2));
1589 fX[1] = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2));
1590 fX[2] = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1));
1591 } else if (iGrid == nx - 1) {
1592 iX0 = iGrid - 2;
1593 iX1 = iGrid;
1594 const double x0 = xAxis[iX0];
1595 const double x1 = xAxis[iX0 + 1];
1596 const double x2 = xAxis[iX0 + 2];
1597 if (x0 == x1 || x0 == x2 || x1 == x2) {
1598 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1599 << " No interpolation.\n";
1600 return false;
1601 }
1602 fX[0] = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2));
1603 fX[1] = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2));
1604 fX[2] = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1));
1605 } else {
1606 iX0 = iGrid - 2;
1607 iX1 = iGrid + 1;
1608 const double x0 = xAxis[iX0];
1609 const double x1 = xAxis[iX0 + 1];
1610 const double x2 = xAxis[iX0 + 2];
1611 const double x3 = xAxis[iX0 + 3];
1612 if (x0 == x1 || x0 == x2 || x0 == x3 ||
1613 x1 == x2 || x1 == x3 || x2 == x3) {
1614 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1615 << " No interpolation.\n";
1616 return false;
1617 }
1618 // Compute the local coordinate for this grid segment.
1619 const double xL = (x - x1) / (x2 - x1);
1620 fX[0] = ((x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))) * (1. - xL);
1621 fX[1] = ((x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))) * (1. - xL) +
1622 ((x - x2) * (x - x3) / ((x1 - x2) * (x1 - x3))) * xL;
1623 fX[2] = ((x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))) * (1. - xL) +
1624 ((x - x1) * (x - x3) / ((x2 - x1) * (x2 - x3))) * xL;
1625 fX[3] = ((x - x1) * (x - x2) / ((x3 - x1) * (x3 - x2))) * xL;
1626 }
1627 }
1628
1629 if (iOrder == 0 || ny == 1) {
1630 // Zeroth order interpolation in y.
1631 // Find the nearest node.
1632 double dist = fabs(y - yAxis[0]);
1633 int iNode = 0;
1634 for (int i = 1; i < ny; i++) {
1635 if (fabs(y - yAxis[i]) < dist) {
1636 dist = fabs(y - yAxis[i]);
1637 iNode = i;
1638 }
1639 }
1640 // Set the summing range.
1641 iY0 = iNode;
1642 iY1 = iNode;
1643 // Establish the shape functions.
1644 fY = {1., 0., 0., 0.};
1645 } else if (iOrder == 1 || ny == 2) {
1646 // First order interpolation in y.
1647 // Find the grid segment containing this point.
1648 int iGrid = 0;
1649 for (int i = 1; i < ny; i++) {
1650 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
1651 iGrid = i;
1652 }
1653 }
1654 // Ensure there won't be divisions by zero.
1655 const double y0 = yAxis[iGrid - 1];
1656 const double y1 = yAxis[iGrid];
1657 if (y1 == y0) {
1658 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1659 return false;
1660 }
1661 // Compute local coordinates.
1662 const double yL = (y - y0) / (y1 - y0);
1663 // Set the summing range.
1664 iY0 = iGrid - 1;
1665 iY1 = iGrid;
1666 // Set the shape functions.
1667 fY = {1. - yL, yL, 0., 0.};
1668 } else if (iOrder == 2) {
1669 // Second order interpolation in y.
1670 // Find the grid segment containing this point.
1671 int iGrid = 0;
1672 for (int i = 1; i < ny; i++) {
1673 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
1674 iGrid = i;
1675 }
1676 }
1677 if (iGrid == 1) {
1678 iY0 = iGrid - 1;
1679 iY1 = iGrid + 1;
1680 const double y0 = yAxis[iY0];
1681 const double y1 = yAxis[iY0 + 1];
1682 const double y2 = yAxis[iY0 + 2];
1683 if (y0 == y1 || y0 == y2 || y1 == y2) {
1684 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1685 << " No interpolation.\n";
1686 return false;
1687 }
1688 fY[0] = (y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2));
1689 fY[1] = (y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2));
1690 fY[2] = (y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1));
1691 } else if (iGrid == ny - 1) {
1692 iY0 = iGrid - 2;
1693 iY1 = iGrid;
1694 const double y0 = yAxis[iY0];
1695 const double y1 = yAxis[iY0 + 1];
1696 const double y2 = yAxis[iY0 + 2];
1697 if (y0 == y1 || y0 == y2 || y1 == y2) {
1698 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1699 << " No interpolation.\n";
1700 return false;
1701 }
1702 fY[0] = (y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2));
1703 fY[1] = (y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2));
1704 fY[2] = (y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1));
1705 } else {
1706 iY0 = iGrid - 2;
1707 iY1 = iGrid + 1;
1708 const double y0 = yAxis[iY0];
1709 const double y1 = yAxis[iY0 + 1];
1710 const double y2 = yAxis[iY0 + 2];
1711 const double y3 = yAxis[iY0 + 3];
1712 if (y0 == y1 || y0 == y2 || y0 == y3 ||
1713 y1 == y2 || y1 == y3 || y2 == y3) {
1714 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1715 << " No interpolation.\n";
1716 return false;
1717 }
1718 // Compute the local coordinate for this grid segment.
1719 const double yL = (y - y1) / (y2 - y1);
1720 fY[0] = ((y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2))) * (1. - yL);
1721 fY[1] = ((y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2))) * (1. - yL) +
1722 ((y - y2) * (y - y3) / ((y1 - y2) * (y1 - y3))) * yL;
1723 fY[2] = ((y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1))) * (1. - yL) +
1724 ((y - y1) * (y - y3) / ((y2 - y1) * (y2 - y3))) * yL;
1725 fY[3] = ((y - y1) * (y - y2) / ((y3 - y1) * (y3 - y2))) * yL;
1726 }
1727 }
1728
1729 if (iOrder == 0 || nz == 1) {
1730 // Zeroth order interpolation in z.
1731 // Find the nearest node.
1732 double dist = fabs(z - zAxis[0]);
1733 int iNode = 0;
1734 for (int i = 1; i < nz; i++) {
1735 if (fabs(z - zAxis[i]) < dist) {
1736 dist = fabs(z - zAxis[i]);
1737 iNode = i;
1738 }
1739 }
1740 // Set the summing range.
1741 iZ0 = iNode;
1742 iZ1 = iNode;
1743 // Establish the shape functions.
1744 fZ = {1., 0., 0., 0.};
1745 } else if (iOrder == 1 || nz == 2) {
1746 // First order interpolation in z.
1747 // Find the grid segment containing this point.
1748 int iGrid = 0;
1749 for (int i = 1; i < nz; i++) {
1750 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1751 iGrid = i;
1752 }
1753 }
1754 const double z0 = zAxis[iGrid - 1];
1755 const double z1 = zAxis[iGrid];
1756 // Ensure there won't be divisions by zero.
1757 if (z1 == z0) {
1758 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1759 return false;
1760 }
1761 // Compute local coordinates.
1762 const double zL = (z - z0) / (z1 - z0);
1763 // Set the summing range.
1764 iZ0 = iGrid - 1;
1765 iZ1 = iGrid;
1766 // Set the shape functions.
1767 fZ = {1. - zL, zL, 0., 0.};
1768 } else if (iOrder == 2) {
1769 // Second order interpolation in z.
1770 // Find the grid segment containing this point.
1771 int iGrid = 0;
1772 for (int i = 1; i < nz; i++) {
1773 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1774 iGrid = i;
1775 }
1776 }
1777 if (iGrid == 1) {
1778 iZ0 = iGrid - 1;
1779 iZ1 = iGrid + 1;
1780 const double z0 = zAxis[iZ0];
1781 const double z1 = zAxis[iZ0 + 1];
1782 const double z2 = zAxis[iZ0 + 2];
1783 if (z0 == z1 || z0 == z2 || z1 == z2) {
1784 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1785 << " No interpolation.\n";
1786 return false;
1787 }
1788 fZ[0] = (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2));
1789 fZ[1] = (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2));
1790 fZ[2] = (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1));
1791 } else if (iGrid == nz - 1) {
1792 iZ0 = iGrid - 2;
1793 iZ1 = iGrid;
1794 const double z0 = zAxis[iZ0];
1795 const double z1 = zAxis[iZ0 + 1];
1796 const double z2 = zAxis[iZ0 + 2];
1797 if (z0 == z1 || z0 == z2 || z1 == z2) {
1798 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1799 << " No interpolation.\n";
1800 return false;
1801 }
1802 fZ[0] = (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2));
1803 fZ[1] = (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2));
1804 fZ[2] = (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1));
1805 } else {
1806 iZ0 = iGrid - 2;
1807 iZ1 = iGrid + 1;
1808 const double z0 = zAxis[iZ0];
1809 const double z1 = zAxis[iZ0 + 1];
1810 const double z2 = zAxis[iZ0 + 2];
1811 const double z3 = zAxis[iZ0 + 3];
1812 if (z0 == z1 || z0 == z2 || z0 == z3 ||
1813 z1 == z2 || z1 == z3 || z2 == z3) {
1814 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1815 << " No interpolation.\n";
1816 return false;
1817 }
1818 // Compute the local coordinate for this grid segment.
1819 const double zL = (z - z1) / (z2 - z1);
1820 fZ[0] = ((z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))) * (1. - zL);
1821 fZ[1] = ((z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))) * (1. - zL) +
1822 ((z - z2) * (z - z3) / ((z1 - z2) * (z1 - z3))) * zL;
1823 fZ[2] = ((z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))) * (1. - zL) +
1824 ((z - z1) * (z - z3) / ((z2 - z1) * (z2 - z3))) * zL;
1825 fZ[3] = ((z - z1) * (z - z2) / ((z3 - z1) * (z3 - z2))) * zL;
1826 }
1827 }
1828
1829 for (int i = iX0; i <= iX1; ++i) {
1830 for (int j = iY0; j <= iY1; ++j) {
1831 for (int k = iZ0; k <= iZ1; ++k) {
1832 f += value[i][j][k] * fX[i - iX0] * fY[j - iY0] * fZ[k - iZ0];
1833 }
1834 }
1835 }
1836 return true;
1837}

Referenced by Garfield::Medium::Interpolate().

◆ Divdif()

double Garfield::Numerics::Divdif ( const std::vector< double > &  f,
const std::vector< double > &  a,
int  nn,
double  x,
int  mm 
)

C++ version of DIVDIF (CERN program library E105) which performs tabular interpolation using symmetrically placed argument points.

Definition at line 1206 of file Numerics.cc.

1207 {
1208
1209 double t[20], d[20];
1210
1211 // Check the arguments.
1212 if (nn < 2) {
1213 std::cerr << "Divdif: Array length < 2.\n";
1214 return 0.;
1215 }
1216 if (mm < 1) {
1217 std::cerr << "Divdif: Interpolation order < 1.\n";
1218 return 0.;
1219 }
1220
1221 // Deal with the case that X is located at first or last point.
1222 const double tol1 = 1.e-6 * fabs(fabs(a[1]) - fabs(a[0]));
1223 const double tol2 = 1.e-6 * fabs(fabs(a[nn-1]) - fabs(a[nn-2]));
1224 if (fabs(x - a[0]) < tol1) return f[0];
1225 if (fabs(x - a[nn - 1]) < tol2) return f[nn - 1];
1226
1227 // Find subscript IX of X in array A.
1228 constexpr int mmax = 10;
1229 const int m = std::min({mm, mmax, nn - 1});
1230 const int mplus = m + 1;
1231 int ix = 0;
1232 int iy = nn + 1;
1233 if (a[0] > a[nn - 1]) {
1234 // Search decreasing arguments.
1235 do {
1236 const int mid = (ix + iy) / 2;
1237 if (x > a[mid - 1]) {
1238 iy = mid;
1239 } else {
1240 ix = mid;
1241 }
1242 } while (iy - ix > 1);
1243 } else {
1244 // Search increasing arguments.
1245 do {
1246 const int mid = (ix + iy) / 2;
1247 if (x < a[mid - 1]) {
1248 iy = mid;
1249 } else {
1250 ix = mid;
1251 }
1252 } while (iy - ix > 1);
1253 }
1254 // Copy reordered interpolation points into (T[I],D[I]), setting
1255 // EXTRA to True if M+2 points to be used.
1256 int npts = m + 2 - (m % 2);
1257 int ip = 0;
1258 int l = 0;
1259 do {
1260 const int isub = ix + l;
1261 if ((1 > isub) || (isub > nn)) {
1262 // Skip point.
1263 npts = mplus;
1264 } else {
1265 // Insert point.
1266 ip++;
1267 t[ip - 1] = a[isub - 1];
1268 d[ip - 1] = f[isub - 1];
1269 }
1270 if (ip < npts) {
1271 l = -l;
1272 if (l >= 0) ++l;
1273 }
1274 } while (ip < npts);
1275
1276 const bool extra = npts != mplus;
1277 // Replace d by the leading diagonal of a divided-difference table,
1278 // supplemented by an extra line if EXTRA is True.
1279 for (l = 1; l <= m; l++) {
1280 if (extra) {
1281 const int isub = mplus - l;
1282 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
1283 }
1284 int i = mplus;
1285 for (int j = l; j <= m; j++) {
1286 const int isub = i - l;
1287 d[i - 1] = (d[i - 1] - d[i - 2]) / (t[i - 1] - t[isub - 1]);
1288 i--;
1289 }
1290 }
1291 // Evaluate the Newton interpolation formula at X, averaging two values
1292 // of last difference if EXTRA is True.
1293 double sum = d[mplus - 1];
1294 if (extra) {
1295 sum = 0.5 * (sum + d[m + 1]);
1296 }
1297 int j = m;
1298 for (l = 1; l <= m; l++) {
1299 sum = d[j - 1] + (x - t[j - 1]) * sum;
1300 j--;
1301 }
1302 return sum;
1303}

Referenced by Garfield::Sensor::ComputeThresholdCrossings(), and Garfield::Medium::Interpolate1D().

◆ LeastSquaresFit()

bool Garfield::Numerics::LeastSquaresFit ( std::function< double(double, const std::vector< double > &)>  f,
std::vector< double > &  par,
std::vector< double > &  epar,
const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  ey,
const unsigned int  nMaxIter,
const double  diff,
double &  chi2,
const double  eps,
const bool  debug,
const bool  verbose 
)

Least-squares minimisation.

Definition at line 1839 of file Numerics.cc.

1845 {
1846
1847 //-----------------------------------------------------------------------
1848 // LSQFIT - Subroutine fitting the parameters A in the routine F to
1849 // the data points (X,Y) using a least squares method.
1850 // Translated from an Algol routine written by Geert Jan van
1851 // Oldenborgh and Rob Veenhof, based on Stoer + Bulirsch.
1852 // VARIABLES : F( . ,A,VAL) : Subroutine to be fitted.
1853 // (X,Y) : Input data.
1854 // D : Derivative matrix.
1855 // R : Difference vector between Y and F(X,A).
1856 // S : Correction vector for A.
1857 // EPSDIF : Used for differentiating.
1858 // EPS : Numerical resolution.
1859 // (Last updated on 23/ 5/11.)
1860 //-----------------------------------------------------------------------
1861
1862 const unsigned int n = par.size();
1863 const unsigned int m = x.size();
1864 // Make sure that the # degrees of freedom < the number of data points.
1865 if (n > m) {
1866 std::cerr << "LeastSquaresFit: Number of parameters to be varied\n"
1867 << " exceeds the number of data points.\n";
1868 return false;
1869 }
1870
1871 // Check the errors.
1872 if (*std::min_element(ey.cbegin(), ey.cend()) <= 0.) {
1873 std::cerr << "LeastSquaresFit: Not all errors are > 0; no fit done.\n";
1874 return false;
1875 }
1876 chi2 = 0.;
1877 // Largest difference.
1878 double diffc = -1.;
1879 // Initialise the difference vector R.
1880 std::vector<double> r(m, 0.);
1881 for (unsigned int i = 0; i < m; ++i) {
1882 // Compute initial residuals.
1883 r[i] = (y[i] - f(x[i], par)) / ey[i];
1884 // Compute initial maximum difference.
1885 diffc = std::max(std::abs(r[i]), diffc);
1886 // And compute initial chi2.
1887 chi2 += r[i] * r[i];
1888 }
1889 if (debug) {
1890 std::cout << " Input data points:\n"
1891 << " X Y Y - F(X)\n";
1892 for (unsigned int i = 0; i < m; ++i) {
1893 std::printf(" %9u %15.8e %15.8e %15.8e\n", i, x[i], y[i], r[i]);
1894 }
1895 std::cout << " Initial values of the fit parameters:\n"
1896 << " Parameter Value\n";
1897 for (unsigned int i = 0; i < n; ++i) {
1898 std::printf(" %9u %15.8e\n", i, par[i]);
1899 }
1900 }
1901 if (verbose) {
1902 std::cout << " MINIMISATION SUMMARY\n"
1903 << " Initial situation:\n";
1904 std::printf(" Largest difference between fit and target: %15.8e\n",
1905 diffc);
1906 std::printf(" Sum of squares of these differences: %15.8e\n",
1907 chi2);
1908 }
1909 // Start optimising loop.
1910 bool converged = false;
1911 double chi2L = 0.;
1912 for (unsigned int iter = 1; iter <= nMaxIter; ++iter) {
1913 // Check the stopping criteria: (1) max norm, (2) change in chi-squared.
1914 if ((diffc < diff) || (iter > 1 && std::abs(chi2L - chi2) < eps * chi2)) {
1915 if (debug || verbose) {
1916 if (diffc < diff) {
1917 std::cout << " The maximum difference stopping criterion "
1918 << "is satisfied.\n";
1919 } else {
1920 std::cout << " The relative change in chi-squared has dropped "
1921 << "below the threshold.\n";
1922 }
1923 }
1924 converged = true;
1925 break;
1926 }
1927 // Calculate the derivative matrix.
1928 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
1929 for (unsigned int i = 0; i < n; ++i) {
1930 const double epsdif = eps * (1. + std::abs(par[i]));
1931 par[i] += 0.5 * epsdif;
1932 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
1933 par[i] -= epsdif;
1934 for (unsigned int j = 0; j < m; ++j) {
1935 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
1936 }
1937 par[i] += 0.5 * epsdif;
1938 }
1939 // Invert the matrix in Householder style.
1940 std::vector<double> colsum(n, 0.);
1941 std::vector<int> pivot(n, 0);
1942 for (unsigned int i = 0; i < n; ++i) {
1943 colsum[i] = std::inner_product(d[i].cbegin(), d[i].cend(),
1944 d[i].cbegin(), 0.);
1945 pivot[i] = i;
1946 }
1947 // Decomposition.
1948 std::vector<double> alpha(n, 0.);
1949 bool singular = false;
1950 for (unsigned int k = 0; k < n; ++k) {
1951 double sigma = colsum[k];
1952 unsigned int jbar = k;
1953 for (unsigned int j = k + 1; j < n; ++j) {
1954 if (sigma < colsum[j]) {
1955 sigma = colsum[j];
1956 jbar = j;
1957 }
1958 }
1959 if (jbar != k) {
1960 // Interchange columns.
1961 std::swap(pivot[k], pivot[jbar]);
1962 std::swap(colsum[k], colsum[jbar]);
1963 std::swap(d[k], d[jbar]);
1964 }
1965 sigma = 0.;
1966 for (unsigned int i = k; i < m; ++i) sigma += d[k][i] * d[k][i];
1967 if (sigma == 0. || sqrt(sigma) < 1.e-8 * std::abs(d[k][k])) {
1968 singular = true;
1969 break;
1970 }
1971 alpha[k] = d[k][k] < 0. ? sqrt(sigma) : -sqrt(sigma);
1972 const double beta = 1. / (sigma - d[k][k] * alpha[k]);
1973 d[k][k] -= alpha[k];
1974 std::vector<double> b(n, 0.);
1975 for (unsigned int j = k + 1; j < n; ++j) {
1976 for (unsigned int i = k; i < n; ++i) b[j] += d[k][i] * d[j][i];
1977 b[j] *= beta;
1978 }
1979 for (unsigned int j = k + 1; j < n; ++j) {
1980 for (unsigned int i = k; i < m; ++i) {
1981 d[j][i] -= d[k][i] * b[j];
1982 colsum[j] -= d[j][k] * d[j][k];
1983 }
1984 }
1985 }
1986 if (singular) {
1987 std::cerr << "LeastSquaresFit: Householder matrix (nearly) singular;\n"
1988 << " no further optimisation.\n"
1989 << " Ensure the function depends on the parameters\n"
1990 << " and try to supply reasonable starting values.\n";
1991 break;
1992 }
1993 // Solve.
1994 for (unsigned int j = 0; j < n; ++j) {
1995 double gamma = 0.;
1996 for (unsigned int i = j; i < m; ++i) gamma += d[j][i] * r[i];
1997 gamma *= 1. / (alpha[j] * d[j][j]);
1998 for (unsigned int i = j; i < m; ++i) r[i] += gamma * d[j][i];
1999 }
2000 std::vector<double> z(n, 0.);
2001 z[n - 1] = r[n - 1] / alpha[n - 1];
2002 for (int i = n - 1; i >= 1; --i) {
2003 double sum = 0.;
2004 for (unsigned int j = i + 1; j <= n; ++j) {
2005 sum += d[j - 1][i - 1] * z[j - 1];
2006 }
2007 z[i - 1] = (r[i - 1] - sum) / alpha[i - 1];
2008 }
2009 // Correction vector.
2010 std::vector<double> s(n, 0.);
2011 for (unsigned int i = 0; i < n; ++i) s[pivot[i]] = z[i];
2012 // Generate some debugging output.
2013 if (debug) {
2014 std::cout << " Correction vector in iteration " << iter << ":\n";
2015 for (unsigned int i = 0; i < n; ++i) {
2016 std::printf(" %5u %15.8e\n", i, s[i]);
2017 }
2018 }
2019 // Add part of the correction vector to the estimate to improve chi2.
2020 chi2L = chi2;
2021 chi2 *= 2;
2022 for (unsigned int i = 0; i < n; ++i) par[i] += s[i] * 2;
2023 for (unsigned int i = 0; i <= 10; ++i) {
2024 if (chi2 <= chi2L) break;
2025 if (std::abs(chi2L - chi2) < eps * chi2) {
2026 if (debug) {
2027 std::cout << " Too little improvement; reduction loop halted.\n";
2028 }
2029 break;
2030 }
2031 chi2 = 0.;
2032 const double scale = 1. / pow(2, i);
2033 for (unsigned int j = 0; j < n; ++j) par[j] -= s[j] * scale;
2034 for (unsigned int j = 0; j < m; ++j) {
2035 r[j] = (y[j] - f(x[j], par)) / ey[j];
2036 chi2 += r[j] * r[j];
2037 }
2038 if (debug) {
2039 std::printf(" Reduction loop %3i: chi2 = %15.8e\n", i, chi2);
2040 }
2041 }
2042 // Calculate the max. norm.
2043 diffc = std::abs(r[0]);
2044 for (unsigned int i = 1; i < m; ++i) {
2045 diffc = std::max(std::abs(r[i]), diffc);
2046 }
2047 // Print some debugging output.
2048 if (debug) {
2049 std::cout << " Values of the fit parameters after iteration "
2050 << iter << "\n Parameter Value\n";
2051 for (unsigned int i = 0; i < n; ++i) {
2052 std::printf(" %9u %15.8e\n", i, par[i]);
2053 }
2054 std::printf(" for which chi2 = %15.8e and diff = %15.8e\n",
2055 chi2, diffc);
2056 } else if (verbose) {
2057 std::printf(" Step %3u: largest deviation = %15.8e, chi2 = %15.8e\n",
2058 iter, diffc, chi2);
2059 }
2060 }
2061 // End of fit, perform error calculation.
2062 if (!converged) {
2063 std::cerr << "LeastSquaresFit: Maximum number of iterations reached.\n";
2064 }
2065 // Calculate the derivative matrix for the final settings.
2066 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
2067 for (unsigned int i = 0; i < n; ++i) {
2068 const double epsdif = eps * (1. + std::abs(par[i]));
2069 par[i] += 0.5 * epsdif;
2070 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
2071 par[i] -= epsdif;
2072 for (unsigned int j = 0; j < m; ++j) {
2073 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
2074 }
2075 par[i] += 0.5 * epsdif;
2076 }
2077 // Calculate the error matrix.
2078 std::vector<std::vector<double> > cov(n, std::vector<double>(n, 0.));
2079 for (unsigned int i = 0; i < n; ++i) {
2080 for (unsigned int j = 0; j < n; ++j) {
2081 cov[i][j] = std::inner_product(d[i].cbegin(), d[i].cend(),
2082 d[j].cbegin(), 0.);
2083 }
2084 }
2085 // Compute the scaling factor for the errors.
2086 double scale = m > n ? chi2 / (m - n) : 1.;
2087 // Invert it to get the covariance matrix.
2088 epar.assign(n, 0.);
2089 if (Garfield::Numerics::CERNLIB::dinv(n, cov) != 0) {
2090 std::cerr << "LeastSquaresFit: Singular covariance matrix; "
2091 << "no error calculation.\n";
2092 } else {
2093 for (unsigned int i = 0; i < n; ++i) {
2094 for (unsigned int j = 0; j < n; ++j) cov[i][j] *= scale;
2095 epar[i] = sqrt(std::max(0., cov[i][i]));
2096 }
2097 }
2098 // Print results.
2099 if (debug) {
2100 std::cout << " Comparison between input and fit:\n"
2101 << " X Y F(X)\n";
2102 for (unsigned int i = 0; i < m; ++i) {
2103 std::printf(" %5u %15.8e %15.8e %15.8e\n", i, x[i], y[i], f(x[i], par));
2104 }
2105 }
2106 if (verbose) {
2107 std::cout << " Final values of the fit parameters:\n"
2108 << " Parameter Value Error\n";
2109 for (unsigned int i = 0; i < n; ++i) {
2110 std::printf(" %9u %15.8e %15.8e\n", i, par[i], epar[i]);
2111 }
2112 std::cout << " The errors have been scaled by a factor of "
2113 << sqrt(scale) << ".\n";
2114 std::cout << " Covariance matrix:\n";
2115 for (unsigned int i = 0; i < n; ++i) {
2116 for (unsigned int j = 0; j < n; ++j) {
2117 std::printf(" %15.8e", cov[i][j]);
2118 }
2119 std::cout << "\n";
2120 }
2121 std::cout << " Correlation matrix:\n";
2122 for (unsigned int i = 0; i < n; ++i) {
2123 for (unsigned int j = 0; j < n; ++j) {
2124 double cor = 0.;
2125 if (cov[i][i] > 0. && cov[j][j] > 0.) {
2126 cor = cov[i][j] / sqrt(cov[i][i] * cov[j][j]);
2127 }
2128 std::printf(" %15.8e", cor);
2129 }
2130 std::cout << "\n";
2131 }
2132 std::cout << " Minimisation finished.\n";
2133 }
2134 return converged;
2135}
Garfield::Numerics::CERNLIB::dinv
int dinv(const int n, std::vector< std::vector< double > > &a)
Replace square matrix A by its inverse.
Definition: Numerics.cc:787
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

Referenced by Garfield::ComponentAnalyticField::MultipoleMoments().

◆ Legendre()

double Garfield::Numerics::Legendre ( const unsigned int  n,
const double  x 
)
inline

Legendre polynomials.

Definition at line 119 of file Numerics.hh.

119 {
120
121 if (std::abs(x) > 1.) return 0.;
122 double p0 = 1.;
123 double p1 = x;
124 if (n == 0) return p0;
125 if (n == 1) return p1;
126 for (unsigned int k = 1; k < n; ++k) {
127 p0 = ((2 * k + 1) * x * p1 - k * p0) / (k + 1);
128 std::swap(p0, p1);
129 }
130 return p1;
131}

Referenced by Garfield::ComponentAnalyticField::MultipoleMoments().