计算分子表面的时候，需要计算直线和球面的交点

https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection

直线和球面的交点：

  ///begin with a segment from r0 to r1
   ///
   ///The points p on the line are described as p = r0 + u*(r1-r0), or, more specifically,
   ///
   ///p[0] = r0[0] + u*(r1[0]-r0[0])
   ///p[1] = r0[1] + u*(r1[1]-r0[1])
   ///p[2] = r0[2] + u*(r1[2]-r0[2])
   ///
   ///The sphere centered at c with radius r is described as
   ///
   /// (x-c[0])^2 + (y-c[1])^2 + (z - c[2])^2 = r^2
   ///
   ///If x,y,z is the point of intersection with the line, then x=p[0], y=p[1], z=p[2], and
   ///we can substitute the expressions above into the sphere equation.  If we evaluate the
   ///squared quantities ( (x-c[0])^2, etc ) and separate out the terms into terms with
   ///u^2, terms with u, and constant terms, we get a quadratic equation of the form
   ///
   /// a*(u^2) + b*u + c = 0
   ///
   ///where
   ///
   /// a = (r1[0]-r0[0])^2 + (r1[1]-r0[1])^2 + (r1[2]-r0[2])^2
   /// b = 2*( ((r1[0]-r0[0])*(r[0]-c[0])) + ((r1[1]-r0[1])*(r[1]-c[1])) + ((r1[2]-r0[2])*(r[2]-c[2])) )
   /// c = c[0]^2 + c[1]^2 + c[2]^2 + r0[0]^2 + r0[1]^2 + r0[2]^2 + 2( (c[0]*r0[0]) + (c[1]*r0[1]) + (c[2]*r0[2]) ) - r^2
   ///
   ///  Since this is quadratic, its solutions are
   ///
   ///       -b +-  sqrt( b^2 - 4*a*c )
   ///  u = -----------------------------
   ///                    2a
   ///
   ///  Which has a descriminant (b^2 - 4*a*c).  
   ///  CASE 1: If the discriminant is less than zero, the segment does not intersect the sphere
   ///  CASE 2: If the discriminant is equal to zero, then the segment intersects the sphere at u = (-b)/(2*a)
   ///  CASE 3: If the discriminant is greater than zero, then the line intersects the sphere at two places.
   
   ///begin by computing the components of the quadratic:

--------------------------------------

Vector3 intersect_Seg_Sphere(Vector3 pt1,Vector3 pt2,Vector3 origin, float r)
{
   Vector3 pt;
   float dist1=pt1.distance(origin);
   float dist2=pt2.distance(origin);
   if( ((dist1 >r)&& (dist2 >r))  || ((dist1 <r) &&(dist2<r) ))
   {
       return pt;
   }
   //https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
   float quad_a,quad_b,quad_c;
   quad_a=(pt2.distance(pt1))*(pt2.distance(pt1));
   quad_b=2*(
       ((pt2.xvalue()-pt1.xvalue())*(pt1.xvalue()-origin.xvalue()))+
       ((pt2.yvalue()-pt1.yvalue())*((pt1.yvalue()-origin.yvalue())))+
       ((pt2.zvalue()-pt1.zvalue())*(pt1.zvalue()-origin.zvalue())));

   float disc=(quad_b*quad_b)-4*quad_a*quad_c;

   if(disc <-.00001){return pt;};
   if((disc<0.00001) &&(disc > -.00001))
   {
       float u=(-quad_b)/(2*quad_a);
       pt.setx(pt1.xvalue()+u*(pt2.xvalue()-pt1.xvalue()));
       pt.sety(pt1.yvalue()+u*(pt2.yvalue()-pt1.yvalue()));
       pt.setz(pt1.zvalue()+u*(pt2.zvalue()-pt1.zvalue()));
       return pt;
   }

   if(disc > 0.00001)
   {
       float u1 = (-quad_b - sqrt( (quad_b*quad_b) - (4*quad_a*quad_c) ) )/(2*quad_a);
       float u2 = (-quad_b + sqrt( (quad_b*quad_b) - (4*quad_a*quad_c) ) )/(2*quad_a);
       //intersect point
       //pt=pt1+u(pt2-pt1)
       //intersect point should local between the pt1 and pt2
       double u;
       if( u1<0 || u1>1 ){ u = u2; }
       if( u2<0 || u2>1 ){ u = u1; }

       pt.setx(pt1.xvalue()+u*(pt2.xvalue()-pt1.xvalue()));
       pt.sety(pt1.yvalue()+u*(pt2.yvalue()-pt1.yvalue()));
       pt.setz(pt1.zvalue()+u*(pt2.zvalue()-pt1.zvalue()));

       return pt;
   }

}

--------------------------------------