﻿ Proto HTML

Represent the ray as:
x = a + t b
where a, b and x are 3-vectors

x' A x + x' B + C = 0
where A is a symmetric 3x3 matrix, B a 3x1 vector, and C a scalar. Note that x' A x is the dot product of x and A x.

Substituting in gives the quadratic equation:

```              2
l2 * t  + l1 * t + l0 = 0
```
where:
```	l2 = b' A b

l1 = 2 b' A a + b' B

l0 = a' A a + b' B + C
```
If the discriminant of this equation is less than 0, the line does not intersect the quadric. Solving the equation and substituting the values of t into the ray equation will give you the points.

```int Torus::line_intersect(float ax, float ay, float az,
float bx, float by, float bz,
int * num_intersections,
float * intersections,
//
// Intersect ray
//          [ax]     [bx]
//	    [ay] + t [by]
//	    [az]     [bz]
//
{

// This struct is syntactic sugar so that a.b,
// (the dot product) looks just right (:-)
struct { float a,b; } a;
a.b = ax*bx + ay*by + az*bz;
a.a = ax*ax + ay*ay + az*az;

// Set up quartic in t:
//
//  4     3     2
// t + A t + B t + C t + D = 0
//
float R2 = R*R;
float K = a.a - r*r - R2;
float A = 4*a.b;
float B = 2*(2*a.b*a.b + K + 2*R2*bz*bz);
float C = 4*(K*a.b + 2*R2*az*bz);
float D = K*K + 4*R2*(az*az - r*r);

// Solve quartic...
double roots[4];
int nroots = solve_quartic(A,B,C,D,roots,SOLVE_ALL);

*num_intersections = 0;
while(nroots--)
{
float t = roots[nroots];
float x = ax + t*bx;
float y = ay + t*by;
float l = R*(M_PI/2 - atan2(y,x));
if (l <= vlength && l >= 0)
intersections[(*num_intersections)++] = t;
}
}
```