LA 5007 Detector Placement 模拟

题意：

给出一束光线(射线)，和一块三角形的棱镜 以及 棱镜的折射率，问光线能否射到X轴上，射到X轴上的坐标是多少。

分析：

其实直接模拟就好了，注意到题目中说不会发生全反射，所以如果射到棱镜中的话就一定能射出来。
一开始判断一下能否经过棱镜折射，不能的话直接算和X轴有没有交点或者交点的坐标。

1. 然后就是根据入射光线T1求入射点P1，注意直线可能和三角形的两条边都有交点，但最近的那个才是入射点。
   找到入射点就求法线，算角度，利用折射公式算折射角，求出折射光线T2。
2. 第二部分其实和上面的过程是一样的，求出射点，根据公式算出射角，最后求得出射光线。
3. 出射光线和X轴求交点，或者没有交点。

上个样例的光路图，仅供参考：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;

const double eps = 1e-8;

int dcmp(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}

struct Point
{
    double x, y;

    void read() { scanf("%lf%lf", &x, &y); }

    void print() { printf(" (%.2f, %.2f) 
", x, y); }

    Point(double x = 0, double y = 0): x(x), y(y) {}
};

typedef Point Vector;

Point operator + (const Point& A, const Point& B) {
    return Point(A.x + B.x, A.y + B.y);
}

Point operator - (const Point& A, const Point& B) {
    return Point(A.x - B.x, A.y - B.y);
}

Point operator * (const Point& A, double p) {
    return Point(A.x * p, A.y * p);
}

Point operator / (const Point& A, double p) {
    return Point(A.x / p, A.y / p);
}

double Dot(const Vector& A, const Vector& B) {
    return A.x * B.x + A.y * B.y;
}

double Cross(const Vector& A, const Vector& B) {
    return A.x * B.y - A.y * B.x;
}

double Length(Vector A) { return sqrt(Dot(A, A)); }

double Distance(Point A, Point B) { return Length(A - B); }

Vector Normal(Vector v) { return Vector(-v.y, v.x); }

Vector Rotate(Vector A, double rad) {
    return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));
}

struct Line
{
    Point P;
    Vector v;
    Line() {}
    Line(Point P, Vector v) : P(P), v(v) {}
};

bool OnSegment(Point P, Point a1, Point a2) {
    Vector v1 = a1 - P, v2 = a2 - P;
    return dcmp(Cross(v1, v2)) == 0 && dcmp(Dot(v1, v2)) < 0;
}

bool isParallel(const Line& A, const Line& B) {
    return dcmp(Cross(A.v, B.v)) == 0;
}

Point GetLineIntersection(Line L1, Line L2, double& t)
{
    Vector u = L1.P - L2.P;
    t = Cross(L2.v, u) / Cross(L1.v, L2.v);
    return L1.P + L1.v * t;
}

double u;
Point p[3], s, t;
Line L[3];

int main()
{
    Point hehe(0, 0), haha(10, 0);
    Line flat(hehe, haha);

    int T; scanf("%d", &T);
    while(T--) {
        s.read(); t.read();
        for(int i = 0; i < 3; i++) p[i].read();
        scanf("%lf", &u);
        L[0] = Line(p[0], p[1] - p[0]);
        L[1] = Line(p[1], p[2] - p[1]);
        L[2] = Line(p[2], p[0] - p[2]);

        Line T(s, t - s);    //入射光线
        int in_id = -1;
        Point in_P;         //入射点
        for(int i = 0; i < 3; i++) {
            double t;
            if(isParallel(T, L[i])) continue;    //直线平行，没有交点
            Point p = GetLineIntersection(T, L[i], t);
            if(dcmp(t) <= 0) continue;
            if(!OnSegment(p, L[i].P, L[i].P + L[i].v)) continue;    //交点不在线段上

            if(in_id == -1 || Distance(s, p) < Distance(s, in_P)) {
                in_P = p; in_id = i;
            }
        }

        if(in_id == -1) {    //没有射到棱镜上
            double t;
            if(dcmp(T.v.y) >= 0) printf("Error
");
            else printf("%.3f
", GetLineIntersection(T, flat, t).x);
            continue;
        }

        //printf("in_Point:");
        //in_P.print();

        Vector normal = Normal(L[in_id].v);    //计算第一次折射的法线
        if(dcmp(Dot(normal, T.v)) < 0) normal.x = -normal.x, normal.y = -normal.y;    //调整法线的方向
        double theta = asin(fabs(Cross(T.v, normal)) / Length(T.v) / Length(normal) / u);//计算折射角
        Line T2; T2.P = in_P;    //折射光线
        if(dcmp(Cross(normal, T.v)) > 0) T2.v = Rotate(normal, theta);
        else T2.v = Rotate(normal, -theta);

        int out_id;
        Point out_P;    //出射点
        for(int i = 0; i < 3; i++) if(i != in_id) {
            if(isParallel(T2, L[i])) continue;
            double t;
            Point p = GetLineIntersection(T2, L[i], t);
            if(dcmp(t) <= 0) continue;
            if(!OnSegment(p, L[i].P, L[i].P + L[i].v)) continue;
            out_P = p;
            out_id = i;
        }

        //printf("out_Point:");
        //out_P.print();

        Vector normal2 = Normal(L[out_id].v);
        if(dcmp(Dot(T2.v, normal2)) < 0) normal2.x = -normal2.x, normal2.y = -normal2.y;
        double theta2 = asin(fabs(Cross(T2.v, normal2)) / Length(T2.v) / Length(normal2) * u);
        Line T3; T3.P = out_P;
        if(dcmp(Cross(normal2, T2.v)) > 0) T3.v = Rotate(normal2, theta2);
        else T3.v = Rotate(normal2, -theta2);

        double t;
        if(dcmp(T3.v.y) >= 0) printf("Error
");
        else printf("%.3f
", GetLineIntersection(T3, flat, t).x);
    }

    return 0;
}