diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindRuleFactory.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindRuleFactory.java
new file mode 100644
index 000000000..23dff1c06
--- /dev/null
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindRuleFactory.java
@@ -0,0 +1,70 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.integration.gauss;
+
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Pair;
+
+/**
+ * Factory that creates Gauss-Chebyshev quadrature rules of the second kind.
+ *
+ * The generated rules approximate integrals of the form:
+ *
+ *
+ * integral from -1 to 1 of f(x) sqrt(1 - x^2) dx
+ *
+ *
+ * using:
+ *
+ *
+ * sum from i = 1 to n of w_i f(x_i)
+ *
+ *
+ * where:
+ *
+ *
+ * x_i = cos(i pi / (n + 1)), i = 1, ..., n
+ * w_i = pi / (n + 1) sin^2(i pi / (n + 1))
+ *
+ */
+public class ChebyshevSecondKindRuleFactory extends AbstractRuleFactory {
+
+ /**
+ * Computes the Gauss-Chebyshev quadrature rule of the second kind.
+ *
+ * @param numberOfPoints order of the rule to be computed
+ * @return nodes and weights of the quadrature rule
+ */
+ @Override
+ protected Pair computeRule(final int numberOfPoints) {
+ final double[] points = new double[numberOfPoints];
+ final double[] weights = new double[numberOfPoints];
+
+ final double scale = FastMath.PI / (numberOfPoints + 1.0);
+
+ for (int i = 0; i < numberOfPoints; i++) {
+ final double angle = (i + 1.0) * scale;
+ final double sin = FastMath.sin(angle);
+ final int index = numberOfPoints - 1 - i;
+
+ points[index] = FastMath.cos(angle);
+ weights[index] = scale * sin * sin;
+ }
+
+ return new Pair<>(points, weights);
+ }
+}
diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindRuleFactory.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindRuleFactory.java
new file mode 100644
index 000000000..7046ed7e3
--- /dev/null
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindRuleFactory.java
@@ -0,0 +1,86 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.integration.gauss;
+
+import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
+import org.hipparchus.util.Pair;
+
+/**
+ * Factory that creates field-based Gauss-Chebyshev quadrature rules of the
+ * second kind.
+ *
+ * The generated rules approximate integrals of the form:
+ *
+ *
+ * integral from -1 to 1 of f(x) sqrt(1 - x^2) dx
+ *
+ *
+ * using:
+ *
+ *
+ * sum from i = 1 to n of w_i f(x_i)
+ *
+ *
+ * where:
+ *
+ *
+ * x_i = cos(i pi / (n + 1)), i = 1, ..., n
+ * w_i = pi / (n + 1) sin^2(i pi / (n + 1))
+ *
+ *
+ * @param type of the field elements
+ */
+public class FieldChebyshevSecondKindRuleFactory>
+ extends FieldAbstractRuleFactory {
+
+ /**
+ * Constructor.
+ *
+ * @param field field to which rule coefficients belong
+ */
+ public FieldChebyshevSecondKindRuleFactory(final Field field) {
+ super(field);
+ }
+
+ /**
+ * Computes the Gauss-Chebyshev quadrature rule of the second kind.
+ *
+ * @param numberOfPoints order of the rule to be computed
+ * @return nodes and weights of the quadrature rule
+ */
+ @Override
+ protected Pair computeRule(final int numberOfPoints) {
+ final T[] points = MathArrays.buildArray(getField(), numberOfPoints);
+ final T[] weights = MathArrays.buildArray(getField(), numberOfPoints);
+
+ final T scale = getField().getZero().newInstance(FastMath.PI / (numberOfPoints + 1.0));
+
+ for (int i = 0; i < numberOfPoints; i++) {
+ final T angle = scale.multiply(i + 1.0);
+ final T sin = FastMath.sin(angle);
+ final int index = numberOfPoints - 1 - i;
+
+ points[index] = FastMath.cos(angle);
+ weights[index] = scale.multiply(sin.square());
+ }
+
+ return new Pair<>(points, weights);
+ }
+}
diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldGaussIntegratorFactory.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldGaussIntegratorFactory.java
index 9f4d83a57..7fcf021de 100644
--- a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldGaussIntegratorFactory.java
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/FieldGaussIntegratorFactory.java
@@ -40,6 +40,8 @@ public class FieldGaussIntegratorFactory> {
private final FieldRuleFactory hermite;
/** Generator of Gauss-Laguerre integrators. */
private final FieldRuleFactory laguerre;
+ /** Generator of Gauss-Chebyshev integrators. */
+ private final FieldRuleFactory chebyshevSecondKind;
/** Simple constructor.
* @param field field to which function argument and value belong
@@ -48,6 +50,7 @@ public FieldGaussIntegratorFactory(final Field field) {
legendre = new FieldLegendreRuleFactory<>(field);
hermite = new FieldHermiteRuleFactory<>(field);
laguerre = new FieldLaguerreRuleFactory<>(field);
+ chebyshevSecondKind = new FieldChebyshevSecondKindRuleFactory<>(field);
}
/**
@@ -120,6 +123,20 @@ public SymmetricFieldGaussIntegrator hermite(int numberOfPoints) {
return new SymmetricFieldGaussIntegrator<>(hermite.getRule(numberOfPoints));
}
+ /**
+ * Creates a Gauss-Chebyshev integrator of the second kind of the given order.
+ * The call to the
+ * {@link FieldGaussIntegrator#integrate(org.hipparchus.analysis.CalculusFieldUnivariateFunction)
+ * integrate} method will perform an integration on the natural interval
+ * {@code [-1 , 1]}.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a Gauss-Chebyshev integrator.
+ */
+ public FieldGaussIntegrator chebyshevSecondKind(int numberOfPoints) {
+ return new FieldGaussIntegrator<>(chebyshevSecondKind.getRule(numberOfPoints));
+ }
+
/**
* Performs a change of variable so that the integration can be performed
* on an arbitrary interval {@code [a, b]}.
diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/GaussIntegratorFactory.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/GaussIntegratorFactory.java
index 4a8fa29b1..1a1008afd 100644
--- a/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/GaussIntegratorFactory.java
+++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/integration/gauss/GaussIntegratorFactory.java
@@ -42,6 +42,8 @@ public class GaussIntegratorFactory {
private final RuleFactory hermite;
/** Generator of Gauss-Laguerre integrators. */
private final RuleFactory laguerre;
+ /** Generator of Gauss-Chebyshev integrators. */
+ private final RuleFactory chebyshevSecondKind;
/** Simple constructor.
*/
@@ -57,6 +59,7 @@ public GaussIntegratorFactory(final int decimalDigits) {
legendreHighPrecision = new ConvertingRuleFactory<>(new FieldLegendreRuleFactory<>(new DfpField(decimalDigits)));
hermite = new HermiteRuleFactory();
laguerre = new LaguerreRuleFactory();
+ chebyshevSecondKind = new ChebyshevSecondKindRuleFactory();
}
/**
@@ -146,6 +149,21 @@ public GaussIntegrator legendreHighPrecision(int numberOfPoints,
lowerBound, upperBound));
}
+ /**
+ * Creates a Gauss-Chebyshev integrator of the second kind of the given order.
+ * The call to the
+ * {@link GaussIntegrator#integrate(org.hipparchus.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the natural interval
+ * {@code [-1 , 1]}.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a GaussChebyshev integrator.
+ */
+ public GaussIntegrator chebyshevSecondKind(int numberOfPoints)
+ throws MathIllegalArgumentException {
+ return new GaussIntegrator(chebyshevSecondKind.getRule(numberOfPoints));
+ }
+
/**
* Creates a Gauss-Hermite integrator of the given order.
* The call to the
diff --git a/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindTest.java b/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindTest.java
new file mode 100644
index 000000000..3328ab38f
--- /dev/null
+++ b/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/ChebyshevSecondKindTest.java
@@ -0,0 +1,329 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.integration.gauss;
+
+import org.hipparchus.special.Erf;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Pair;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+/**
+ * Tests Gauss-Chebyshev quadrature rules of the second kind against reference
+ * data generated by the {@code GenGCQuad} MATLAB function from NASA's CARA
+ * Analysis Tools.
+ *
+ * @see NASA CARA Analysis Tools
+ */
+class ChebyshevSecondKindTest {
+
+ private static final double TOLERANCE = 1.0e-15;
+
+ static final double[][] GEN_GC_QUAD_16 = {
+ {-9.829730996839018e-01, 1.837495178165702e-01, 6.773407894247465e-03},
+ {-9.324722294043558e-01, 3.612416661871530e-01, 1.331615550646369e-02},
+ {-8.502171357296142e-01, 5.264321628773557e-01, 1.940543741387547e-02},
+ {-7.390089172206590e-01, 6.736956436465573e-01, 2.483389042441458e-02},
+ {-6.026346363792563e-01, 7.980172272802396e-01, 2.941665508151885e-02},
+ {-4.457383557765382e-01, 8.951632913550623e-01, 3.299767083121100e-02},
+ {-2.736629900720829e-01, 9.618256431728190e-01, 3.545499047709023e-02},
+ {-9.226835946330189e-02, 9.957341762950346e-01, 3.670493294584622e-02},
+ {9.226835946330202e-02, 9.957341762950345e-01, 3.670493294584622e-02},
+ {2.736629900720828e-01, 9.618256431728190e-01, 3.545499047709023e-02},
+ {4.457383557765383e-01, 8.951632913550623e-01, 3.299767083121100e-02},
+ {6.026346363792564e-01, 7.980172272802395e-01, 2.941665508151885e-02},
+ {7.390089172206591e-01, 6.736956436465573e-01, 2.483389042441458e-02},
+ {8.502171357296142e-01, 5.264321628773557e-01, 1.940543741387547e-02},
+ {9.324722294043558e-01, 3.612416661871530e-01, 1.331615550646369e-02},
+ {9.829730996839018e-01, 1.837495178165702e-01, 6.773407894247465e-03}
+ };
+
+ static final double[][] GEN_GC_QUAD_64 = {
+ {-9.988322268323265e-01, 4.831337952550822e-02, 4.657833967754515e-04},
+ {-9.953316347176486e-01, 9.651392091451513e-02, 9.304789348454766e-04},
+ {-9.895063994510511e-01, 1.444890495692213e-01, 1.393001296249116e-03},
+ {-9.813701261394134e-01, 1.921267173537084e-01, 1.852270238580170e-03},
+ {-9.709418174260519e-01, 2.393156642875583e-01, 2.307213117943439e-03},
+ {-9.582458291091662e-01, 2.859456783986892e-01, 2.756767394164220e-03},
+ {-9.433118132577431e-01, 3.319078531285286e-01, 3.199883112400166e-03},
+ {-9.261746489577764e-01, 3.770948416883209e-01, 3.635525355359408e-03},
+ {-9.068743608505453e-01, 4.214011077725294e-01, 4.062676660397876e-03},
+ {-8.854560256532098e-01, 4.647231720437687e-01, 4.480339395850453e-03},
+ {-8.619696668800491e-01, 5.069598538135907e-01, 4.887538091045943e-03},
+ {-8.364701380102265e-01, 5.480125073546703e-01, 5.283321714564022e-03},
+ {-8.090169943749473e-01, 5.877852522924733e-01, 5.666765895413191e-03},
+ {-7.796743540632223e-01, 6.261851975383138e-01, 6.036975081942061e-03},
+ {-7.485107481711009e-01, 6.631226582407955e-01, 6.393084633441722e-03},
+ {-7.155989607441211e-01, 6.985113652489370e-01, 6.734262839554183e-03},
+ {-6.810158587867969e-01, 7.322686665977737e-01, 7.059712862770465e-03},
+ {-6.448422127361704e-01, 7.643157205458485e-01, 7.368674599481501e-03},
+ {-6.071625078187112e-01, 7.945776797137543e-01, 7.660426455235353e-03},
+ {-5.680647467311557e-01, 8.229838658936565e-01, 7.934287030054488e-03},
+ {-5.276402441061325e-01, 8.494679351215212e-01, 8.189616709876989e-03},
+ {-4.859834132426061e-01, 8.739680326265179e-01, 8.425819160404840e-03},
+ {-4.431915455992412e-01, 8.964269372957038e-01, 8.642342719870315e-03},
+ {-3.993645835656953e-01, 9.167921953165825e-01, 8.838681687467587e-03},
+ {-3.546048870425355e-01, 9.350162426854148e-01, 9.014377504440393e-03},
+ {-3.090169943749473e-01, 9.510565162951535e-01, 9.169019825067272e-03},
+ {-2.627073781985868e-01, 9.648755533435515e-01, 9.302247475042994e-03},
+ {-2.157841967678060e-01, 9.764410788292721e-01, 9.413749295017889e-03},
+ {-1.683570413470385e-01, 9.857260809316509e-01, 9.503264867324932e-03},
+ {-1.205366802553229e-01, 9.927088740980540e-01, 9.570585124197263e-03},
+ {-7.243480016176228e-02, 9.973731496914912e-01, 9.615552836055650e-03},
+ {-2.416374523613207e-02, 9.997080140801929e-01, 9.638062978725454e-03},
+ {2.416374523613242e-02, 9.997080140801929e-01, 9.638062978725454e-03},
+ {7.243480016176240e-02, 9.973731496914912e-01, 9.615552836055650e-03},
+ {1.205366802553232e-01, 9.927088740980540e-01, 9.570585124197263e-03},
+ {1.683570413470386e-01, 9.857260809316508e-01, 9.503264867324930e-03},
+ {2.157841967678064e-01, 9.764410788292720e-01, 9.413749295017888e-03},
+ {2.627073781985870e-01, 9.648755533435515e-01, 9.302247475042994e-03},
+ {3.090169943749475e-01, 9.510565162951535e-01, 9.169019825067272e-03},
+ {3.546048870425358e-01, 9.350162426854147e-01, 9.014377504440393e-03},
+ {3.993645835656957e-01, 9.167921953165823e-01, 8.838681687467584e-03},
+ {4.431915455992413e-01, 8.964269372957038e-01, 8.642342719870315e-03},
+ {4.859834132426062e-01, 8.739680326265179e-01, 8.425819160404840e-03},
+ {5.276402441061328e-01, 8.494679351215211e-01, 8.189616709876989e-03},
+ {5.680647467311559e-01, 8.229838658936564e-01, 7.934287030054487e-03},
+ {6.071625078187113e-01, 7.945776797137542e-01, 7.660426455235352e-03},
+ {6.448422127361706e-01, 7.643157205458483e-01, 7.368674599481499e-03},
+ {6.810158587867972e-01, 7.322686665977735e-01, 7.059712862770463e-03},
+ {7.155989607441211e-01, 6.985113652489370e-01, 6.734262839554183e-03},
+ {7.485107481711012e-01, 6.631226582407951e-01, 6.393084633441719e-03},
+ {7.796743540632225e-01, 6.261851975383136e-01, 6.036975081942058e-03},
+ {8.090169943749475e-01, 5.877852522924731e-01, 5.666765895413190e-03},
+ {8.364701380102267e-01, 5.480125073546700e-01, 5.283321714564020e-03},
+ {8.619696668800493e-01, 5.069598538135903e-01, 4.887538091045939e-03},
+ {8.854560256532099e-01, 4.647231720437685e-01, 4.480339395850450e-03},
+ {9.068743608505454e-01, 4.214011077725293e-01, 4.062676660397874e-03},
+ {9.261746489577766e-01, 3.770948416883203e-01, 3.635525355359402e-03},
+ {9.433118132577432e-01, 3.319078531285284e-01, 3.199883112400164e-03},
+ {9.582458291091662e-01, 2.859456783986892e-01, 2.756767394164220e-03},
+ {9.709418174260520e-01, 2.393156642875579e-01, 2.307213117943434e-03},
+ {9.813701261394134e-01, 1.921267173537084e-01, 1.852270238580170e-03},
+ {9.895063994510511e-01, 1.444890495692213e-01, 1.393001296249116e-03},
+ {9.953316347176486e-01, 9.651392091451513e-02, 9.304789348454766e-04},
+ {9.988322268323266e-01, 4.831337952550593e-02, 4.657833967754294e-04}
+ };
+
+ @Test
+ void matchesGenGCQuadSixteenPointTableAfterPcCircleWeightConversion() {
+ assertMatchesGenGCQuadTable(GEN_GC_QUAD_16);
+ }
+
+ @Test
+ void matchesGenGCQuadSixtyFourPointTableAfterPcCircleWeightConversion() {
+ assertMatchesGenGCQuadTable(GEN_GC_QUAD_64);
+ }
+
+ /**
+ * Verifies that the generated quadrature rule matches one of the reference
+ * tables in the {@code GenGCQuad} MATLAB function after converting the
+ * standard second-kind weights to the table's PcCircle weight convention.
+ *
+ * @param genGcQuadTable reference table from {@code GenGCQuad}
+ */
+ private static void assertMatchesGenGCQuadTable(final double[][] genGcQuadTable) {
+ final int order = genGcQuadTable.length;
+
+ final ChebyshevSecondKindRuleFactory factory = new ChebyshevSecondKindRuleFactory();
+ final Pair rule = factory.getRule(order);
+
+ final double[] secondKindPoints = rule.getFirst();
+ final double[] secondKindWeights = rule.getSecond();
+
+ Assertions.assertEquals(order, secondKindPoints.length);
+ Assertions.assertEquals(order, secondKindWeights.length);
+
+ for (int tableIndex = 0; tableIndex < order; tableIndex++) {
+
+ final double expectedX = genGcQuadTable[tableIndex][0];
+ final double expectedY = genGcQuadTable[tableIndex][1];
+ final double expectedPcCircleWeight = genGcQuadTable[tableIndex][2];
+
+ final double x = secondKindPoints[tableIndex];
+ final double y = FastMath.sqrt(1.0 - x * x);
+
+ /*
+ * Standard second-kind weight:
+ *
+ * W = pi / (n + 1) * sin^2(theta)
+ *
+ * PcCircle-customized table weight:
+ *
+ * wGC = pi / (n + 1) * sin(theta) / sqrt(8*pi)
+ *
+ * Since y = sin(theta):
+ *
+ * wGC = W / y / sqrt(8*pi)
+ */
+ final double pcCircleWeight = secondKindWeights[tableIndex] / y / FastMath.sqrt(8.0 * FastMath.PI);
+
+ Assertions.assertEquals(expectedX, x, TOLERANCE, "xGC mismatch at table index " + tableIndex);
+ Assertions.assertEquals(expectedY, y, TOLERANCE, "yGC mismatch at table index " + tableIndex);
+ Assertions.assertEquals(expectedPcCircleWeight, pcCircleWeight, TOLERANCE,
+ "wGC mismatch at table index " + tableIndex);
+ }
+ }
+
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForModerateOffsetWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.35,
+ 0.20,
+ 1.10,
+ 0.75,
+ 0.45,
+ 16,
+ 1.06055364621781625e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForModerateOffsetWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.35,
+ 0.20,
+ 1.10,
+ 0.75,
+ 0.45,
+ 64,
+ 1.06055364621781625e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForNearCenterWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.05,
+ 0.08,
+ 0.80,
+ 1.25,
+ 0.30,
+ 16,
+ 4.37336618192281021e-02
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForNearCenterWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.05,
+ 0.08,
+ 0.80,
+ 1.25,
+ 0.30,
+ 64,
+ 4.37336618192281021e-02
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForLargerRadiusWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.90,
+ 0.40,
+ 1.50,
+ 0.65,
+ 0.70,
+ 16,
+ 1.56273946267817959e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForLargerRadiusWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.90,
+ 0.40,
+ 1.50,
+ 0.65,
+ 0.70,
+ 64,
+ 1.56273946267818015e-01
+ );
+ }
+
+ private static void assertMatchesMatlabReference(final double xm,
+ final double zm,
+ final double sx,
+ final double sz,
+ final double hbr,
+ final int order,
+ final double expected) {
+ final double actual = integratePcCircleStyle(xm, zm, sx, sz, hbr, order);
+
+ Assertions.assertEquals(expected, actual, TOLERANCE,
+ "PcCircle-style second-kind quadrature mismatch for order " + order);
+ }
+
+ /**
+ * Uses the quadrature rule in the PcCircle convention to test probability of
+ * collision calculations against MATLAB results for simple reference cases.
+ *
+ * @param xm projected miss distance along the x-axis
+ * @param zm projected miss distance along the z-axis
+ * @param sx standard deviation along the x-axis
+ * @param sz standard deviation along the z-axis
+ * @param hbr hard-body radius
+ * @param order quadrature rule order
+ * @return probability of collision estimate
+ */
+ private static double integratePcCircleStyle(final double xm,
+ final double zm,
+ final double sx,
+ final double sz,
+ final double hbr,
+ final int order) {
+ final GaussIntegratorFactory factory = new GaussIntegratorFactory();
+ final GaussIntegrator integrator = factory.chebyshevSecondKind(order);
+
+ final double sqrt2 = FastMath.sqrt(2.0);
+ final double dx = sqrt2 * sx;
+ final double dz = sqrt2 * sz;
+
+ final double integral = integrator.integrate(u -> {
+ final double y = FastMath.sqrt(1.0 - u * u);
+ final double x = xm + hbr * u;
+ final double hx = hbr * y;
+
+ final double exponential = FastMath.exp(-FastMath.pow(x / dx, 2.0));
+ final double errorFunctionDifference = Erf.erf((zm - hx) / dz, (zm + hx) / dz);
+
+ /*
+ * PcCircle's GenGCQuad convention uses:
+ *
+ * wGC = pi / (n + 1) * y / sqrt(8*pi)
+ *
+ * while the standard second-kind rule uses:
+ *
+ * W = pi / (n + 1) * y^2
+ *
+ * Therefore the integrand passed to GaussIntegrator must include:
+ *
+ * 1 / (y * sqrt(8*pi))
+ *
+ * so that:
+ *
+ * W * F(u) / (y * sqrt(8*pi)) == wGC * F(u)
+ */
+ return exponential * errorFunctionDifference / (y * FastMath.sqrt(8.0 * FastMath.PI));
+ });
+
+ return hbr / sx * integral;
+ }
+}
\ No newline at end of file
diff --git a/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindTest.java b/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindTest.java
new file mode 100644
index 000000000..e86c73f87
--- /dev/null
+++ b/hipparchus-core/src/test/java/org/hipparchus/analysis/integration/gauss/FieldChebyshevSecondKindTest.java
@@ -0,0 +1,255 @@
+/*
+ * Licensed to the Hipparchus project under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * https://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.hipparchus.analysis.integration.gauss;
+
+import org.hipparchus.special.Erf;
+import org.hipparchus.util.Binary64;
+import org.hipparchus.util.Binary64Field;
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.Pair;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+/**
+ * Tests field-based Gauss-Chebyshev quadrature rules of the second kind against
+ * reference data generated by the {@code GenGCQuad} MATLAB function from NASA's
+ * CARA Analysis Tools.
+ *
+ * @see NASA CARA Analysis Tools
+ */
+class FieldChebyshevSecondKindTest {
+
+ private static final double TOLERANCE = 1.0e-15;
+
+ @Test
+ void matchesGenGCQuadSixteenPointTableAfterPcCircleWeightConversion() {
+ assertMatchesGenGCQuadTable(ChebyshevSecondKindTest.GEN_GC_QUAD_16);
+ }
+
+ @Test
+ void matchesGenGCQuadSixtyFourPointTableAfterPcCircleWeightConversion() {
+ assertMatchesGenGCQuadTable(ChebyshevSecondKindTest.GEN_GC_QUAD_64);
+ }
+
+ /**
+ * Verifies that the generated quadrature rule matches one of the reference
+ * tables in the {@code GenGCQuad} MATLAB function after converting the
+ * standard second-kind weights to the table's PcCircle weight convention.
+ *
+ * @param genGcQuadTable reference table from {@code GenGCQuad}
+ */
+ private static void assertMatchesGenGCQuadTable(final double[][] genGcQuadTable) {
+ final int order = genGcQuadTable.length;
+
+ final FieldChebyshevSecondKindRuleFactory factory =
+ new FieldChebyshevSecondKindRuleFactory<>(Binary64Field.getInstance());
+ final Pair rule = factory.getRule(order);
+
+ final Binary64[] secondKindPoints = rule.getFirst();
+ final Binary64[] secondKindWeights = rule.getSecond();
+
+ Assertions.assertEquals(order, secondKindPoints.length);
+ Assertions.assertEquals(order, secondKindWeights.length);
+
+ for (int tableIndex = 0; tableIndex < order; tableIndex++) {
+ final double expectedX = genGcQuadTable[tableIndex][0];
+ final double expectedY = genGcQuadTable[tableIndex][1];
+ final double expectedPcCircleWeight = genGcQuadTable[tableIndex][2];
+
+ final Binary64 x = secondKindPoints[tableIndex];
+ final Binary64 y = FastMath.sqrt(Binary64.ONE.subtract(x.square()));
+
+ /*
+ * Standard second-kind weight:
+ *
+ * W = pi / (n + 1) * sin^2(theta)
+ *
+ * PcCircle-customized table weight:
+ *
+ * wGC = pi / (n + 1) * sin(theta) / sqrt(8*pi)
+ *
+ * Since y = sin(theta):
+ *
+ * wGC = W / y / sqrt(8*pi)
+ */
+ final Binary64 pcCircleWeight =
+ secondKindWeights[tableIndex].divide(y).divide(FastMath.sqrt(8.0 * FastMath.PI));
+
+ Assertions.assertEquals(expectedX, x.getReal(), TOLERANCE, "xGC mismatch at table index " + tableIndex);
+ Assertions.assertEquals(expectedY, y.getReal(), TOLERANCE, "yGC mismatch at table index " + tableIndex);
+ Assertions.assertEquals(expectedPcCircleWeight, pcCircleWeight.getReal(), TOLERANCE,
+ "wGC mismatch at table index " + tableIndex);
+ }
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForModerateOffsetWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.35,
+ 0.20,
+ 1.10,
+ 0.75,
+ 0.45,
+ 16,
+ 1.06055364621781625e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForModerateOffsetWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.35,
+ 0.20,
+ 1.10,
+ 0.75,
+ 0.45,
+ 64,
+ 1.06055364621781625e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForNearCenterWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.05,
+ 0.08,
+ 0.80,
+ 1.25,
+ 0.30,
+ 16,
+ 4.37336618192281021e-02
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForNearCenterWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.05,
+ 0.08,
+ 0.80,
+ 1.25,
+ 0.30,
+ 64,
+ 4.37336618192281021e-02
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForLargerRadiusWithOrder16() {
+ assertMatchesMatlabReference(
+ 0.90,
+ 0.40,
+ 1.50,
+ 0.65,
+ 0.70,
+ 16,
+ 1.56273946267817959e-01
+ );
+ }
+
+ @Test
+ void matchesMatlabPcCircleIntegrandForLargerRadiusWithOrder64() {
+ assertMatchesMatlabReference(
+ 0.90,
+ 0.40,
+ 1.50,
+ 0.65,
+ 0.70,
+ 64,
+ 1.56273946267818015e-01
+ );
+ }
+
+ private static void assertMatchesMatlabReference(final double xm,
+ final double zm,
+ final double sx,
+ final double sz,
+ final double hbr,
+ final int order,
+ final double expected) {
+ final Binary64 actual = integratePcCircleStyle(
+ new Binary64(xm),
+ new Binary64(zm),
+ new Binary64(sx),
+ new Binary64(sz),
+ new Binary64(hbr),
+ order
+ );
+
+ Assertions.assertEquals(expected, actual.getReal(), TOLERANCE,
+ "PcCircle-style field second-kind quadrature mismatch for order " + order);
+ }
+
+ /**
+ * Uses the quadrature rule in the PcCircle convention to test probability of
+ * collision calculations against MATLAB results for simple reference cases.
+ *
+ * @param xm projected miss distance along the x-axis
+ * @param zm projected miss distance along the z-axis
+ * @param sx standard deviation along the x-axis
+ * @param sz standard deviation along the z-axis
+ * @param hbr hard-body radius
+ * @param order quadrature rule order
+ * @return probability of collision estimate
+ */
+ private static Binary64 integratePcCircleStyle(final Binary64 xm,
+ final Binary64 zm,
+ final Binary64 sx,
+ final Binary64 sz,
+ final Binary64 hbr,
+ final int order) {
+ final FieldGaussIntegratorFactory factory =
+ new FieldGaussIntegratorFactory<>(Binary64Field.getInstance());
+ final FieldGaussIntegrator integrator = factory.chebyshevSecondKind(order);
+
+ final Binary64 sqrt2 = FastMath.sqrt(new Binary64(2.0));
+ final Binary64 dx = sx.multiply(sqrt2);
+ final Binary64 dz = sz.multiply(sqrt2);
+
+ final Binary64 integral = integrator.integrate(u -> {
+ final Binary64 y = FastMath.sqrt(Binary64.ONE.subtract(u.square()));
+ final Binary64 x = xm.add(hbr.multiply(u));
+ final Binary64 hx = hbr.multiply(y);
+
+ final Binary64 exponential = FastMath.exp(x.divide(dx).square().negate());
+ final Binary64 errorFunctionDifference = Erf.erf(zm.subtract(hx).divide(dz),
+ zm.add(hx).divide(dz));
+
+ /*
+ * PcCircle's GenGCQuad convention uses:
+ *
+ * wGC = pi / (n + 1) * y / sqrt(8*pi)
+ *
+ * while the standard second-kind rule uses:
+ *
+ * W = pi / (n + 1) * y^2
+ *
+ * Therefore the integrand passed to FieldGaussIntegrator must include:
+ *
+ * 1 / (y * sqrt(8*pi))
+ *
+ * so that:
+ *
+ * W * F(u) / (y * sqrt(8*pi)) == wGC * F(u)
+ */
+ return exponential.multiply(errorFunctionDifference)
+ .divide(y.multiply(FastMath.sqrt(8.0 * FastMath.PI)));
+ });
+
+ return hbr.divide(sx).multiply(integral);
+ }
+}
\ No newline at end of file