aboutsummaryrefslogtreecommitdiffstats
path: root/buch/papers/kugel/figures/povray/spherecurve.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'buch/papers/kugel/figures/povray/spherecurve.cpp')
-rw-r--r--buch/papers/kugel/figures/povray/spherecurve.cpp292
1 files changed, 292 insertions, 0 deletions
diff --git a/buch/papers/kugel/figures/povray/spherecurve.cpp b/buch/papers/kugel/figures/povray/spherecurve.cpp
new file mode 100644
index 0000000..8ddf5e5
--- /dev/null
+++ b/buch/papers/kugel/figures/povray/spherecurve.cpp
@@ -0,0 +1,292 @@
+/*
+ * spherecurve.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <string>
+#include <iostream>
+
+inline double sqr(double x) { return x * x; }
+
+/**
+ * \brief Class for 3d vectors (also used as colors)
+ */
+class vector {
+ double X[3];
+public:
+ vector() { X[0] = X[1] = X[2] = 0; }
+ vector(double a) { X[0] = X[1] = X[2] = a; }
+ vector(double x, double y, double z) {
+ X[0] = x; X[1] = y; X[2] = z;
+ }
+ vector(double theta, double phi) {
+ double s = sin(theta);
+ X[0] = cos(phi) * s;
+ X[1] = sin(phi) * s;
+ X[2] = cos(theta);
+ }
+ vector(const vector& other) {
+ for (int i = 0; i < 3; i++) {
+ X[i] = other.X[i];
+ }
+ }
+ vector operator+(const vector& other) const {
+ return vector(X[0] + other.X[0],
+ X[1] + other.X[1],
+ X[2] + other.X[2]);
+ }
+ vector operator*(double l) const {
+ return vector(X[0] * l, X[1] * l, X[2] * l);
+ }
+ double operator*(const vector& other) const {
+ double s = 0;
+ for (int i = 0; i < 3; i++) {
+ s += X[i] * other.X[i];
+ }
+ return s;
+ }
+ double norm() const {
+ double s = 0;
+ for (int i = 0; i < 3; i++) {
+ s += sqr(X[i]);
+ }
+ return sqrt(s);
+ }
+ vector normalize() const {
+ double l = norm();
+ return vector(X[0]/l, X[1]/l, X[2]/l);
+ }
+ double max() const {
+ return std::max(X[0], std::max(X[1], X[2]));
+ }
+ double l0norm() const {
+ double l = 0;
+ for (int i = 0; i < 3; i++) {
+ if (fabs(X[i]) > l) {
+ l = fabs(X[i]);
+ }
+ }
+ return l;
+ }
+ vector l0normalize() const {
+ double l = l0norm();
+ vector result(X[0]/l, X[1]/l, X[2]/l);
+ return result;
+ }
+ const double& operator[](int i) const { return X[i]; }
+ double& operator[](int i) { return X[i]; }
+};
+
+/**
+ * \brief Derived 3d vector class implementing color
+ *
+ * The constructor in this class converts a single value into a
+ * color on a suitable gradient.
+ */
+class color : public vector {
+public:
+ static double utop;
+ static double ubottom;
+ static double green;
+public:
+ color(double u) {
+ u = (u - ubottom) / (utop - ubottom);
+ if (u > 1) {
+ u = 1;
+ }
+ if (u < 0) {
+ u = 0;
+ }
+ u = pow(u,2);
+ (*this)[0] = u;
+ (*this)[1] = green * u * (1 - u);
+ (*this)[2] = 1-u;
+ double l = l0norm();
+ for (int i = 0; i < 3; i++) {
+ (*this)[i] /= l;
+ }
+ }
+};
+
+double color::utop = 12;
+double color::ubottom = -31;
+double color::green = 0.5;
+
+/**
+ * \brief Surface model
+ *
+ * This class contains the definitions of the functions to plot
+ * and the parameters to
+ */
+class surfacefunction {
+ static vector axes[6];
+
+ double _a;
+ double _A;
+
+ double _umin;
+ double _umax;
+public:
+ double a() const { return _a; }
+ double A() const { return _A; }
+
+ double umin() const { return _umin; }
+ double umax() const { return _umax; }
+
+ surfacefunction(double a, double A) : _a(a), _A(A), _umin(0), _umax(0) {
+ }
+
+ double f(double z) {
+ return A() * exp(a() * (sqr(z) - 1));
+ }
+
+ double g(double z) {
+ return -f(z) * 2*a() * ((2*a()*sqr(z) + (3-2*a()))*sqr(z) - 1);
+ }
+
+ double F(const vector& v) {
+ double s = 0;
+ for (int i = 0; i < 6; i++) {
+ s += f(axes[i] * v);
+ }
+ return s / 6;
+ }
+
+ double G(const vector& v) {
+ double s = 0;
+ for (int i = 0; i < 6; i++) {
+ s += g(axes[i] * v);
+ }
+ return s / 6;
+ }
+protected:
+ color farbe(const vector& v) {
+ double u = G(v);
+ if (u < _umin) {
+ _umin = u;
+ }
+ if (u > _umax) {
+ _umax = u;
+ }
+ return color(u);
+ }
+};
+
+static double phi = (1 + sqrt(5)) / 2;
+static double sl = sqrt(sqr(phi) + 1);
+vector surfacefunction::axes[6] = {
+ vector( 0. , -1./sl, phi/sl ),
+ vector( 0. , 1./sl, phi/sl ),
+ vector( 1./sl, phi/sl, 0. ),
+ vector( -1./sl, phi/sl, 0. ),
+ vector( phi/sl, 0. , 1./sl ),
+ vector( -phi/sl, 0. , 1./sl )
+};
+
+/**
+ * \brief Class to construct the plot
+ */
+class surface : public surfacefunction {
+ FILE *outfile;
+
+ int _phisteps;
+ int _thetasteps;
+ double _hphi;
+ double _htheta;
+public:
+ int phisteps() const { return _phisteps; }
+ int thetasteps() const { return _thetasteps; }
+ double hphi() const { return _hphi; }
+ double htheta() const { return _htheta; }
+ void phisteps(int s) { _phisteps = s; _hphi = 2 * M_PI / s; }
+ void thetasteps(int s) { _thetasteps = s; _htheta = M_PI / s; }
+
+ surface(const std::string& filename, double a, double A)
+ : surfacefunction(a, A) {
+ outfile = fopen(filename.c_str(), "w");
+ phisteps(400);
+ thetasteps(200);
+ }
+
+ ~surface() {
+ fclose(outfile);
+ }
+
+private:
+ void triangle(const vector& v0, const vector& v1, const vector& v2) {
+ fprintf(outfile, " mesh {\n");
+ vector c = (v0 + v1 + v2) * (1./3.);
+ vector color = farbe(c.normalize());
+ vector V0 = v0 * (1 + F(v0));
+ vector V1 = v1 * (1 + F(v1));
+ vector V2 = v2 * (1 + F(v2));
+ fprintf(outfile, "\ttriangle {\n");
+ fprintf(outfile, "\t <%.6f,%.6f,%.6f>,\n",
+ V0[0], V0[2], V0[1]);
+ fprintf(outfile, "\t <%.6f,%.6f,%.6f>,\n",
+ V1[0], V1[2], V1[1]);
+ fprintf(outfile, "\t <%.6f,%.6f,%.6f>\n",
+ V2[0], V2[2], V2[1]);
+ fprintf(outfile, "\t}\n");
+ fprintf(outfile, "\tpigment { color rgb<%.4f,%.4f,%.4f> }\n",
+ color[0], color[1], color[2]);
+ fprintf(outfile, "\tfinish { metallic specular 0.5 }\n");
+ fprintf(outfile, " }\n");
+ }
+
+ void northcap() {
+ vector v0(0, 0, 1);
+ for (int i = 1; i <= phisteps(); i++) {
+ fprintf(outfile, " // northcap i = %d\n", i);
+ vector v1(htheta(), (i - 1) * hphi());
+ vector v2(htheta(), i * hphi());
+ triangle(v0, v1, v2);
+ }
+ }
+
+ void southcap() {
+ vector v0(0, 0, -1);
+ for (int i = 1; i <= phisteps(); i++) {
+ fprintf(outfile, " // southcap i = %d\n", i);
+ vector v1(M_PI - htheta(), (i - 1) * hphi());
+ vector v2(M_PI - htheta(), i * hphi());
+ triangle(v0, v1, v2);
+ }
+ }
+
+ void zone() {
+ for (int j = 1; j < thetasteps() - 1; j++) {
+ for (int i = 1; i <= phisteps(); i++) {
+ fprintf(outfile, " // zone j = %d, i = %d\n",
+ j, i);
+ vector v0( j * htheta(), (i-1) * hphi());
+ vector v1((j+1) * htheta(), (i-1) * hphi());
+ vector v2( j * htheta(), i * hphi());
+ vector v3((j+1) * htheta(), i * hphi());
+ triangle(v0, v1, v2);
+ triangle(v1, v2, v3);
+ }
+ }
+ }
+public:
+ void draw() {
+ northcap();
+ southcap();
+ zone();
+ }
+};
+
+/**
+ * \brief main function
+ */
+int main(int argc, char *argv[]) {
+ surface S("spherecurve.inc", 5, 10);
+ color::green = 1.0;
+ S.draw();
+ std::cout << "umin: " << S.umin() << std::endl;
+ std::cout << "umax: " << S.umax() << std::endl;
+ return EXIT_SUCCESS;
+}
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-18 15:52:37 +0200