zjoi 2007 particle 粒子运动 计算几何

题意：阿Q博士正在观察一个圆形器皿中的粒子运动。不妨建立一个平面直角坐标系，圆形器皿的圆心坐标为(x0, y0)，半径为R。器皿中有若干个粒子，假设第i个粒子在时刻0的位置为(xi, yi)，速度为(vxi,vyi)（注：这是一个速度向量，若没有发生碰撞，t时刻的位置应该是(xi + t * vxi, yi + t * vyi) ）。假设所有粒子的运动互不干扰；若某个粒子在某个时刻碰到了器皿壁，将发生完全弹性碰撞，即速度方向按照碰撞点的切线镜面反射，且速度大小不变（如图）。认为碰撞是瞬间完成的。  尽管碰撞不会影响粒子的速率，但是粒子却会受到一定的伤害，所以若某一个粒子碰撞了k次器皿壁，那么在第k次碰撞时它便会消亡。 出于研究的需要，阿Q博士希望知道从时刻0到所有粒子都消亡这段时间内，所有粒子之间的最近距离是什么。你能帮助他么？

对于所有的数据，2 ≤N ≤100。1≤k ≤100。

思路：计算几何

数据规模很小 每2个点分别计算

计算碰撞时间 对于每一段时间计算最小值

时间复杂度 n*n*k

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define MAXN 110
#define EPS 1e-8
struct point
{
    double x,y;
    point(){}
    point(double x0,double y0): x(x0),y(y0) {}
};
struct crash
{
    double t;
    point start,v;
}a[MAXN][MAXN];
point start[MAXN],v[MAXN],circle;
int n,k;
double r,ans=987654321;
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 double &A,const point &B)
{
    return point(A*B.x,A*B.y);
}
double cross(point A,point B)
{
    return A.x*B.y-B.x*A.y;
}
double length(point A)
{
    return sqrt(A.x*A.x+A.y*A.y);
}
double find_intersect(int i,int j)
{
    point temp=a[i][j].start-circle;
    double A,B,C;
    A=a[i][j].v.x*a[i][j].v.x+a[i][j].v.y*a[i][j].v.y;
    B=2*a[i][j].v.x*temp.x+2*a[i][j].v.y*temp.y;
    C=temp.x*temp.x+temp.y*temp.y-r*r;
    return (-B+sqrt(B*B-4*A*C))/2/A;
}
void make_crash(int i)
{
    int j;
    a[i][0].t=0;
    double t,sinn,cosn,sin2n,cos2n;
    point temp,v;
    for(j=0;j<k;j++)
    {
        t=find_intersect(i,j);
        a[i][j+1].t=a[i][j].t+t;
        temp.x=a[i][j].start.x+t*a[i][j].v.x;
        temp.y=a[i][j].start.y+t*a[i][j].v.y;
        a[i][j+1].start=temp;
        sinn=cross(a[i][j].start-temp,circle-temp)/r/length(a[i][j].start-temp);
        cosn=sqrt(1-sinn*sinn);
        sin2n=2*sinn*cosn;
        cos2n=cosn*cosn-sinn*sinn;
        v.x=a[i][j].v.x*cos2n-sin2n*a[i][j].v.y;
        v.y=a[i][j].v.x*sin2n+cos2n*a[i][j].v.y;
        v.x=-v.x; v.y=-v.y;
        a[i][j+1].v=v;
    }
    a[i][k+1].t=987654321;
}
double calc_time(point start_i,point start_j,point v_i,point v_j, double t)
{
    point v=v_i-v_j;
    point temp=start_i-start_j;
    double A,B,C,ans=987654321;
    A=v.x*v.x+v.y*v.y;
    B=2*v.x*temp.x+2*v.y*temp.y;
    C=temp.x*temp.x+temp.y*temp.y;
    if(A*t*t+B*t+C<ans) ans=A*t*t+B*t+C;
    if(C<ans) ans=C;
    if(A*B<=0&&A!=0&&B/(-2*A)<=t)
        ans=(4*A*C-B*B)/4/A;
    return sqrt(ans);
}    
double solve(int i,int j)
{
    int p,q;
    double ans=987654321,start_time,temp;
    point start_i,start_j;
    p=q=0;
    while(p<k||q<k)
    {
        start_time=max(a[i][p].t,a[j][q].t);
        start_i=a[i][p].start+(start_time-a[i][p].t)*a[i][p].v;
        start_j=a[j][q].start+(start_time-a[j][q].t)*a[j][q].v;
        if(a[i][p+1].t<a[j][q+1].t) 
        {
            temp=calc_time(start_i,start_j,a[i][p].v,a[j][q].v,a[i][p+1].t-start_time);
            if(temp<ans) 
                ans=temp;
            p++;
            continue;
        }
        else 
        {
            temp=calc_time(start_i,start_j,a[i][p].v,a[j][q].v,a[j][q+1].t-start_time);
            if(temp<ans) ans=temp;
            q++;
            continue;
        }
    }
    return ans;
}

int main()
{
    double temp;
    freopen("particle.in","r",stdin);
    freopen("particle.out","w",stdout);
    scanf("%lf%lf%lf",&circle.x,&circle.y,&r);
    scanf("%d%d",&n,&k);
    int i,j;
    for(i=1;i<=n;i++)
        scanf("%lf%lf%lf%lf",&a[i][0].start.x,&a[i][0].start.y,&a[i][0].v.x,&a[i][0].v.y);
    for(i=1;i<=n;i++) make_crash(i);
    for(i=1;i<=n;i++)
    for(j=i+1;j<=n;j++)
        ans=min(ans,solve(i,j));
    printf("%.3lf\n",ans);
    return 0;
}