1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

#include <iostream>
#include <cmath>
using namespace std;
//======================================================================
struct Pt { double x, y, z; };
Pt operator + ( Pt p, Pt q ) { return { p.x + q.x, p.y + q.y, p.z + q.z }; }
Pt operator  ( Pt p, Pt q ) { return { p.x  q.x, p.y  q.y, p.z  q.z }; }
Pt operator * ( double r, Pt p ) { return { r * p.x, r * p.y, r * p.z }; }
Pt operator / ( Pt p, double r ) { return { p.x / r, p.y / r, p.z / r }; }
double dot( Pt p, Pt q ) { return p.x * q.x + p.y * q.y + p.z * q.z; }
ostream &operator << ( ostream &out, const Pt &p ) { return out << p.x << " " << p.y << " " << p.z; }
istream &operator >> ( istream &in , Pt &p ) { return in >> p.x >> p.y >> p.z; }
//======================================================================
int main()
{
Pt centre, origin, direction;
double radius;
cout << "Enter sphere centre (vector): "; cin >> centre;
cout << "Enter sphere radius: " ; cin >> radius;
cout << "Enter ray origin (vector): " ; cin >> origin;
cout << "Enter ray direction (vector): "; cin >> direction;
// t satisfies a quadratic
double a = dot( direction, direction );
double b = 2 * dot( direction, origin  centre );
double c = dot( origin  centre, origin  centre )  radius * radius;
double discriminant = b * b  4 * a * c;
if ( discriminant < 0 )
{
cout << "No intersection.\n";
}
else if ( discriminant > 0 )
{
double t = ( b + sqrt( discriminant ) ) / ( 2 * a );
double t2 = b / a  t;
if ( abs( t2 ) < abs( t ) ) t = t2;
cout << "Closest intersection at " << origin + t * direction << '\n';
}
else
{
cout << "Touches sphere at " << origin + (0.5 * b / a) * direction << '\n';
}
}
 