aboutsummaryrefslogtreecommitdiffstats
path: root/buch/chapters/110-elliptisch/images/jacobi.cpp
blob: 438fe4894a59fdfd153b2fa7e97855358fa84c50 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
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;
}