为防止广告，目前nocow只有登录用户能够创建新页面。如要创建页面请先登录/注册（新用户需要等待1个小时才能正常使用该功能）。
Sgu/110
来自NOCOW
< Sgu
简单的三维计算几何.需要的功能有:直线和圆的交点, 镜面反射.
//P*M #include<cstdio> #include<cstring> #include<cmath> const double zero = 1e-8;//这题其实不用zero struct triple{ double x, y, z; triple(){ x = 0; y = 0; z = 0; } triple(double xx, double xy, double xz){ x = xx; y = xy; z = xz; } }; double sqr(double a){ return a*a; } typedef triple point; struct line{ triple s, d; line(){ } line(triple xt1, triple xt2){ s = xt1; d= xt2; } }; struct circle{ triple h; double r; }; circle cir[100]; line lin1; point co1, co2; bool flag; int n; triple operator +(triple t1, triple t2){ return triple(t1.x+t2.x, t1.y+t2.y, t1.z+t2.z); } triple operator -(triple t1, triple t2){ return triple(t1.x-t2.x, t1.y-t2.y, t1.z-t2.z); } triple operator -=(triple &t1, triple t2){ t1 = t1 - t2; return t1; } triple operator *(double r, triple t){ return triple(t.x*r, t.y*r, t.z*r); } double operator *(triple t1, triple t2){ return t1.x*t2.x+t1.y*t2.y+t1.z*t2.z; } double dis(point p1, point p2){ return sqrt(sqr(p1.x-p2.x)+sqr(p1.y - p2.y)+sqr(p1.z-p2.z)); } double dis(point p, line lin1){ double t = (lin1.d.x*(p.x-lin1.s.x)+lin1.d.y*(p.y-lin1.s.y)+lin1.d.z*(p.z-lin1.s.z))/(sqr(lin1.d.x)+sqr(lin1.d.y)+sqr(lin1.d.z)); return dis(p, t*lin1.d+lin1.s); } void cross_over(circle cir1, line lin1, point &co1, point &co2){//直接联立方程求解射线和圆的交点 double a=sqr(lin1.d.x)+sqr(lin1.d.y)+sqr(lin1.d.z), b=2*(lin1.d.x*(lin1.s.x-cir1.h.x)+lin1.d.y*(lin1.s.y-cir1.h.y)+lin1.d.z*(lin1.s.z-cir1.h.z)), c = sqr(lin1.s.x-cir1.h.x)+sqr(lin1.s.y-cir1.h.y)+sqr(lin1.s.z-cir1.h.z)-sqr(cir1.r); double root1 = (-b+sqrt(b*b-4*a*c))/2/a; double root2 = (-b-sqrt(b*b-4*a*c))/2/a; co1 = triple(lin1.s.x+root1*lin1.d.x, lin1.s.y+root1*lin1.d.y, lin1.s.z+root1*lin1.d.z); co2 = triple(lin1.s.x+root2*lin1.d.x, lin1.s.y+root2*lin1.d.y, lin1.s.z+root2*lin1.d.z); } triple mirror(triple ray, triple normal){ double t = ((normal.x*ray.x)+(normal.y*ray.y)+(normal.z*ray.z))/(normal.x*normal.x+normal.y*normal.y+normal.z*normal.z); triple along = t*normal; triple vertical = ray - along; return vertical - along;//镜面,不是对称 } int main(){ freopen("sgu110.in", "r", stdin); freopen("sgu110.out", "w", stdout); scanf("%d\n", &n); for(int i = 1; i <= n; i++){ scanf("%lf %lf %lf %lf", &cir[i].h.x, &cir[i].h.y, &cir[i].h.z, &cir[i].r); } scanf("%lf %lf %lf %lf %lf %lf", &lin1.s.x, &lin1.s.y, &lin1.s.z, &lin1.d.x, &lin1.d.y, &lin1.d.z); lin1.d -= lin1.s; int bak = -1; for(int i = 1; i <= 10; i++){ flag = false; triple comin(100000, 100000, 100000); int cirmin; for(int j = 1; j <= n; j++){ if(dis(cir[j].h, lin1) - cir[j].r<zero){ cross_over(cir[j], lin1, co1, co2); if((co1-lin1.s) * (lin1.d) > -zero && bak != j){ flag = true; if(dis(lin1.s, co1) - dis(lin1.s, comin)< zero){comin = co1; cirmin = j;} if(dis(lin1.s, co2) - dis(lin1.s, comin) < zero){comin = co2; cirmin = j;} } } } if(flag){ if(i==1){printf("%d", cirmin);}else{printf(" %d", cirmin);} bak = cirmin; lin1 = line(comin, mirror(lin1.d, comin-cir[cirmin].h)); }else{ printf("\n"); break; } } flag = false; triple comin(100000, 100000, 100000); int cirmin; for(int j = 1; j <= n; j++){ if(dis(cir[j].h, lin1) - cir[j].r<zero){ cross_over(cir[j], lin1, co1, co2); if((co1-lin1.s) * (lin1.d) > -zero && bak != j){ flag = true; if(dis(lin1.s, co1) - dis(lin1.s, comin)<zero){comin = co1; cirmin = j;} if(dis(lin1.s, co2) - dis(lin1.s, comin)<zero){comin = co2; cirmin = j;} } } } if(flag){ printf(" etc.\n"); }else{ } fclose(stdin); fclose(stdout); return 0; }
//by hza #include<cstdio> #include<cmath> const int MAX=100; const double INF=1e10; const double EPS=1e-10; int answer[MAX]; struct ball { double x,y,z,r; }b[MAX]; int n; double sqr(double x){return x*x;} struct vector { double x,y,z; void print(){printf("%lf %lf %lf\n",x,y,z);} double dist(){return sqrt(sqr(x)+sqr(y)+sqr(z));} void change() { double len=dist(); x/=len;y/=len;z/=len; } void change(double k) { double len=dist(); x/=len;y/=len;z/=len; x*=k;y*=k;z*=k; } void operator =(const ball& b){x=b.x;y=b.y;z=b.z;} void operator =(const vector& b){x=b.x;y=b.y;z=b.z;} }; vector operator + (const vector& a,const vector& b) { vector next; next.x=a.x+b.x; next.y=a.y+b.y; next.z=a.z+b.z; return next; } struct node { double x,y,z; void operator =(const node& b){x=b.x;y=b.y;z=b.z;} void operator =(const ball& b){x=b.x;y=b.y;z=b.z;} void operator =(const vector& b){x=b.x;y=b.y;z=b.z;} void print(){printf("%lf %lf %lf\n",x,y,z);} }first,second; void get_vector(vector& now,node &first,node &last) { now.x=last.x-first.x; now.y=last.y-first.y; now.z=last.z-first.z; } double vector_ball(node now,vector v,ball ba) { double a,b,c; a=sqr(v.x)+sqr(v.y)+sqr(v.z); b=2*(now.x-ba.x)*v.x+2*(now.y-ba.y)*v.y+2*(now.z-ba.z)*v.z; c=sqr(now.x-ba.x) +sqr(now.y-ba.y) +sqr(now.z-ba.z)-sqr(ba.r); if(b*b<4*a*c)return INF; double x1=(-b-sqrt(b*b-4*a*c))/a*0.5; double x2=(-b+sqrt(b*b-4*a*c))/a*0.5; if(x1>0)return x1; if(x2>0)return x2; return INF; } node get_node(node& begin,vector& v,double len) { v.change(); node a; a.x=begin.x+len*v.x; a.y=begin.y+len*v.y; a.z=begin.z+len*v.z; return a; } double dianji(vector& a,vector& b){return a.x*b.x+a.y*b.y+a.z*b.z;} int main() { #ifndef ONLINE_JUDGE freopen("110.in","r",stdin);freopen("110.out","w",stdout); #endif int i; scanf("%d",&n); for(i=1;i<=n;++i)scanf("%lf %lf %lf %lf",&b[i].x,&b[i].y,&b[i].z,&b[i].r); scanf("%lf %lf %lf %lf %lf %lf",&first.x,&first.y,&first.z,&second.x,&second.y,&second.z); vector now,nv,in;//nv法向量 in入射光线 node begin=first,sec,center; get_vector(now,first,second); int k=0; double t,nearest=INF; double d; for(k=1;k<=11;++k) { nearest=INF; now.change(); for(i=1;i<=n;++i) if(i!=answer[k-1]) { t=vector_ball(begin,now,b[i]); if(t<nearest){nearest=t;answer[k]=i;} } if(fabs(nearest-INF)<EPS)break; t=nearest; sec=get_node(begin,now,t); center=b[answer[k]]; get_vector(nv,center,sec); get_vector(in,begin,sec); in.change(); d=fabs(dianji(in,nv))/nv.dist();d*=2; nv.change(d); now=nv+in; begin=sec; } for(i=1;i<=11;++i) { if(answer[i]==0) break; if(i<=10) printf("%d ",answer[i]); else printf("etc."); } printf("\n"); }
