题目大意:在平面上有三个没有公共部分的圆,求平面上一点使得到三个圆的切线的夹角相等。
题解:根据题意易知,要求的这个点和每个圆的圆心以及切点构成的三角形是相似的。因此该点是以三个圆心为圆心,半径之比等于三个圆半径之比的三个圆的交点。可以先求两个圆的交点,再看这两个点和第三个圆的关系,因为这个关系是单调的,所以可以二分三个圆的半径,使得它们恰好交于一点。
求圆的交点可以直接解方程,比较繁琐,这里可以用向量法来解,但是要分几种情况讨论。
设圆心距为
O0O1
情况1
max(r1,r2)≤O0O1≤r1+r2
首先可以求出
∠BO1O0
,令其等于
θ
,可求
cosθ=r21+O0O21−r222r1O0O1
,则
QB=r1sinθ
,
O0Q
和
O1Q
可求,
Q
坐标可由圆心坐标推出,
情况2
这种情况和上面类似,只是向量
O0Q
和
O0O1
反向。
情况3
r1+r2<O0O1
或
O0O1+min(r1,r2)<max(r1,r2)
此时无交点。
精度是深坑…感觉用了一坨数值函数分各种情况讨论还不如解方程…
在判是否在圆上的时候,可以用相对误差来判断,可能会好一些。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const double eps2 = 1e-5;
const double inf = 1e12;
struct Point{
double x,y;
Point(double x,double y):x(x),y(y){}
Point(){}
double length2(){
return x*x+y*y;
}
double length(){
return sqrt(length2());
}
Point operator-(const Point &a)const{
return Point(x-a.x,y-a.y);
}
Point operator+(const Point &a)const{
return Point(x+a.x,y+a.y);
}
Point operator*(const double a)const{
return Point(x*a,y*a);
}
Point unity(){
Point ret;
double l = length();
ret.x =x/l;
ret.y =y/l;
return ret;
}
};
struct Cir{
double x,y,r;
bool operator<(const Cir &a)const{
return r<a.r;
}
bool judge(const Point &a){
double xx,yy;
xx = a.x-x;
yy = a.y-y;
double tmp = xx*xx+yy*yy;
double rr = r*r;
//printf("xx=%f yy=%f tmp=%f rr=%f eps = %f eps2=%f\n",xx,yy,tmp,rr,fabs(tmp-rr),tmp/rr);
if(fabs(tmp/rr-1)<eps2) return true;
return false;
}
bool in(const Point &a){
if((a.x-x)*(a.x-x)+(a.y-y)*(a.y-y)-r*r<0) return true;
return false;
}
};
Cir c[3];
Point Ver(const Point &a){
Point ret;
if(fabs(a.x)<eps){
ret.x = 1;
ret.y = 0;
}
else{
ret.y = 1;
ret.x = -a.y/a.x;
}
ret = ret.unity();
return ret;
}
int main(){
for(int i = 0;i < 3;i++){
scanf("%lf%lf%lf",&c[i].x,&c[i].y,&c[i].r);
}
//for(int i = 0;i < 3;i++){
// printf("O%d=%f,%f\n",i,c[i].x,c[i].y);
//}
sort(c,c+3);
double L,R,mid;
L = c[0].r;
R = 2000;
c[1].r /= c[0].r;
c[2].r /= c[0].r;
c[0].r = 1;
mid = (L+R)/2;
Point OO1 = Point(c[1].x-c[0].x,c[1].y-c[0].y);
Point O0 = Point(c[0].x,c[0].y);
//printf("OO1=%f,%f e(OO1)=%f,%f\n",OO1.x,OO1.y,OO1.unity().x,OO1.unity().y);
Point p1,p2;
//printf("L %f R %f\n",L,R);
int cnt = 1;
Cir cc = c[2];
while(R-L>eps){
double r[3];
for(int i = 0;i < 3;i++) r[i] = c[i].r * mid;
cc.r = r[2];
if(r[0]+r[1]<OO1.length()){
L = mid;
mid = (L+R)/2;
//printf("i=%d L=%f R=%f mid=%f\n\n",cnt++,L,R,mid);
continue;
}
else if(r[0]+OO1.length()<r[1]){
R = mid;
mid = (L+R)/2;
//printf("i=%d L=%f R=%f mid=%f\n\n",cnt++,L,R,mid);
continue;
}
double cossita = (r[1]*r[1]+OO1.length2()-r[0]*r[0])/(2*r[1]*OO1.length());
int flag;
if(OO1.length2()+r[0]*r[0]<r[1]*r[1]) flag = -1;
else flag = 1;
//printf("sin(sita)=%f cos(sita)=%f\n",sin(sita),cos(sita));
Point Q = O0 + OO1.unity()*sqrt(r[0]*r[0]-r[1]*r[1]*(1-cossita*cossita))*flag;
//printf("Q=%f,%f\n",Q.x,Q.y);
Point v = Ver(OO1)*r[1]*sin(acos(cossita));
p1 = Q+v;
p2 = Q-v;
if(cc.in(p1)||cc.in(p2)){
R = mid;
mid = (L+R)/2;
//printf("i=%d L=%f R=%f mid=%f p1=%f,%f p2=%f,%f\n\n",cnt++,L,R,mid,p1.x,p1.y,p2.x,p2.y);
}
else{
L = mid;
mid = (L+R)/2;
//printf("i=%d L=%f R=%f mid=%f p1=%f,%f p2=%f,%f\n\n",cnt++,L,R,mid,p1.x,p1.y,p2.x,p2.y);
}
}
//printf("mid %f\n",mid);
//printf("cc=%f,%f,%f\n",cc.x,cc.y,cc.r);
//printf("p1=%f,%f p2=%f,%f\n",p1.x,p1.y,p2.x,p2.y);
if(cc.judge(p1)){
printf("%.5f %.5f\n",p1.x,p1.y);
return 0;
}
else if(cc.judge(p2)){
printf("%.5f %.5f\n",p2.x,p2.y);
return 0;
}
return 0;
}