推荐理由:
细节较多的计算几何,适合一做。
Describe
在平面上有三个没有公共部分的圆,求平面上一点使得到三个圆的切线的夹角相等。
Solution
容易得出:
稍加转换:
那么我们只要找到一个点 ,满足他到
个点的距离符合给定比例。
我们设出
平方一下,然后乘以分母。可以得到
我们把 看做常数,那么原方程可以化成一个
有点常识的人就可以知道这是圆的一般方程,然后再分 种情况讨论。
那么 的轨迹是一条直线。
只需要把三个圆两两向交,再分开判断就行了。
code
#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;
const bool LINE = true;
const bool CIRCLE = false;
const double EPS = 1e-8;
struct Point { double x, y;};
struct Node {
bool type; double r;
Point O, p[2];
} node[3];
int x[3], y[3], r[3];
double dis(const Point & p1, const Point & p2) {
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
void getget(int i, int j, int k) {
if(r[i] == r[j]) { node[k].type = LINE;
node[k].p[0].x = (x[i] + x[j]) / 2.0;
node[k].p[0].y = (y[i] + y[j]) / 2.0;
double dx = x[j] - x[i], dy = y[j] - y[i];
double temp = dx; dx = -dy, dy = temp;
node[k].p[1].x = node[k].p[0].x + dx;
node[k].p[1].y = node[k].p[0].y + dy;
} else { node[k].type = CIRCLE;
if(r[i] > r[j]) swap(i, j);
double dx = x[j] - x[i], dy = y[j] - y[i];
double d = sqrt(dx * dx + dy * dy);
node[k].p[0].x = x[i] + dx * r[i] / (r[i] + r[j]);
node[k].p[0].y = y[i] + dy * r[i] / (r[i] + r[j]);
node[k].p[1].x = x[i] - dx * r[i] / (r[j] - r[i]);
node[k].p[1].y = y[i] - dy * r[i] / (r[j] - r[i]);
node[k].r = dis(node[k].p[0], node[k].p[1]) / 2.0;
node[k].O.x = (node[k].p[0].x + node[k].p[1].x) / 2.0;
node[k].O.y = (node[k].p[0].y + node[k].p[1].y) / 2.0;
}
} bool zero(double x) { return ((x > 0 ? x : -x) < EPS);}
bool work(const Point & u1, const Point & u2, const Point & v1, const Point & v2) {
return zero((u1.x - u2.x) * (v1.y - v2.y) - (v1.x - v2.x) * (u1.y - u2.y));
}
Point solve(const Point & u1, const Point & u2, const Point & v1, const Point & v2) { Point ret = u1;
double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y - v1.y) * (v1.x - v2.x)) / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x));
ret.x += (u2.x - u1.x) * t, ret.y += (u2.y - u1.y) * t;
return ret;
}
void yuan(const Point & c, double r, const Point & l1, const Point & l2, Point & p1, Point & p2) {
Point p = c;
p.x += l1.y - l2.y;
p.y += l2.x - l1.x;
p = solve(p, c, l1, l2);
double t = sqrt(r * r - dis(p, c) * dis(p, c)) / dis(l1, l2);
p1.x = p.x + (l2.x - l1.x) * t;
p1.y = p.y + (l2.y - l1.y) * t;
p2.x = p.x - (l2.x - l1.x) * t;
p2.y = p.y - (l2.y - l1.y) * t;
}
void find(const Point & c1, double r1, const Point & c2, double r2, Point & p1, Point & p2) {
Point u, v;
double t = (1 + (r1 * r1 - r2 * r2) / dis(c1, c2) / dis(c1, c2)) / 2;
u.x = c1.x + (c2.x - c1.x) * t;
u.y = c1.y + (c2.y - c1.y) * t;
v.x = u.x + c1.y - c2.y;
v.y = u.y - c1.x + c2.x;
yuan(c1, r1, u, v, p1, p2);
}
double xmult(const Point & p1, const Point & p2, const Point & p0) {
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
double Dis(const Point & p, const Point & l1, const Point & l2) {
return fabs(xmult(p, l1, l2)) / dis(l1, l2);
}
int Get(int i, int j, Point p[2]) {
if(node[i].type == LINE && node[j].type == LINE) {
if(work(node[i].p[0], node[i].p[1], node[j].p[0], node[j].p[1])) {
return 0;
}
p[0] = solve(node[i].p[0], node[i].p[1], node[j].p[0], node[j].p[1]);
return 1;
}
if(node[i].type == CIRCLE && node[j].type == CIRCLE) {
if(dis(node[i].O, node[j].O) - EPS >= node[i].r + node[j].r) return 0;
if(dis(node[i].O, node[j].O) + EPS <= fabs(node[i].r - node[j].r)) return 0;
find(node[i].O, node[i].r, node[j].O, node[j].r, p[0], p[1]);
return 2;
} if(node[i].type == CIRCLE) swap(i, j);
if(Dis(node[j].O, node[i].p[0], node[i].p[1]) - EPS >= node[j].r) return 0;
yuan(node[j].O, node[j].r, node[i].p[0], node[i].p[1], p[0], p[1]);
return 2;
}
bool checkp(const Point &p) {
if(node[2].type == LINE)
return Dis(p, node[2].p[0], node[2].p[1]) <= EPS;
else return fabs(dis(p, node[2].O) - node[2].r) <= EPS;
}
signed main() { Point pp;
for (int i = 0; i < 3; ++ i) {
scanf("%d %d %d", x + i, y + i, r + i);
} pp.x = x[0], pp.y = y[0];
for (int i = 0; i < 3; ++ i) {
getget(i, (i + 1) % 3, i);
} Point p[2]; int np = Get(0, 1, p);
if(np == 1) printf("%.5lf %.5lf\n", p[0].x, p[0].y);
else if(np == 2) {
if(dis(p[0], pp) < dis(p[1], pp))
printf("%.5lf %.5lf\n", p[0].x, p[0].y);
else printf("%.5lf %.5lf\n", p[1].x, p[1].y);
}
return 0;
} 
京公网安备 11010502036488号