计算几何与圆心和球有关的计算问题

struct Line
{
	Vector v;
	Point p,b;
	Point point(double t)
	{
		return p+v*t;
	}
	Line(Point a,Point b,Vector v) : p(a),b(b),v(v) {}
	Line(Point p,Vector v) : p(p),v(v){}
};

struct Circle
{
	double r;
	Point  c;
	Circle(Point c = 0,double r = 0):c(c),r(r){};
	Point point(double a)
	{
		return Point(c.x+cos(a)*r,c.y+sin(a)*r);
	}
};


int getLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point> & sol)
{
	double a = L.v.x,b = L.p.x-C.c.x,c = L.v.y,d = L.p.y-C.c.y;
	double e = a*a + c*c, f = 2*a*b + 2*c*d , g = b*b+d*d-C.r*C.r;
	double delta = f*f - 4.0*e*g;
	if(dcmp(delta) < 0) return 0;
	if(dcmp(delta) == 0)
	{
		t1 = t2 = (-f  ) / (2*e);
		sol.push_back(L.point(t1));
		return 1;
	}
	else
	{
		t1 = (-f - sqrt(delta)) / (2*e);
		sol.push_back(L.point(t1));
		t2 = (-f + sqrt(delta)) / (2*e);
		sol.push_back(L.point(t2));
		return 2;
	}
}

int getCircleCircleIntersection(Circle C1,Circle C2,vector<Point> & sol)
{
	double d = Length(C1.c - C2.c);
	if(dcmp(d) == 0)
	{
		if(dcmp(C1.r-C2.r) == 0) return -1;
		return 0;
	}
	if(dcmp(C1.r + C2.r - d) < 0) return 0;
	if(dcmp(fabs(C1.r - C2.r ) - d)   > 0) return 0;

	double a = angle(C2.c-C1.c);
	double da = acos( (C1.r*C1.r + d*d - C2.r*C2.r) / (2 * C1.r*d) );
	Point p1 = C1.point(a - da),p2 = C1.point(a + da);
	sol.push_back(p1);
	if(p1 == p2)
		return 1;
	sol.push_back(p2);
	return 2;
}

int getTangents(Point p,Circle C,Vector * v)
{
	Vector u = C.c - p;
	double dist = Length(u);
	if(dist < C.r) return 0;
	else if(dcmp(dist-C.r) == 0) 
	{
		v[0] = Rotate(u,PI/2);
		return 1;
	}
	else
	{
		double ang = asin(C.r / dist);
		v[0] = Rotate(u,-ang);
		v[1] = Rotate(u,+ang);
		return 2;
	}
}

int getTangents(Circle A,Circle B,Point *a,Point *b)
{
	int cnt = 0;
	if(A.r < B.r) {swap(A,B),swap(a,b);}
	int d2 = (A.c.x-B.c.x)*(A.c.x-B.c.x) + (A.c.y-B.c.y)* (A.c.y-B.c.y);
	int rdiff = A.r - B.r;
	int rsum = A.r + B.r;
	if(d2 < rdiff*rdiff) return 0;

	double base = atan2(B.c.y - A.c.y,B.c.x - A.c.x);
	if(d2 == 0 && A.r == B.r) return -1;
	if(d2 == rdiff*rdiff)
	{
		a[cnt] = A.point(base),b[cnt] = B.point(base); cnt++;
		return 1;
	}

	double ang = acos((A.r-B.r) / sqrt(d2) );
	a[cnt] = A.point(base + ang);b[cnt] = B.point(base + ang);cnt++;
	a[cnt] = A.point(base - ang);b[cnt] = B.point(base - ang);cnt++;
	if(d2 == rsum*rsum)
	{
		a[cnt] = A.point(base);b[cnt] = B.point(base+PI);cnt++;
	}
	else if(d2 > rsum*rsum)
	{
		double ang = acos((A.r+B.r) / sqrt(d2) );
		a[cnt] = A.point(base+ang); b[cnt] = B.point(PI+base+ang);cnt++;
		a[cnt] = A.point(base-ang); b[cnt] = B.point(PI+base-ang);cnt++;
	}
	return cnt;
}

Point GetLineIntersection(Point P,Vector v,Point Q,Vector w)
{
    Vector u = P - Q;
    double t = Cross(w,u)/Cross(v,w);
    return P + v*t;
}

Circle CircumscribedCircle(Point a,Point b,Point c)
{
	Point a2 = Point((a.x+c.x) / 2.0,(a.y+c.y) / 2.0);
	Point b2 = Point((a.x+b.x) / 2.0,(a.y+b.y) / 2.0);
	
	Vector tmp1 = a - a2;
	tmp1 = Rotate(tmp1,PI/2); 

	Vector tmp2 = a - b2;
	tmp2 = Rotate(tmp2,PI/2);

	Point t = GetLineIntersection(a2,tmp1,b2,tmp2);
	double r = Length(a-t);
	return Circle(t,r);
}

Circle InscribedCircle(Point p1,Point p2,Point p3)
{
	double a = Length(p1-p2);
	double b = Length(p2-p3);
	double c = Length(p3-p1);
	Point p = (p3*a+p2*c+p1*b) / (a+b+c);
	return Circle(p,DistanceToLine(p,p1,p2));
}
  
double torad(double deg){ return deg / 180 * PI;}
void get_coord(double R,double lat,double lng double &x,double &y,double&z)
{
	double lat = torad(lat);
	double lng = torad(lng);
	x = R*cos(lat)*sin(lng);
	y = R*cos(lat)*sin(lng);
	z = R*sin(lat);
}