标签:c style class blog code java
这道题是计算几何,这是写的第一道计算几何,主要是难在如何求入射光线的反射光线。
我们可以用入射光线 - 入射光线在法线(交点到圆心的向量)上的投影*2 来计算反射光线,自己画一个图,非常清晰明了。
具体到程序里,我们可以 v2 = v1 - fa / Length(fa) * 2 * ( Dot(v1, fa) / Length(fa)) 来求,简单来说就是用内积(点积)求出入射光线在法线上的长度,然后用法线的单位向量乘这个长度就可以了。
#include <cstdio> #include <cstring> #include <cstdlib> #include <iostream> #include <algorithm> #include <cmath> #define N 55 #define inf 0x7f7f7f7f using namespace std; struct Point3 { double x, y, z; Point3(double x=0, double y=0, double z=0):x(x), y(y), z(z) { } }; typedef Point3 Vector3; const double eps = 1e-9; int dcmp(double x) { if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } int n; double r[N]; Point3 Begin, circle[N]; Vector3 Direct; Vector3 operator + (Vector3 A, Vector3 B) { return Vector3(A.x+B.x, A.y+B.y, A.z+B.z); } Vector3 operator - (Point3 A, Point3 B) { return Vector3(A.x-B.x, A.y-B.y, A.z-B.z); } Vector3 operator * (Vector3 A, double p) { return Vector3(A.x*p, A.y*p, A.z*p); } Vector3 operator / (Vector3 A, double p) { return Vector3(A.x/p, A.y/p, A.z/p); } double min(double x, double y) { return x < y ? x : y; } double Dot (Vector3 A, Vector3 B) { return A.x*B.x + A.y*B.y + A.z*B.z; } double Length (Vector3 A) { return sqrt(Dot(A, A)); } double Angle (Vector3 A, Vector3 B) { return acos(Dot(A, B) / Length(A) / Length(B)); } Vector3 Cross (Vector3 A, Vector3 B) { return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.y*B.z - A.z*B.y); } double disPointLine(Point3 P, Point3 A, Point3 B) { Vector3 v1, v2; v1 = B - A; v2 = P - A; return Length(Cross(v1, v2)) / Length(v1); } int main() { scanf("%d", &n); double x, y, z; for (int i = 1; i <= n; ++i) { scanf("%lf%lf%lf%lf", &x, &y, &z, &r[i]); circle[i] = Point3(x, y, z); } scanf("%lf%lf%lf", &x, &y, &z); Begin = Point3(x, y ,z); scanf("%lf%lf%lf", &x, &y, &z); Direct = Point3(x, y, z) - Begin; int nowcircle = 0; for (int w = 1; w <= 11; ++w) { int nextcircle = 0; double mindis = inf; for (int i = 1; i <= n; ++i) { if (i == nowcircle) continue; double a, b, c; Point3 nc = circle[i]; a = Direct.x*Direct.x + Direct.y*Direct.y + Direct.z*Direct.z; b = 2 * ((Begin.x-nc.x) * Direct.x + (Begin.y-nc.y)*Direct.y + (Begin.z-nc.z)*Direct.z); c = (Begin.x-nc.x)*(Begin.x-nc.x) + (Begin.y-nc.y)*(Begin.y-nc.y) + (Begin.z-nc.z)*(Begin.z-nc.z) - r[i]*r[i]; double delta; delta = b*b - 4*a*c; if (delta < 0) continue; double ans1, ans2; ans1 = (-b + sqrt(delta)) / (2*a); ans2 = (-b - sqrt(delta)) / (2*a); if (dcmp(ans1) < 0 && dcmp(ans2) < 0) continue; else if (dcmp(ans1) < 0 && ans2 < mindis) { mindis = ans2; nextcircle = i; } else if (dcmp(ans2) < 0 && ans1 < mindis) { mindis = ans1; nextcircle = i; } else if (min(ans1, ans2) < mindis) { mindis = min(ans1, ans2); nextcircle = i; } } if (!nextcircle) break; else if (w == 11) { printf(" etc."); break; } else { if (w == 1) printf("%d", nextcircle); else printf(" %d", nextcircle); Point3 jiao; Vector3 v1, v2, fa; jiao = Begin + Direct*mindis; v1 = Direct; fa = circle[nextcircle] - jiao; v2 = (fa / Length(fa)) * (2 * Dot(v1, fa) / Length(fa)); nowcircle = nextcircle; Begin = jiao; Direct = v1 - v2; } } printf("\n"); }
sgu 110 Dungeon,布布扣,bubuko.com
标签:c style class blog code java
原文地址:http://www.cnblogs.com/handsomeJian/p/3765456.html