From 24e0ebb5475b275af76900be4cbe5d7a8df0ef06 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andreas=20M=C3=BCller?= <andreas.mueller@ost.ch>
Date: Mon, 18 Oct 2021 08:51:09 +0200
Subject: improve elliptic integrals

---
 buch/chapters/110-elliptisch/images/rechteck.cpp |  53 +++++++++++++++++++++--
 buch/chapters/110-elliptisch/images/rechteck.pdf | Bin 70950 -> 91630 bytes
 2 files changed, 50 insertions(+), 3 deletions(-)

(limited to 'buch/chapters/110-elliptisch/images')

diff --git a/buch/chapters/110-elliptisch/images/rechteck.cpp b/buch/chapters/110-elliptisch/images/rechteck.cpp
index 3ece11c..c3f6f1b 100644
--- a/buch/chapters/110-elliptisch/images/rechteck.cpp
+++ b/buch/chapters/110-elliptisch/images/rechteck.cpp
@@ -243,6 +243,7 @@ void	curvedrawer::operator()(const curvetracer::curve_t& curve) {
 static struct option	longopts[] = {
 { "outfile",	required_argument,	NULL,	'o' },
 { "k",		required_argument,	NULL,	'k' },
+{ "deltax",	required_argument,	NULL,	'd' },
 { NULL,		0,			NULL,	 0  }
 };
 
@@ -250,13 +251,18 @@ static struct option	longopts[] = {
  * \brief Main function
  */
 int	main(int argc, char *argv[]) {
-	double	k = 0.6;
+	double	k = 0.625;
+	double	deltax = 0.2;
 
 	int	c;
 	int	longindex;
 	std::string	outfilename;
-	while (EOF != (c = getopt_long(argc, argv, "o:k:", longopts, &longindex)))
+	while (EOF != (c = getopt_long(argc, argv, "o:k:d:", longopts,
+		&longindex)))
 		switch (c) {
+		case 'd':
+			deltax = std::stod(optarg);
+			break;
 		case 'o':
 			outfilename = std::string(optarg);
 			break;
@@ -296,7 +302,7 @@ int	main(int argc, char *argv[]) {
 
 	// "circles"
 	std::complex<double>	dir(0.01, 0);
-	for (double im = 0.2; im < 3; im += 0.2) {
+	for (double im = deltax; im < 3; im += deltax) {
 		std::complex<double>	startz(0, im);
 		std::complex<double>	startw = ct.startpoint(startz);
 		curvetracer::curve_t	curve = ct.trace(startz, dir,
@@ -345,7 +351,28 @@ int	main(int argc, char *argv[]) {
 		(*cdp)(curvetracer::mirrory(curve));
 	}
 
+	// arguments between 1 and 1/k
+	{
+		for (double x0 = 1 + deltax; x0 < 1/k; x0 += deltax) {
+			double	y0 = sqrt(1-1/(x0*x0))/kprime;
+			//std::cout << "y0 = " << y0 << std::endl;
+			double	y = gsl_sf_ellint_F(asin(y0), kprime,
+				GSL_PREC_DOUBLE);
+			std::complex<double>	startz(x0);
+			std::complex<double>	startw(xmax, y);
+			curvetracer::curve_t	curve = ct.trace(startz, dir,
+				startw, 1000);
+			cdp->color("red");
+			(*cdp)(curve);
+			(*cdp)(curvetracer::mirrorx(curve));
+			cdp->color("blue");
+			(*cdp)(curvetracer::mirror(curve));
+			(*cdp)(curvetracer::mirrory(curve));
+		}
+	}
+
 	// argument 1/k
+#if 0
 	{
 		std::complex<double>	startz(1/k);
 		std::complex<double>	startw(xmax, ymax);
@@ -358,6 +385,26 @@ int	main(int argc, char *argv[]) {
 		(*cdp)(curvetracer::mirror(curve));
 		(*cdp)(curvetracer::mirrory(curve));
 	}
+#endif
+
+	// arguments larger than 1/k
+	{
+		double	x0 = 1;
+		while (x0 <= 1/k + 0.0001) { x0 += deltax; }
+		for (; x0 < 4; x0 += deltax) {
+			std::complex<double>	startz(x0);
+			std::complex<double>	startw(gsl_sf_ellint_F(
+				asin(1/(k*x0)), k, GSL_PREC_DOUBLE), ymax);
+		curvetracer::curve_t	curve = ct.trace(startz, dir,
+			startw, 1000);
+		cdp->color("red");
+		(*cdp)(curve);
+		(*cdp)(curvetracer::mirrorx(curve));
+		cdp->color("blue");
+		(*cdp)(curvetracer::mirror(curve));
+		(*cdp)(curvetracer::mirrory(curve));
+		}
+	}
 
 	// border
 	(*cdp->out()) << "\\def\\xmax{" << xmax << "}" << std::endl;
diff --git a/buch/chapters/110-elliptisch/images/rechteck.pdf b/buch/chapters/110-elliptisch/images/rechteck.pdf
index 5fa07cb..845a550 100644
Binary files a/buch/chapters/110-elliptisch/images/rechteck.pdf and b/buch/chapters/110-elliptisch/images/rechteck.pdf differ
-- 
cgit v1.2.1