Merge pull request #5 from AndreasFMueller/master

update

5 files changed, 554 insertions, 1 deletions

diff --git a/buch/papers/kugel/images/Makefile b/buch/papers/kugel/images/Makefile
index e8bf919..4226dab 100644
--- a/buch/papers/kugel/images/Makefile
+++ b/buch/papers/kugel/images/Makefile
@@ -3,7 +3,7 @@
 #
 # (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
 #
-all: curvature.jpg
+all: curvature.jpg spherecurve.jpg
 
 curvature.inc:	curvgraph.m
 	octave curvgraph.m
@@ -13,3 +13,18 @@ curvature.png:	curvature.pov curvature.inc
 
 curvature.jpg:	curvature.png
 	convert curvature.png -density 300 -units PixelsPerInch curvature.jpg
+
+spherecurve2.inc:	spherecurve.m
+	octave spherecurve.m
+
+spherecurve.png:	spherecurve.pov spherecurve.inc
+	povray +A0.1 +W1080 +H1080 +Ospherecurve.png spherecurve.pov
+
+spherecurve.jpg:	spherecurve.png
+	convert spherecurve.png -density 300 -units PixelsPerInch spherecurve.jpg
+
+spherecurve:	spherecurve.cpp
+	g++ -o spherecurve -g -Wall -O spherecurve.cpp
+
+spherecurve.inc:	spherecurve
+	./spherecurve
diff --git a/buch/papers/kugel/images/spherecurve.cpp b/buch/papers/kugel/images/spherecurve.cpp
new file mode 100644
index 0000000..8ddf5e5
--- /dev/null
+++ b/buch/papers/kugel/images/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;
+}
diff --git a/buch/papers/kugel/images/spherecurve.m b/buch/papers/kugel/images/spherecurve.m
new file mode 100644
index 0000000..99d5c9a
--- /dev/null
+++ b/buch/papers/kugel/images/spherecurve.m
@@ -0,0 +1,160 @@
+#
+# spherecurve.m
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+global a;
+a = 5;
+global A;
+A = 10;
+
+phisteps = 400;
+hphi = 2 * pi / phisteps;
+thetasteps = 200;
+htheta = pi / thetasteps;
+
+function retval = f(z)
+	global a;
+	global A;
+	retval = A * exp(a * (z^2 - 1));
+endfunction
+
+function retval = g(z)
+	global a;
+	retval = -f(z) * 2 * a * (2 * a * z^4 + (3 - 2*a) * z^2 - 1);
+	#                                                              2
+	#                   - a     2  4             2   2          a z
+	#(%o6)          - %e    (4 a  z  + (6 a - 4 a ) z  - 2 a) %e
+endfunction
+
+phi = (1 + sqrt(5)) / 2;
+
+global axes;
+axes = [
+     0,   0,   1,  -1, phi, -phi;
+     1,  -1, phi, phi,   0,    0;
+   phi, phi,   0,   0,   1,    1;
+];
+axes = axes / (sqrt(phi^2+1));
+
+function retval = kugel(theta, phi)
+	retval = [
+	    cos(phi) * sin(theta);
+	    sin(phi) * sin(theta);
+	               cos(theta)
+	];
+endfunction
+
+function retval = F(v)
+	global axes;
+	s = 0;
+	for i = (1:6)
+		z = axes(:,i)' * v;
+		s = s + f(z);
+	endfor
+	retval = s / 6;
+endfunction
+
+function retval = F2(theta, phi)
+	v = kugel(theta, phi);
+	retval = F(v);
+endfunction
+
+function retval = G(v)
+	global axes;
+	s = 0;
+	for i = (1:6)
+		s = s + g(axes(:,i)' * v);
+	endfor
+	retval = s / 6;
+endfunction
+
+function retval = G2(theta, phi)
+	v = kugel(theta, phi);
+	retval = G(v);
+endfunction
+
+function retval = cnormalize(u)
+	utop = 11;
+	ubottom = -30;
+	retval = (u - ubottom)  / (utop - ubottom);
+	if (retval > 1) 
+		retval = 1;
+	endif
+	if (retval < 0)
+		retval = 0;
+	endif
+endfunction
+
+global umin;
+umin = 0;
+global umax;
+umax = 0;
+
+function color = farbe(v)
+	global umin;
+	global umax;
+	u = G(v);
+	if (u < umin) 
+		umin = u;
+	endif
+	if (u > umax)
+		umax = u;
+	endif
+	u = cnormalize(u);
+	color = [ u, 0.5, 1-u ];
+	color = color/max(color);
+endfunction
+
+function dreieck(fn, v0, v1, v2)
+	fprintf(fn, "    mesh {\n");
+	c = (v0 + v1 + v2) / 3;
+	c = c / norm(c);
+	color = farbe(c);
+	v0 = v0 * (1 + F(v0));
+	v1 = v1 * (1 + F(v1));
+	v2 = v2 * (1 + F(v2));
+	fprintf(fn, "\ttriangle {\n");
+	fprintf(fn, "\t    <%.6f,%.6f,%.6f>,\n", v0(1,1), v0(3,1), v0(2,1));
+	fprintf(fn, "\t    <%.6f,%.6f,%.6f>,\n", v1(1,1), v1(3,1), v1(2,1));
+	fprintf(fn, "\t    <%.6f,%.6f,%.6f>\n",  v2(1,1), v2(3,1), v2(2,1));
+	fprintf(fn, "\t}\n");
+	fprintf(fn, "\tpigment { color rgb<%.4f,%.4f,%.4f> }\n",
+		color(1,1), color(1,2), color(1,3));
+	fprintf(fn, "\tfinish { metallic specular 0.5 }\n");
+	fprintf(fn, "    }\n");
+endfunction
+
+fn = fopen("spherecurve2.inc", "w");
+
+	for i = (1:phisteps)
+		# Polkappe nord
+		v0 = [ 0; 0; 1 ];
+		v1 = kugel(htheta, (i-1) * hphi);
+		v2 = kugel(htheta,  i    * hphi);
+		fprintf(fn, "    // i = %d\n", i);
+		dreieck(fn, v0, v1, v2);
+
+		# Polkappe sued
+		v0 = [ 0; 0; -1 ];
+		v1 = kugel(pi-htheta, (i-1) * hphi);
+		v2 = kugel(pi-htheta,  i    * hphi);
+		dreieck(fn, v0, v1, v2);
+	endfor
+
+	for j = (1:thetasteps-2)
+		for i = (1:phisteps)
+			v0 = kugel( j    * htheta, (i-1) * hphi);
+			v1 = kugel((j+1) * htheta, (i-1) * hphi);
+			v2 = kugel( j    * htheta,  i    * hphi);
+			v3 = kugel((j+1) * htheta,  i    * hphi);
+			fprintf(fn, "    // i = %d, j = %d\n", i, j);
+			dreieck(fn, v0, v1, v2);
+			dreieck(fn, v1, v2, v3);
+		endfor
+	endfor
+
+fclose(fn);
+
+umin
+umax
diff --git a/buch/papers/kugel/images/spherecurve.maxima b/buch/papers/kugel/images/spherecurve.maxima
new file mode 100644
index 0000000..1e9077c
--- /dev/null
+++ b/buch/papers/kugel/images/spherecurve.maxima
@@ -0,0 +1,13 @@
+/*
+ * spherecurv.maxima
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+f: exp(-a * sin(theta)^2);
+
+g: ratsimp(diff(sin(theta) * diff(f, theta), theta)/sin(theta));
+g: subst(z, cos(theta), g);
+g: subst(sqrt(1-z^2), sin(theta), g);
+ratsimp(g);
+
+f: ratsimp(subst(sqrt(1-z^2), sin(theta), f));
diff --git a/buch/papers/kugel/images/spherecurve.pov b/buch/papers/kugel/images/spherecurve.pov
new file mode 100644
index 0000000..b1bf4b8
--- /dev/null
+++ b/buch/papers/kugel/images/spherecurve.pov
@@ -0,0 +1,73 @@
+//
+// curvature.pov
+//
+// (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+//
+
+#version 3.7;
+#include "colors.inc"
+
+global_settings {
+	assumed_gamma 1
+}
+
+#declare imagescale = 0.13;
+
+camera {
+	location <10, 10, -40>
+	look_at <0, 0, 0>
+	right x * imagescale
+	up y * imagescale
+}
+
+light_source {
+	<-10, 10, -40> color White
+	area_light <1,0,0> <0,0,1>, 10, 10
+	adaptive 1
+	jitter
+}
+
+sky_sphere {
+	pigment {
+		color rgb<1,1,1>
+	}
+}
+
+//
+// draw an arrow from <from> to <to> with thickness <arrowthickness> with
+// color <c>
+//
+#macro arrow(from, to, arrowthickness, c)
+#declare arrowdirection = vnormalize(to - from);
+#declare arrowlength = vlength(to - from);
+union {
+	sphere {
+		from, 1.1 * arrowthickness
+	}
+	cylinder {
+		from,
+		from + (arrowlength - 5 * arrowthickness) * arrowdirection,
+		arrowthickness
+	}
+	cone {
+		from + (arrowlength - 5 * arrowthickness) * arrowdirection,
+		2 * arrowthickness,
+		to,
+		0
+	}
+	pigment {
+		color c
+	}
+	finish {
+		specular 0.9
+		metallic
+	}
+}
+#end
+
+arrow(<-2.7,0,0>, <2.7,0,0>, 0.03, White)
+arrow(<0,-2.7,0>, <0,2.7,0>, 0.03, White)
+arrow(<0,0,-2.7>, <0,0,2.7>, 0.03, White)
+
+#include "spherecurve.inc"
+