From 05d75b0f467b2535db538ecaee461cf0c8b637d1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andreas=20M=C3=BCller?= <andreas.mueller@ost.ch>
Date: Mon, 27 Jun 2022 20:17:16 +0200
Subject: add stuff for elliptic filters

---
 buch/chapters/110-elliptisch/agm/sn.cpp | 52 +++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)
 create mode 100644 buch/chapters/110-elliptisch/agm/sn.cpp

(limited to 'buch/chapters/110-elliptisch/agm/sn.cpp')

diff --git a/buch/chapters/110-elliptisch/agm/sn.cpp b/buch/chapters/110-elliptisch/agm/sn.cpp
new file mode 100644
index 0000000..ff2ed17
--- /dev/null
+++ b/buch/chapters/110-elliptisch/agm/sn.cpp
@@ -0,0 +1,52 @@
+/*
+ * ns.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <gsl/gsl_sf_ellint.h>
+#include <gsl/gsl_sf_elljac.h>
+
+static const int	N = 10;
+
+inline long double	sqrl(long double x) {
+	return x * x;
+}
+
+int	main(int argc, char *argv[]) {
+	long double	u = 0.6;
+	long double	k = 0.9;
+	long double	kprime = sqrt(1 - sqrl(k));
+	
+	long double	a[N], b[N], x[N+1];
+	a[0] = 1;
+	b[0] = kprime;
+
+	for (int n = 0; n < N-1; n++) {
+		printf("a[%d] = %22.18Lf  b[%d] = %22.18Lf\n", n, a[n], n, b[n]);
+		a[n+1] = (a[n] + b[n]) / 2;
+		b[n+1] = sqrtl(a[n] * b[n]);
+	}
+
+	x[N] = sinl(u * a[N-1]);
+	printf("x[%d] = %22.18Lf\n", N, x[N]);
+
+	for (int n = N - 1; n >= 0; n--) {
+		x[n] = 2 * a[n] * x[n+1] / (a[n] + b[n] + (a[n] - b[n]) * sqrl(x[n+1]));
+		printf("x[%2d] = %22.18Lf\n", n, x[n]);
+	}
+	
+	printf("sn(%7.4Lf, %7.4Lf) = %20.24Lf\n", u, k, x[0]);
+
+	double	sn, cn, dn;
+	double	m = sqrl(k);
+	gsl_sf_elljac_e((double)u, m, &sn, &cn, &dn);
+	printf("sn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, sn);
+	printf("cn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, cn);
+	printf("dn(%7.4Lf, %7.4Lf) = %20.24f\n", u, k, dn);
+
+	return EXIT_SUCCESS;
+}
-- 
cgit v1.2.1