add image code

4 files changed, 153 insertions, 22 deletions

diff --git a/buch/papers/kugel/images/Makefile b/buch/papers/kugel/images/Makefile
index 8efa228..e8bf919 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.png
+all: curvature.jpg
 
 curvature.inc:	curvgraph.m
 	octave curvgraph.m
@@ -11,3 +11,5 @@ curvature.inc:	curvgraph.m
 curvature.png:	curvature.pov curvature.inc
 	povray +A0.1 +W1920 +H1080 +Ocurvature.png curvature.pov
 
+curvature.jpg:	curvature.png
+	convert curvature.png -density 300 -units PixelsPerInch curvature.jpg
diff --git a/buch/papers/kugel/images/curvature.maxima b/buch/papers/kugel/images/curvature.maxima
new file mode 100644
index 0000000..6313642
--- /dev/null
+++ b/buch/papers/kugel/images/curvature.maxima
@@ -0,0 +1,6 @@
+
+f: exp(-r^2/sigma^2)/sigma;
+laplacef: ratsimp(diff(r * diff(f,r), r) / r);
+f: exp(-r^2/(2*sigma^2))/(sqrt(2)*sigma);
+laplacef: ratsimp(diff(r * diff(f,r), r) / r);
+
diff --git a/buch/papers/kugel/images/curvature.pov b/buch/papers/kugel/images/curvature.pov
index 3535488..9dbaa86 100644
--- a/buch/papers/kugel/images/curvature.pov
+++ b/buch/papers/kugel/images/curvature.pov
@@ -11,17 +11,17 @@ global_settings {
 	assumed_gamma 1
 }
 
-#declare imagescale = 0.1;
+#declare imagescale = 0.09;
 
 camera {
-	location <40, 10, -20>
+	location <10, 10, -40>
 	look_at <0, 0, 0>
 	right 16/9 * x * imagescale
 	up y * imagescale
 }
 
 light_source {
-	<10, 10, -40> color White
+	<-10, 10, -40> color White
 	area_light <1,0,0> <0,0,1>, 10, 10
 	adaptive 1
 	jitter
@@ -70,3 +70,69 @@ arrow(<0,-1,0>, <0,1,0>, 0.01, White)
 arrow(<0,0,-2.1>, <0,0,2.1>, 0.01, White)
 
 #include "curvature.inc"
+
+#declare sigma = 1;
+#declare N0 = 0.5;
+#declare funktion = function(r) {
+	(exp(-r*r/(sigma*sigma)) / sigma
+	-
+	exp(-r*r/(2*sigma*sigma)) / (sqrt(2)*sigma)) / N0
+};
+#declare hypot = function(xx, yy) { sqrt(xx*xx+yy*yy) };
+
+#declare Funktion = function(x,y) { funktion(hypot(x+1,y)) - funktion(hypot(x-1,y)) };
+#macro punkt(xx,yy)
+	<xx, Funktion(xx, yy), yy>
+#end
+
+#declare griddiameter = 0.006;
+union {
+	#declare xmin = -3;
+	#declare xmax = 3;
+	#declare ymin = -2;
+	#declare ymax = 2;
+
+
+	#declare xstep = 0.2;
+	#declare ystep = 0.02;
+	#declare xx = xmin;
+	#while (xx < xmax + xstep/2)
+		#declare yy = ymin;
+		#declare P = punkt(xx, yy);
+		#while (yy < ymax - ystep/2)
+			#declare yy = yy + ystep;
+			#declare Q = punkt(xx, yy);
+			sphere { P, griddiameter }
+			cylinder { P, Q, griddiameter }
+			#declare P = Q;
+		#end
+		sphere { P, griddiameter }
+		#declare xx = xx + xstep;
+	#end
+
+	#declare xstep = 0.02;
+	#declare ystep = 0.2;
+	#declare yy = ymin;
+	#while (yy < ymax + ystep/2)
+		#declare xx = xmin;
+		#declare P = punkt(xx, yy);
+		#while (xx < xmax - xstep/2)
+			#declare xx = xx + xstep;
+			#declare Q = punkt(xx, yy);
+			sphere { P, griddiameter }
+			cylinder { P, Q, griddiameter }
+			#declare P = Q;
+		#end
+		sphere { P, griddiameter }
+		#declare yy = yy + ystep;
+	#end
+
+	pigment {
+		color rgb<0.8,0.8,0.8>
+	}
+	finish {
+		metallic
+		specular 0.8
+	}
+}
+
diff --git a/buch/papers/kugel/images/curvgraph.m b/buch/papers/kugel/images/curvgraph.m
index 96ca4b1..b83c877 100644
--- a/buch/papers/kugel/images/curvgraph.m
+++ b/buch/papers/kugel/images/curvgraph.m
@@ -13,23 +13,34 @@ sigma2 = 1;
 global s;
 s = 1;
 
