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
|
/*
* airy.cpp
*
* (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
*/
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <string>
#include <getopt.h>
#include <gsl/gsl_sf_hyperg.h>
bool debug = false;
struct option options[] = {
{ "debug", no_argument, 0, 'd' },
{ "outfile", required_argument, 0, 'o' },
{ "a", required_argument, 0, 'a' },
{ "b", required_argument, 0, 'b' },
{ "s", required_argument, 0, 's' },
{ NULL, 0, 0, 0 }
};
double h0f1(double c, double x) {
if (debug)
fprintf(stderr, "%s:%d: c = %.5f, x = %.5f\n",
__FILE__, __LINE__, c, x);
double s = 0;
int k = 0;
double term = 1.;
s += term;
int counter = 0;
while (fabs(term) > 0.000000000001) {
k++;
term *= x / ((c+k-1) * k);
s += term;
counter++;
}
if (debug)
fprintf(stderr, "x = %f, terms = %d\n", x, counter);
return s;
}
double f1(double x) {
// unfortunately, gsl_sf_hyperg_0F1 does not work if c<1, because
// it uses a relation to the bessel functions
//return gsl_sf_hyperg_0F1(2/3, x*x*x/9.);
return h0f1(2./3., x*x*x/9.);
}
double f2(double x) {
return x * gsl_sf_hyperg_0F1(4./3., x*x*x/9.);
}
double f3(double x) {
return x * h0f1(4./3., x*x*x/9.);
}
void plot(FILE *outfile, const char *name, double (*f)(double x),
double a, double b, int steps) {
fprintf(outfile, "\\def\\%s{\n", name);
fprintf(outfile, "({%.5f*\\dx},{%.5f*\\dy})", a, f(a));
double h = (b-a)/steps;
for (int i = 0; i <= steps; i++) {
double x = a + h*i;
if (debug)
fprintf(stderr, "x = %.5f\n", x);
fprintf(outfile, "\n\t--({%.5f*\\dx},{%.5f*\\dy})",
x, f(x));
}
fprintf(outfile, "}\n");
}
int main(int argc, char *argv[]) {
std::string outfilename("airypaths.tex");
int c;
int longindex;
double a = -8;
double b = 2.5;
int steps = 500;
while (EOF != (c = getopt_long(argc, argv, "a:b:do:s:",
options, &longindex)))
switch (c) {
case 'a':
a = std::stod(optarg);
break;
case 'b':
b = std::stod(optarg);
break;
case 'd':
debug = true;
break;
case 'o':
outfilename = std::string(optarg);
break;
case 's':
steps = std::stoi(optarg);
break;
}
if (debug)
fprintf(stderr, "%s:%d: outfile: '%s', a = %.4f, b = %.4f, steps = %d\n",
__FILE__, __LINE__, outfilename.c_str(), a, b, steps);
FILE *outfile = fopen(outfilename.c_str(), "w");
plot(outfile, "yonepath", f1, a, b, steps);
plot(outfile, "ytwopath", f2, a, b, steps);
plot(outfile, "ythreepath", f3, a, b, steps);
fclose(outfile);
return EXIT_SUCCESS;
}
|
