diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/specialfunctions/Elliptic.java b/src/main/java/com/wildbitsfoundry/etk4j/math/specialfunctions/Elliptic.java new file mode 100644 index 0000000..d87b692 --- /dev/null +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/specialfunctions/Elliptic.java @@ -0,0 +1,161 @@ +package com.wildbitsfoundry.etk4j.math.specialfunctions; + +public final class Elliptic { + + /** + * Complete Elliptic integral of the first kind ({@code K}). + * @param k The parameter of the elliptic integral. The absolute value of {@code K(k2} has to be less + * than 1. + * @return The complete elliptic integral of the first kind {@code K(k)}. + * @see Elliptic Integrals + */ + public static double completeEllipticIntegralFirstKind(double k) { + double sum; + double term; + double above; + double below; + sum = 1; + term = 1; + + above = 1; + below = 2; + + for (int i = 1; i <= 100; i++) { + term *= above / below; + sum += Math.pow(k, i) * Math.pow(term, 2); + above += 2; + below += 2; + } + sum *= 0.5 * Math.PI; + return sum; + } + + /** + * Complete Elliptic integral of the second kind ({@code E}). + * @param k The parameter of the elliptic integral. The absolute value of {@code E(k2} has to be less + * than 1. + * @return The complete elliptic integral of the second kind {@code E(k)}. + * @see Elliptic Integrals + */ + public static double completeEllipticIntegralSecondKind(double k) { + double sum; + double term; + double above; + double below; + sum = 1; + term = 1; + + above = 1; + below = 2; + + for (int i = 1; i <= 100; i++) { + term *= above / below; + sum -= Math.pow(k, i) * Math.pow(term, 2) / above; + above += 2; + below += 2; + } + sum *= 0.5 * Math.PI; + return sum; + } + + /** + * Incomplete Elliptic integral of the first kind ({@code Kinc}). + * @param angle The angle in rad/s. The angle should be between {@code 0} and {@code π/2} + * @param k The parameter is taken as {@code k2}. + * @return The incomplete elliptic integral of the first kind {@code Kinc(k)}. + * @see Elliptic Integrals + */ + public static double incompleteEllipticIntegralFirstKind(double angle, double k) { + double result; + result = Math.sin(angle) * RF(Math.pow(Math.cos(angle), 2), 1 - k * Math.pow(Math.sin(angle), 2), 1); + return result; + } + + /** + * Incomplete Elliptic integral of the second kind ({@code Einc}). + * @param angle The angle in rad/s. The angle should be between {@code 0} and {@code π/2} + * @param k The parameter is taken as {@code k2}. + * @return The incomplete elliptic integral of the first kind {@code Einc(k)}. + * @see Elliptic Integrals + */ + public static double incompleteEllipticIntegralSecondKind(double angle, double k) { + double y = 1 - k * Math.pow(Math.sin(angle), 2); + return Math.sin(angle) * RF(Math.pow(Math.cos(angle), 2), y, 1) + (-1.0 / 3.0) * k * Math.pow(Math.sin(angle), + (3.0)) * RD(Math.pow(Math.cos(angle), 2), y, 1); + } + + /** + * Carlson's Elliptic integral of the first kind + * @see Series Expansion + */ + private static double RF(double x, double y, double z) { + double dx, dy, dz; + double lambda; + double n = 1.0 / 3.0; + double mean; + double tmp; + do { + lambda = Math.sqrt(x * y) + Math.sqrt(y * z) + Math.sqrt(z * x); + x = 0.25 * (x + lambda); + y = 0.25 * (y + lambda); + z = 0.25 * (z + lambda); + mean = (x + y + z) * n; + tmp = 1 / mean; + dx = (mean - x) * tmp; + dy = (mean - y) * tmp; + dz = (mean - z) * tmp; + } while (Math.max(Math.max(Math.abs(dx), Math.abs(dy)), Math.abs(dz)) > 1e-7); + double e2 = dx * dy - dz * dz; + double e3 = dx * dy * dz; + double c1 = 1.0 / 24.0; + double c2 = 0.1; + double c3 = 3.0 / 44.0; + double c4 = 1.0 / 14.0; + + double result = 1.0 + (c1 * e2 - c2 - c3 * e3) * e2 + c4 * e3; + return result / Math.sqrt(mean); + } + + /** + * Carlson's Elliptic integral of the second kind + * @see Series Expansion + */ + private static double RD(double x, double y, double z) { + double dx, dy, dz; + double lambda; + double mu; + double muInv; + double sum = 0.0; + double pow4 = 1.0; + + do { + lambda = Math.sqrt(x * y) + Math.sqrt(y * z) + Math.sqrt(z * x); + sum += pow4 / (Math.sqrt(z) * (z + lambda)); + + pow4 *= 0.25; + + x = 0.25 * (x + lambda); + y = 0.25 * (y + lambda); + z = 0.25 * (z + lambda); + mu = (x + y + 3.0 * z) * 0.2; + muInv = 1.0 / mu; + + dx = 1 - x * muInv; + dy = 1 - y * muInv; + dz = 1 - z * muInv; + } while (Math.max(Math.max(Math.abs(dx), Math.abs(dy)), Math.abs(dz)) > 1e-7); + double C1 = 3.0 / 14.0; + double C2 = 1.0 / 6.0; + double C3 = 9.0 / 22.0; + double C4 = 3.0 / 26.0; + double EA = dx * dy; + double EB = dz * dz; + double EC = EA - EB; + double ED = EA - 6.0 * EB; + double EF = ED + EC + EC; + double S1 = ED * (-C1 + 0.25 * C3 * ED - 1.50 * C4 * dz * EF); + double S2 = dz * (C2 * EF + dz * (-C3 * EC + dz * C4 * EA)); + + return 3.0 * sum + pow4 * (1.0 + S1 + S2) / (mu * Math.sqrt(mu)); + } +} \ No newline at end of file diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/specialfunctions/EllipticTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/specialfunctions/EllipticTest.java new file mode 100644 index 0000000..5605b4c --- /dev/null +++ b/src/test/java/com/wildbitsfoundry/etk4j/math/specialfunctions/EllipticTest.java @@ -0,0 +1,42 @@ +package com.wildbitsfoundry.etk4j.math.specialfunctions; + +import org.junit.Test; + +import static org.junit.Assert.assertEquals; + +public class EllipticTest { + + @Test + public void testCompleteOfFirstKind() { + double k = 0.1; + double expected = 1.6124413487202198; + double y = Elliptic.completeEllipticIntegralFirstKind(k); + + assertEquals(expected, y, 1e-12); + } + + @Test + public void testIncompleteOfFirstKind() { + double expected = 1.2261911708835167; + double y = Elliptic.incompleteEllipticIntegralFirstKind(1.0, 1.0); + + assertEquals(expected, y, 1e-12); + } + + @Test + public void testCompleteOfSecondKind() { + double k = 0.1; + double expected = 1.5307576368977631; + double y = Elliptic.completeEllipticIntegralSecondKind(k); + + assertEquals(expected, y, 1e-12); + } + + @Test + public void testIncompleteOfSecondKind() { + double expected = 0.8414709848078963; + double y = Elliptic.incompleteEllipticIntegralSecondKind(1.0, 1.0); + + assertEquals(expected, y, 1e-12); + } +}