+global cmax;
+cmax = 0.9;
+global cmin;
+cmin = -0.9;
+
+global Cmax;
+global Cmin;
+Cmax = 0;
+Cmin = 0;
+
 xmin = -3;
 xmax = 3;
-xsteps = 1000;
+xsteps = 200;
 hx = (xmax - xmin) / xsteps;
 
 ymin = -2;
 ymax = 2;
-ysteps = 1000;
+ysteps = 200;
 hy = (ymax - ymin) / ysteps;
 
 function retval = f0(r)
 	global sigma2;
-	retval = exp(-r^2/sigma2)/sigma2 - exp(-r^2/(2*sigma2))/(sqrt(2)*sigma2);
+	retval = exp(-r^2/sigma2)/sqrt(sigma2) - exp(-r^2/(2*sigma2))/(sqrt(2*sigma2));
 end
 
 global N0;
-N0 = f0(0);
+N0 = f0(0)
+N0 = 0.5;
 
 function retval = f1(x,y)
 	global N0;
@@ -44,10 +55,10 @@ endfunction
 function retval = curvature0(r)
 	global sigma2;
 	retval = (
-		(2*sigma2-r^2)*exp(-r^2/(2*sigma2))
+		-4*(sigma2-r^2)*exp(-r^2/sigma2)
 		+
-		4*(r^2-sigma2)*exp(-r^2/sigma2)
-	) / (sigma2^2);
+		(2*sigma2-r^2)*exp(-r^2/(2*sigma2))
+	) / (sigma2^(5/2));
 endfunction
 
 function retval = curvature1(x, y)
@@ -56,28 +67,74 @@ endfunction
 
 function retval = curvature(x, y)
 	global s;
-	retval = curvature1(x+s, y) + curvature1(x-s, y);
+	retval = curvature1(x+s, y) - curvature1(x-s, y);
 endfunction
 
 function retval = farbe(x, y)
+	global Cmax;
+	global Cmin;
+	global cmax;
+	global cmin;
 	c = curvature(x, y);
-	retval = c * ones(1,3);
+	if (c < Cmin) 
+		Cmin = c
+	endif
+	if (c > Cmax) 
+		Cmax = c
+	endif
+	u = (c - cmin) / (cmax - cmin);
+	if (u > 1)
+		u = 1;
+	endif
+	if (u < 0)
+		u = 0;
+	endif
+	color = [ u, 0.5, 1-u ];
+	color = color/max(color);
+	color(1,4) = c/2;
+	retval = color;
 endfunction
 
-fn = fopen("curvature.inc", "w");
+function dreieck(fn, A, B, C)
+	fprintf(fn, "\ttriangle {\n");
+	fprintf(fn, "\t    <%.4f,%.4f,%.4f>,\n", A(1,1), A(1,3), A(1,2));
+	fprintf(fn, "\t    <%.4f,%.4f,%.4f>,\n", B(1,1), B(1,3), B(1,2));
+	fprintf(fn, "\t    <%.4f,%.4f,%.4f>\n",  C(1,1), C(1,3), C(1,2));
+	fprintf(fn, "\t}\n");
+endfunction
 
+function viereck(fn, punkte)
+	color = farbe(mean(punkte(:,1)), mean(punkte(:,2)));
+	fprintf(fn, "    mesh {\n");
+	dreieck(fn, punkte(1,:), punkte(2,:), punkte(3,:));
+	dreieck(fn, punkte(2,:), punkte(3,:), punkte(4,:));
+	fprintf(fn, "\tpigment { color rgb<%.4f,%.4f,%.4f> } // %.4f\n",
+		color(1,1), color(1,2), color(1,3), color(1,4));
+	fprintf(fn, "    }\n");
+endfunction
+
+fn = fopen("curvature.inc", "w");
+punkte = zeros(4,3);
 for ix = (0:xsteps)
 	x = xmin + ix * hx;
+	punkte(1,1) = x;      
+	punkte(2,1) = x;
+	punkte(3,1) = x + hx;
+	punkte(4,1) = x + hx;
 	for iy = (0:ysteps)
 		y = ymin + iy * hy;
-		fprintf(fn, "sphere { <%.4f, %.4f, %.4f>, 0.01\n",
-			x, f(x, y), y);
-		color = farbe(x, y);
-		fprintf(fn, "pigment { color rgb<%.4f,%.4f,%.4f> }\n",
-			color(1,1), color(1,2), color(1,3));
-		fprintf(fn, "finish { metallic specular 0.5 }\n");
-		fprintf(fn, "}\n");
+		punkte(1,2) = y;
+		punkte(2,2) = y + hy;
+		punkte(3,2) = y;
+		punkte(4,2) = y + hy;
+		for i = (1:4)
+			punkte(i,3) = f(punkte(i,1), punkte(i,2));
+		endfor
+		viereck(fn, punkte);
 	end
 end
-
+#fprintf(fn, "    finish { metallic specular 0.5 }\n");
 fclose(fn);
+
+printf("Cmax = %.4f\n", Cmax);
+printf("Cmin = %.4f\n", Cmin);