summaryrefslogtreecommitdiffstats
path: root/libdimension/cylinder.c
diff options
context:
space:
mode:
authorTavian Barnes <tavianator@gmail.com>2010-10-26 15:42:13 -0400
committerTavian Barnes <tavianator@gmail.com>2010-10-26 15:42:13 -0400
commit0075e0c37d9c33ed00e4308e6444b61b204327ba (patch)
tree4c112619d917fc8374715bde8ce4e4f0e6723892 /libdimension/cylinder.c
parent1fba91c6fe9115be67929ce1e247dd759a21fcd1 (diff)
downloaddimension-0075e0c37d9c33ed00e4308e6444b61b204327ba.tar.xz
Add numerical polynomial solver based on Uspensky's algorithm.
Diffstat (limited to 'libdimension/cylinder.c')
-rw-r--r--libdimension/cylinder.c34
1 files changed, 21 insertions, 13 deletions
diff --git a/libdimension/cylinder.c b/libdimension/cylinder.c
index e9db50c..ba016cf 100644
--- a/libdimension/cylinder.c
+++ b/libdimension/cylinder.c
@@ -69,19 +69,27 @@ dmnsn_cylinder_intersection_fn(const dmnsn_object *cylinder, dmnsn_line line,
/* Solve (x0 + nx*t)^2 + (z0 + nz*t)^2
== (((r2 - r1)*(y0 + ny*t) + r1 + r2)/2)^2 */
- double a, b, c, t;
- a = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
- b = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
- - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
- c = l.x0.x*l.x0.x + l.x0.z*l.x0.z
- - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
-
- if (b*b - 4.0*a*c >= 0.0) {
- t = (-b - sqrt(b*b - 4.0*a*c))/(2.0*a);
- dmnsn_vector p = dmnsn_line_point(l, t);
-
- if (t < 0.0 || p.y <= -1.0 || p.y >= 1.0) {
- t = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a);
+ double poly[3], x[2];
+ poly[2] = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
+ poly[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
+ - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
+ poly[0] = l.x0.x*l.x0.x + l.x0.z*l.x0.z
+ - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4;
+
+ size_t n = dmnsn_solve_polynomial(poly, 2, x);
+
+ if (n > 0) {
+ double t = x[0];
+ dmnsn_vector p;
+ if (n == 2) {
+ t = dmnsn_min(t, x[1]);
+ p = dmnsn_line_point(l, t);
+
+ if (p.y <= -1.0 || p.y >= 1.0) {
+ t = dmnsn_max(x[0], x[1]);
+ p = dmnsn_line_point(l, t);
+ }
+ } else {
p = dmnsn_line_point(l, t);
}