Category : Geometry
Solution:
First we need to find the intersection between a sphere and a line.
a line is denoted by an origin and a direction vector.
a sphere is denoted a centre and radius. To find the intersection ,lets consider the equation of sphere.
Let p be a point on the sphere. Then the following equation holds
where c is the centre of the sphere and r is radius.
Now if the line intersects the sphere in point p, then we can write p as
where o is origin of the line and d is direction unit vector of the line. x is the scaling factor.
putting p in equation (1)
This is a quadratic equation of x in form
where
if discriminant >= 0 , then the line intersects the sphere.
We can find x.
We will only take the value of x which is greater than ZERO (why ? )
After finding x , We can find the intersecting point.
Now we need to rotate the direction vector. This is shown in the following ways Reflection
/*
*************************
Id : Matrix.code
Task:
Date: 2016-02-04
**************************
*/
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
/*------- Constants---- */
#define Long long long
#define ull unsigned long long
#define mod 1000000007
#define MEMSET_INF 63
#define MEM_VAL 1061109567
#define forn(i,n) for( int i=0 ; i < n ; i++ )
#define mp(i,j) make_pair(i,j)
#define lop(i,a,b) for( int i = (a) ; i < (b) ; i++)
#define pb(a) push_back((a))
#define SZ(a) (int) a.size()
#define all(x) (x).begin(),(x).end()
#define gc getchar_unlocked
#define PI acos(-1.0)
#define INF 1<<29
#define EPS 1e-12
#define Fr first
#define Sc second
#define Sz size()
#define lc ((n)<<1)
#define rc ((n)<<1|1)
#define db(x) cout << #x << " -> " << x << endl;
#define Di(n) int n;si(n)
#define Dii(a,b) int a,b;si(a);si(b)
#define Diii(a,b,c) int a,b,c;si(a);si(b);si(c)
#define Si(n) si(n)
#define Sii(a,b) si(a);si(b)
#define Siii(a,b,c) si(a);si(b);si(c)
#define min(a,b) ((a)>(b) ? (b) : (a) )
#define max(a,b) ((a)>(b) ? (a):(b))
/*---- short Cuts ------- */
#define ms(ara_name,value) memset(ara_name,value,sizeof(ara_name))
typedef pair<int, int> ii;
typedef vector<int> vi ;
typedef vector<Long> vl;
/*------ template functions ------ */
#ifndef getchar_unlocked
#define getchar_unlocked getchar
#endif
template<class T> inline void si(T &x){register int c = gc();x = 0;int neg = 0;for(;((c<48 | c>57) && c != '-');c = gc());
if(c=='-') {neg=1;c=gc();}for(;c>47 && c<58;c = gc()) {x = (x<<1) + (x<<3) + c - 48;}if(neg) x=-x;}
Long bigmod(Long p,Long e,Long M){Long ret = 1;for(; e > 0; e >>= 1){if(e & 1) ret = (ret * p) % M;p = (p * p) % M;} return ret;}
Long gcd(Long a,Long b){if(b==0)return a;return gcd(b,a%b);}
Long modinverse(Long a,Long M){return bigmod(a,M-2,M);}
void io(){freopen("/Users/MyMac/Desktop/in.txt","r",stdin);}
/*************************** END OF TEMPLATE ****************************/
#define Z(x) (abs(x) < EPS)
#define ZN(x) ( x < 0 || Z(x))
#define ZP(x) ( x > 0 || Z(x))
const int N = 1005;
typedef double T;
struct PT{
T x,y,z;
PT(){}
PT(T x,T y,T z):x(x),y(y),z(z){}
void scan(){cin >> x >> y >> z;}
PT operator + (const PT &p) const{ return PT(x+p.x,y+p.y,z+p.z);}
PT operator - (const PT &p) const{ return PT(x-p.x,y-p.y,z-p.z);}
PT operator * (const double c) const {return PT(x*c,y*c,z*c);}
double abs(){ return sqrt(x*x+y*y+z*z);}
void normal(){
double d = this->abs();
if(fabs(d) < EPS) return;
x/=d;
y/=d;
z/=d;
}
void print() { printf("(%f %f %f)",x,y,z);}
};
double dot(PT a,PT b){ return a.x*b.x+a.y*b.y+a.z*b.z;}
struct Sp{
PT cen;
double rad;
Sp(){}
};
vector<Sp> vs;
int main()
{
//freopen("in.txt","r",stdin);
int n;
while(scanf("%d",&n)==1) {
if(n==0) break;
PT org(0,0,0);
PT dir;
dir.scan();
dir.normal();
Sp sp;
for(int i= 0;i < n ;i ++) {
sp.cen.scan();
cin >> sp.rad;
vs.pb(sp);
}
PT ans;
while(1) {
vector<pair<double,int> > Map;
for(int i= 0;i < n ;i ++) {
double a = dot(dir,dir);
double b = 2 * dot(dir, org - vs[i].cen);
double c = dot(org-vs[i].cen , org - vs[i].cen) - vs[i].rad * vs[i].rad;
double disc = b * b - 4 * a * c;
if(ZP(disc) ) {
double x1 = -b - sqrt(disc);
double x2 = -b + sqrt(disc);
x1/=(2*a); x2/=(2*a);
vector<double> k;
if(x1>EPS) k.pb(x1);
if(x2 > EPS) k.pb(x2);
sort(all(k));
if(k.size()) {
Map.pb(mp(k[0],i));
}
}
}
if(Map.size() == 0 ) break;
sort(all(Map));
double x = Map[0].first;
double col = Map[0].second;
org = org + dir * x;
PT nor = org - vs[col].cen;nor.normal();
dir = dir - nor * 2 * dot(dir,nor);
ans = org;
Map.clear();
}
printf("%f %f %f\n",ans.x, ans.y,ans.z);
vs.clear();
}
return 0;
}
Code for Fun 