1 files changed, 154 insertions, 0 deletions
diff --git a/buch/chapters/110-elliptisch/images/jacobi.cpp b/buch/chapters/110-elliptisch/images/jacobi.cpp
new file mode 100644
index 0000000..438fe48
--- /dev/null
+++ b/buch/chapters/110-elliptisch/images/jacobi.cpp
@@ -0,0 +1,154 @@
+/*
+ * jacobi.cpp
+ *
+ * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+#include <gsl/gsl_sf_elljac.h>
+#include <vector>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <getopt.h>
+
+namespace jacobi {
+
+static struct option	longopts[] {
+{ "k",		required_argument,	NULL,	'k' },
+{ "m",		required_argument,	NULL,	'm' },
+{ "outfile",	required_argument,	NULL,	'o' },
+{ "infile",	required_argument,	NULL,	'i' },
+{ "umin",	required_argument,	NULL,	'f' },
+{ "umax",	required_argument,	NULL,	't' },
+{ "steps",	required_argument,	NULL,	's' },
+{ NULL,		0,			NULL,	 0  }
+};
+
+/**
+ * \brief A class to contain a single plot
+ */
+class jacobiplot : public std::vector<std::pair<double,double> > {
+	std::string	punkt(std::pair<double,double> x) const {
+		char	buffer[128];
+		snprintf(buffer, sizeof(buffer), "({%.4f*\\dx},{%.4f*\\dy})",
+			x.first, x.second);
+		return std::string(buffer);
+	}
+public:
+	jacobiplot() { }
+	void	write(std::ostream& out, const std::string& name) const {
+		out << "\\def\\" << name << "{" << std::endl;
+		auto	i = begin();
+		out << punkt(*i);
+		for (i++; i != end(); i++) {
+			out << std::endl << " -- " << punkt(*i);
+		}
+		out << std::endl << "}" << std::endl;
+	}
+};
+
+/**
+ * \brief A class containing all three basic functions
+ */
+class jacobiplots {
+	jacobiplot	_sn;
+	jacobiplot	_cn;
+	jacobiplot	_dn;
+	double	_m;
+public:
+	const jacobiplot&	sn() const { return _sn; }
+	const jacobiplot&	cn() const { return _cn; }
+	const jacobiplot&	dn() const { return _dn; }
+	jacobiplots(double m) : _m(m) { }
+	double	m() const { return _m; }
+	void	add(double umin, double umax, int steps) {
+		_sn.clear();
+		_cn.clear();
+		_dn.clear();
+		double	h = (umax - umin) / steps;
+		for (int i = 0; i <= steps; i++) {
+			double	u = umin + i * h;
+			double	snvalue, cnvalue, dnvalue;
+			gsl_sf_elljac_e(u, m(), &snvalue, &cnvalue, &dnvalue);
+			_sn.push_back(std::make_pair(u, snvalue));
+			_cn.push_back(std::make_pair(u, cnvalue));
+			_dn.push_back(std::make_pair(u, dnvalue));
+		}
+	}
+	void	write(std::ostream& out, const std::string& name) const {
+		out << "% name=" << name << ", m=" << m() << std::endl;
+		out << "\\def\\m" << name << "{" << m() << "}" << std::endl;
+		_sn.write(out, std::string("sn") + name);
+		_cn.write(out, std::string("cn") + name);
+		_dn.write(out, std::string("dn") + name);
+	}
+};
+
+int	main(int argc, char *argv[]) {
+	int	c;
+	int	longindex;
+	std::ostream	*out = &std::cout;
+	std::string	infilename;
+	int	steps = 100;
+	double	umin = 0;
+	double	umax = 10;
+	while (EOF != (c = getopt_long(argc, argv, "o:i:f:t:s:", longopts,
+		&longindex)))
+		switch (c) {
+		case 'o':
+			out = new std::ofstream(optarg);
+			(*out) << "%" << std::endl;
+			(*out) << "% " << optarg << std::endl;
+			(*out) << "%" << std::endl;
+			break;
+		case 'i':
+			infilename = std::string(optarg);
+			break;
+		case 'f':
+			umin = std::stod(optarg);
+			break;
+		case 't':
+			umax = std::stod(optarg);
+			break;
+		case 's':
+			steps = std::stoi(optarg);
+			break;
+		}
+
+	if (infilename.size() == 0) {
+		while ((optind + 1) < argc) {
+			std::string	name(argv[optind++]);
+			double	m = std::stod(argv[optind++]);
+			jacobiplots	j(m);
+			j.add(umin, umax, steps);
+			j.write(*out, name);
+		}
+	} else {
+		std::ifstream	infile(infilename);
+		while (!infile.eof()) {
+			std::string	name;
+			double	m;
+			infile >> name;
+			infile >> m;
+			jacobiplots	j(m);
+			j.add(umin, umax, steps);
+			j.write(*out, name);
+		}
+	}
+
+	return EXIT_SUCCESS;
+}
+
+} // namespace jacobi
+
+int	main(int argc, char *argv[]) {
+	try {
+		return jacobi::main(argc, argv);
+	} catch (const std::exception& x) {
+		std::cerr << "cannot build sn/cn/dn file: " << x.what();
+		std::cerr << std::endl;
+	}
+	return EXIT_FAILURE;
+}