Merge branch 'master' of https://github.com/AndreasFMueller/SeminarSpezielleFunktionen

4 files changed, 258 insertions, 0 deletions

diff --git a/buch/chapters/110-elliptisch/experiments/KK.pdf b/buch/chapters/110-elliptisch/experiments/KK.pdf
new file mode 100644
index 0000000..13a2739
--- /dev/null
+++ b/buch/chapters/110-elliptisch/experiments/KK.pdf
diff --git a/buch/chapters/110-elliptisch/experiments/KK.tex b/buch/chapters/110-elliptisch/experiments/KK.tex
new file mode 100644
index 0000000..a3ae425
--- /dev/null
+++ b/buch/chapters/110-elliptisch/experiments/KK.tex
@@ -0,0 +1,66 @@
+%
+% KK.tex -- template for standalon tikz images
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\documentclass[tikz]{standalone}
+\usepackage{amsmath}
+\usepackage{times}
+\usepackage{txfonts}
+\usepackage{pgfplots}
+\usepackage{csvsimple}
+\usetikzlibrary{arrows,intersections,math}
+\begin{document}
+\def\skala{1}
+\begin{tikzpicture}[>=latex,thick,scale=\skala]
+
+\def\dx{10}
+\def\dy{3}
+\input{KKpath.tex}
+
+\draw[->] (-0.1,0) -- (10.3,0) coordinate[label={$k$}];
+\draw[->] (0,-0.1) -- (0,{2*\dy+0.3}) coordinate[label={right:$y$}];
+
+\node at (3,{1.2*\dy}) {$\displaystyle y = \frac{K(k)}{K(\!\sqrt{1-k^2})}$};
+
+\begin{scope}
+\clip (0,0) rectangle (10,{2*\dy});
+\draw[color=red,line width=1.4pt] \KKpath;
+\end{scope}
+
+\draw[line width=0.2pt] (10,0) -- (10,{2*\dy});
+
+\foreach \y in {0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0}{
+	\draw (-0.05,{\y*\dy}) -- (0.05,{\y*\dy});
+	\node at (0,{\y*\dy}) [left] {$\y\mathstrut$};
+}
+
+\foreach \k in {1,...,9}{
+	\draw ({\k*\dx/10},-0.05) -- ({\k*\dx/10},0.05);
+	\node at ({\k*\dx/10},0) [below] {$0.\k\mathstrut$};
+}
+\node at (0,0) [below] {$0\mathstrut$};
+\node at (10,0) [below] {$1\mathstrut$};
+
+\draw[color=blue] ({\knull*\dx},0) -- ({\knull*\dx},{\KKnull*\dy});
+\foreach \y in {1,2,3,4}{
+	\draw[color=blue]
+		({\knull*\dx-0.05},{\y*\KKnull*\dy/5})
+		--
+		({\knull*\dx+0.05},{\y*\KKnull*\dy/5});
+}
+\draw[color=black,line width=0.1pt] (0,{\KKnull*\dy}) -- ({\knull*\dx},{\KKnull*\dy});
+\draw[color=black,line width=0.1pt] (0,{\KKnull*\dy/5}) -- ({\kone*\dx},{\KKnull*\dy/5});
+\node at ({0.6*\dx},{\KKnull*\dy}) [above] {$y=1.7732$};
+\node at ({0.6*\dx},{\KKnull*\dy/5}) [above] {$y=0.3546$};
+\draw[color=blue] ({\kone*\dx},0) -- ({\kone*\dx},{\KKnull*\dy/5});
+\draw[color=blue] ({\kone*\dx},{\KKnull*\dy/5}) -- ({\knull*\dx},{\KKnull*\dy/5});
+\fill[color=blue] ({\kone*\dx},{\KKnull*\dy/5}) circle[radius=0.05];
+\fill[color=blue] ({\knull*\dx},{\KKnull*\dy/5}) circle[radius=0.05];
+\fill[color=blue] ({\knull*\dx},{\KKnull*\dy}) circle[radius=0.05];
+\node[color=blue] at ({\knull*\dx},0) [left,rotate=90] {$k=0.97\mathstrut$};
+\node[color=blue] at ({\kone*\dx},0) [left,rotate=90] {$k_1=0.0477$};
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/buch/chapters/110-elliptisch/experiments/KN.cpp b/buch/chapters/110-elliptisch/experiments/KN.cpp
new file mode 100644
index 0000000..1dcca9e
--- /dev/null
+++ b/buch/chapters/110-elliptisch/experiments/KN.cpp
@@ -0,0 +1,177 @@
+/*
+ * KN.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <getopt.h>
+#include <vector>
+#include <gsl/gsl_sf_elljac.h>
+#include <gsl/gsl_sf_ellint.h>
+
+namespace KN {
+
+bool	debug = false;
+
+static struct option	longopts[] {
+{ "debug",	no_argument,		NULL,	'd' },
+{ "N",		required_argument,	NULL,	'N' },
+{ "outfile",	required_argument,	NULL,	'o' },
+{ "min",	required_argument,	NULL,	'm' },
+{ NULL,		0,			NULL,	 0  }
+};
+
+double	KprimeK(double k) {
+	double	kprime = sqrt(1-k*k);
+	if (debug)
+		printf("%s:%d: k = %f, k' = %f\n", __FILE__, __LINE__, k, kprime);
+	double	v
+		=
+		gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE)
+		/
+		gsl_sf_ellint_Kcomp(kprime, GSL_PREC_DOUBLE)
+		;
+	if (debug)
+		printf("%s:%d: KprimeK(k = %f) = %f\n", __FILE__, __LINE__, k, v);
+	return v;
+}
+
+static const int	L = 100000000;
+static const double	h = 1. / L;
+
+double	Kd(double k) {
+	double	m = 0;
+	if (k < h) {
+		m = 2 * (KprimeK(k) - KprimeK(k / 2)) / k;
+	} else if (k > 1-h) {
+		m = 2 * (KprimeK((1 + k) / 2) - KprimeK(k)) / (1 - k);
+		
+	} else {
+		m = L * (KprimeK(k + h) - KprimeK(k));
+	}
+	if (debug)
+		printf("%s:%d: Kd(%f) = %f\n", __FILE__, __LINE__, k, m);
+	return m;
+}
+
+double	k1(double y) {
+	if (debug)
+		printf("%s:%d: Newton for y = %f\n", __FILE__, __LINE__, y);
+	double	kn = 0.5;
+	double	delta = 1;
+	int	n = 0;
+	while ((fabs(delta) > 0.000001) && (n < 10)) {
+		double	yn = KprimeK(kn);
+		if (debug)
+			printf("%s:%d: k%d = %f, y%d = %f\n", __FILE__, __LINE__, n, kn, n, yn);
+		delta = (yn - y) / Kd(kn);
+		if (debug)
+			printf("%s:%d: delta = %f\n", __FILE__, __LINE__, delta);
+		double	kneu = kn - delta;
+		if (kneu <= 0) {
+			kneu = kn / 4;
+		}
+		if (kneu >= 1) {
+			kneu = (3 + kn) / 4;
+		}
+		kn = kneu;
+		if (debug)
+			printf("%s:%d: kneu = %f, kn = %f\n", __FILE__, __LINE__, kneu, kn);
+		n++;
+	}
+	if (debug)
+		printf("%s:%d: Newton result: k = %f\n", __FILE__, __LINE__, kn);
+	return kn;
+}
+
+double	k1(int N, double k) {
+	return k1(KprimeK(k) / N);
+}
+
+/**
+ * \brief Main function for the slcl program
+ */
+int	main(int argc, char *argv[]) {
+	int	longindex;
+	int	c;
+	int	N = 5;
+	double	kmin = 0.01;
+	std::string	outfilename;
+	while (EOF != (c = getopt_long(argc, argv, "d:N:o:m:",
+		longopts, &longindex)))
+		switch (c) {
+		case 'd':
+			debug = true;
+			break;
+		case 'N':
+			N = std::stoi(optarg);
+			break;
+		case 'o':
+			outfilename = std::string(optarg);
+			break;
+		case 'm':
+			kmin = std::stod(optarg);
+			break;
+		}
+
+	double	d = 0.01;
+	if (outfilename.size() > 0) {
+		FILE	*fn = fopen(outfilename.c_str(), "w");
+		fprintf(fn, "\\def\\KKpath{ ");
+		double	k = d;
+		fprintf(fn, " (0,0)");
+		double	k0 = k/16;
+		while (k0 < k) {
+			fprintf(fn, "\n\t-- ({%.4f*\\dx},{%.4f*\\dy})", k0, KprimeK(k0));
+			k0 *= 2;
+		}
+		while (k < 1-0.5*d) {
+			fprintf(fn, "\n\t-- ({%.4f*\\dx},{%.4f*\\dy})", k, KprimeK(k));
+			k += d;
+		}
+		fprintf(fn, "}\n");
+
+		k0 = 0.97;
+		fprintf(fn, "\\def\\knull{%.4f}\n", k0);
+		double	KK = KprimeK(k0);
+		fprintf(fn, "\\def\\KKnull{%.4f}\n", KK);
+		fprintf(fn, "\\def\\kone{%.4f}\n", k1(N, k0));
+
+		fclose(fn);
+		return EXIT_SUCCESS;
+	}
+
+	for (double k = kmin; k < (1 - d/2); k += d) {
+		if (debug)
+			printf("%s:%d: k = %f\n", __FILE__, __LINE__, k);
+		double	y = KprimeK(k);
+		double	k0 = k1(y);
+		double	kone = k1(N, k0);
+		printf("g(%4.2f) = %10.6f,", k, y);
+		printf("    g'(%.2f) = %10.6f,", k, Kd(k));
+		printf("    g^{-1} = %10.6f,", k0);
+		printf("    k1 = %10.6f,", kone);
+		printf("    g(k1) = %10.6f\n", KprimeK(kone));
+	}
+
+	return EXIT_SUCCESS;
+}
+
+} // namespace KN
+
+int	main(int argc, char *argv[]) {
+	try {
+		return KN::main(argc, argv);
+	} catch (const std::exception& e) {
+		std::cerr << "terminated by exception: " << e.what();
+		std::cerr << std::endl;
+	} catch (...) {
+		std::cerr << "terminated by unknown exception" << std::endl;
+	}
+	return EXIT_FAILURE;
+}
diff --git a/buch/chapters/110-elliptisch/experiments/Makefile b/buch/chapters/110-elliptisch/experiments/Makefile
new file mode 100644
index 0000000..fac4fbc
--- /dev/null
+++ b/buch/chapters/110-elliptisch/experiments/Makefile
@@ -0,0 +1,15 @@
+#
+# Makefile
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschlue
+#
+all:	KK.pdf
+
+KN:	KN.cpp
+	g++ -O -Wall -std=c++11 KN.cpp -o KN `pkg-config --cflags gsl` `pkg-config --libs gsl`
+
+KKpath.tex:	KN
+	./KN --outfile KKpath.tex
+
+KK.pdf:	KK.tex KKpath.tex
+	pdflatex KK.tex