From 98b6050e8369bfac26c57d7e6dee847ba2e0c074 Mon Sep 17 00:00:00 2001
From: Werner Almesberger <werner@openmoko.org>
Date: Sat, 25 Sep 2010 03:15:11 -0300
Subject: [PATCH] Fix our trigonometry: the z0 inclination is expressed as a
 tangent, not a sine.

- solidify/project.c (read_face_data, save_face_data): fx/fy represent a
  tangent, not a sine
- solidify/overlap.c (merge_matrix): corrected shrinkage formula to use
  tangent instead of sine
---
 solidify/overlap.c | 15 +++++++++++----
 solidify/project.c |  6 +++---
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/solidify/overlap.c b/solidify/overlap.c
index e47e13a..4ecec71 100644
--- a/solidify/overlap.c
+++ b/solidify/overlap.c
@@ -169,12 +169,19 @@ static void merge_matrix(struct matrix *m, const struct solid *s,
 	m->b[1] += f->cy;
 
 	/*
-	 * Apply shrinkage caused by rotation out of z0. We use that
-	 * cos a = sqrt(1-sin^2 a)
+	 * Apply shrinkage caused by rotation out of z0.
+	 * We need to divide by x = cos a. We have f = tan a.
+	 * With sin^2 a + cos^2 a = 1, we get
+	 *
+	 * f = sqrt(1-cos^2 a)/cos a
+	 *   = sqrt(1-x^2)/x
+	 * f^2 = 1/x^2-1
+	 * 1/(f^2+1) = x^2
+	 * cos a = sqrt(1/(f^2+1))
 	 */
 
-	f_x = 1.0/sqrt(1-f->fx*f->fx);
-	f_y = 1.0/sqrt(1-f->fy*f->fy);
+	f_x = sqrt(f->fx*f->fx+1);
+	f_y = sqrt(f->fy*f->fy+1);
 
 	m->a[0][0] *= f_x;
 	m->a[0][1] *= f_x;
diff --git a/solidify/project.c b/solidify/project.c
index 2f803b8..792d3fc 100644
--- a/solidify/project.c
+++ b/solidify/project.c
@@ -71,11 +71,11 @@ static void read_face_data(FILE *file, struct face *f)
 
 	if (fscanf(file, "%f", &v) != 1)
 		return;
-	f->fx = sin(v/180*M_PI);
+	f->fx = tan(v/180*M_PI);
 
 	if (fscanf(file, "%f", &v) != 1)
 		return;
-	f->fy = sin(v/180*M_PI);
+	f->fy = tan(v/180*M_PI);
 
 	if (fscanf(file, "%f", &v) != 1)
 		return;
@@ -152,7 +152,7 @@ static void save_face_data(FILE *file, const char *name, const struct face *f)
 		a = 180-a;
 	if (fprintf(file, "%g %g %g\n%g %g %g\n",
 	    f->z_ref*f->z_step,
-	    asin(f->fx)/M_PI*180, asin(f->fy)/M_PI*180,
+	    atan(f->fx)/M_PI*180, atan(f->fy)/M_PI*180,
 	    a, f->m.b[0]*f->x_step, f->m.b[1]*f->y_step) < 0) {
 		perror(name);
 		exit(1);