推荐理由:
细节较多的计算几何,适合一做。
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; }