ray tracing.cpp

基于snell定律的射线追踪

//ray tracing
#include"math.h"
#include"iostream"
#include"fstream"
#include "vector"
using namespace std;
class CPointd
{
public:
	double x;
	double y;
	double V;
	CPointd(){x=0,y=0;}
};
int main()
{
     int i,j,xi,yj,NS,NG,n=0,N=500;
	 char opfname[20];ifstream ifile;
	 double f=1.0,df=2,e=100000;int a=0,b=0;
	 double temp,tangle0,tangle1,Pa,v;
	 CPointd *p,p0,p1,p2;bool judge=true;
	 vector <CPointd>  pp;
	 cin>>opfname;
     ifile.open(opfname,ios::in|ios::binary);
	 ifile.read((char*)&NS,4);
     ifile.read((char*)&NG,4);
     p=new CPointd[NS*NG];
	 for(i=0;i<NS;i++)
		 for(j=0;j<NG;j++)
		 {
               ifile.read((char*)&((p+i*NG+j)->x),8);
               ifile.read((char*)&temp,8);
               ifile.read((char*)&temp,8);
               ifile.read((char*)&((p+i*NG+j)->y),8);
               ifile.read((char*)&((p+i*NG+j)->V),8);
			   (p+i*NG+j)->y=-(p+i*NG+j)->y;  
			   (p+i*NG+j)->V=3000;//velocity
		 }	 
     cin>>xi>>yj;
     int kx=0,ky=yj;
	 while(fabs(e)>0.000000001&&n<=N)
	 {   
	     pp.clear();
		 judge=true;
		 kx=0,ky=yj;
		 p0.x=0;
		 p0.y=(p+ky)->y;
		 p2.x=f*(p+xi*NG)->x;
		 p2.y=0;
		 tangle1=(p2.x-p0.x)/(p2.y-p0.y);
		 Pa=(p+kx*NG+ky)->V*sqrt(pow(p2.x,2)+pow(p0.y,2))/p2.x;
		 p1.x=(p+kx*NG+ky-1)->x,p1.y=(p+kx*NG+ky-1)->y; 
		 tangle0=(p1.x-p0.x)/(p1.y-p0.y);
		 pp.push_back(p0);
		 while(judge==true)
		 { 
			 if(tangle1<=tangle0)
			 {
				 if(ky>=1)
				 { 
					 ky--;
					 p2.y=(p+ky)->y;
				 }
				 else 
				 {
					 p2.y=0;	judge=false;
				 }
				 p2.x=(p2.y-p0.y)*tangle1+p0.x;
				 if(fabs(tangle1-tangle0)<0.0000000000001)
				 {	 kx++;
				     if(kx>xi)
					    judge=false;
				 }
			 }
			 else
			 {
				 p2.x=(p+kx*NG)->x;
				 p2.y=(p2.x-p0.x)/tangle1+p0.y; 
				 kx++;
				 if(kx>xi)
					 judge=false;
			 }
			 if(ky>=0)
			 {p1.x=(p+kx*NG)->x,p1.y=(p+ky-1)->y;}
			 else 
			 {
				 p1.x=(p+kx*NG)->x,p1.y=0;
			 }
			 p0.x=p2.x,p0.y=p2.y;
			 v=(p+kx*NG+ky)->V;			
			 tangle0=(p1.x-p0.x)/(p1.y-p0.y);
			 if(sqrt(pow(Pa/v,2)-1)!=0)
				 tangle1=1.0/sqrt(pow(Pa/v,2)-1);
			 else 
			 {
			     judge=false;
			 }
             pp.push_back(p0);
		 }
		 if(p0.x<(p+xi*NG)->x&&sqrt(pow(Pa/v,2)-1)>0)
		 {		 
			if(a!=0)
			{
				df=df/2;
			}
			f=f+df;
			b=1;
			
		 }
		 else if(p0.y<0)
		 {
			if(b==0)
				f=f/2.0;
			else 
			{df=df/2;
			f=f-df;}
			 a=1;	 
		 }
		 n++;
		 e=sqrt(pow(p0.x-(p+xi*NG)->x,2)+pow(p0.y,2));
	 }
     if(fabs(e)<0.000000001)
	 {
		 cout<<pp.size()<<endl;
	     for(i=0;i<pp.size();i++)
		 cout<<pp[i].x<<"      "<<-pp[i].y<<endl;
	 }
	 else if(fabs(e)>0.000000001)
		 cout<<"这种速度分布两点不可进行透射追踪！"<<endl;
	 p=NULL;
	return 0;
}