У меня есть найти решение: здесь код в C#:
Double alpha = Math.Asin(sphere.Radius/(Math.Sqrt(Math.Pow(sphere.Origin.X-ray.Origin.X,2)+Math.Pow(sphere.Origin.Y-ray.Origin.Y,2)+Math.Pow(sphere.Origin.Z-ray.Origin.Z,2))));
Double beta = Math.Acos((ray.Direction * (sphere.Origin - ray.Origin))/(ray.Direction.Length * (sphere.Origin - ray.Origin).Length));
ray.HitParam = VypA(sphere.Origin - ray.Origin, beta, sphere.Radius)/ray.Direction.Length;
Vector4 g = ray.Origin + ray.Direction * ray.HitParam;
ray.HitNormal = (g - sphere.Origin).Normalized;
////the VypA function
public static Double VypA(Vector4 b, Double beta, Double radius)
{
return b.Length * (Math.Cos(beta) - Math.Sqrt(((radius * radius)/(b.Length2) - (Math.Sin(beta) * (Math.Sin(beta))))));
}