diff options
Diffstat (limited to '')
| -rw-r--r-- | buch/chapters/110-elliptisch/experiments/KN.cpp | 177 |
1 files changed, 177 insertions, 0 deletions
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; +} |