/* 兩點射線追蹤，設計者：胡治權  原理：參考楊文采《地球物理反演理論》，快速精準射線追蹤*/
/*功能：完成兩點射線追蹤，并將射線在剖面上的交點寫入磁盤文件*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PAI 3.1415926

struct dotinform
{ 
	float x;
	float z;
	float nx;
	float nz;
	float tsum;
};/*存放點信息*/

struct jie
{ 
	float s;
	char ch;
	float val;
};/*存放解信息*/
int deltv(float x,float z,float nx,float nz,float *DVX,float *DVZ,float *vp,float *v,float fx,float dx,int xn,
		  float fz ,float dz,int zn);
float raytrace (float *X,float *Z,float *NX,float *NZ,float *T, float v,float dvx,float dvz,float dx,float dz,float fx,float fz,
			   struct dotinform *p,int *k,float xend, float zend);
void main()
{
	FILE *infp,*outfpxz,*outfpp;
	dotinform p[5000];
	char ch[20];
	int ii,ij,is;
	int vtype;
	float *v;
	float fangle ,endangle,dangle;
	float fsx,dsx;
	float fx,dx;
	float fz,dz;
	int   xn,zn,nsx;
    int   k=1,vflag;
    float xend,zend,n,sign;
	float x,z,t,va=0,dvxa=0,dvza=0,nx,nz,alpha,sx;
   
	printf("*******************射線追蹤(設計：胡治權)***********************\n");
    printf("***選擇輸入的速度文件類型：1. 文本文件  2. 二進制文件\n");
    do{
	   scanf("%d",&vtype);
	   if(vtype != 1 && vtype != 2)
		  printf("輸入錯誤，請重新輸入\n");
	   }while(vtype !=1 && vtype != 2);
	do{
        printf("***輸入速度文件名：");
        scanf("%s",ch);
        if((infp = fopen(ch,"r")) == NULL)
			printf("無法打開指定的速度文件，請核對后重新輸入！\n");
	}while(infp == NULL);
	printf("***輸入速度網格剖分參數***\n");
	printf("橫向起始坐標(數據類型:float):fx=");
	scanf("%f",&fx);
	do{
	   printf("橫向采樣間隔(數據類型:float>0):dx=");
	   scanf("%f",&dx);
	   if(dx <= 0) printf("dx<=0輸入錯誤，請重新輸入！");
	   }while(dx<=0);
	do{
	   printf("橫向采樣點數(數據類型:int>0):xn=");
	   scanf("%d",&xn);
	   if(xn<=0)printf("xn<=0輸入錯誤，請重新輸入！");
	   }while(xn<=0);
    printf("縱向起始坐標(數據類型:float):默認設置fz=0\n");
    fz = 0;
	do{
	    printf("縱向采樣間隔(數據類型:float>0):dz=");
	    scanf("%f",&dz);
		if(dz<=0)printf("dz<=0輸入錯誤，請重新輸入！");
	}while(dz<=0);
	do{
	    printf("縱向采樣點數(數據類型:int>0):zn=");
		scanf("%d",&zn);
		if(zn<=0)printf("zn<=0輸入錯誤，請重新輸入！");
	}while(zn<=0);
    
	xend = fx +(xn-1)*dx;
	zend = fz +(zn-1)*dz;
	printf("***輸入炮參數***\n");
	printf("請輸入第一炮的橫坐標(數據類型float):fsx=");
	do{
	   scanf("%f",&fsx);
       if(fsx<fx || fsx >xend) printf("炮點不在剖面上，請重新輸入！");
	   }while(fsx <fx || fsx >xend);
   do{
	    printf("炮點采樣間隔(數據類型:float>0):dsx=");
	    scanf("%f",&dsx);
		if(dsx<=0)printf("dsx<=0輸入錯誤，請重新輸入！");
	}while(dsx<=0);
	do{
	    printf("炮點采樣點數(數據類型:int>0):nsx=");
		scanf("%d",&nsx);
		if(nsx<=0)printf("nsx<=0輸入錯誤，請重新輸入！");
	}while(nsx<=0);
	printf("***角度參數***\n");
	do{
		printf("輸入角度起始值(°):fangle=");
		scanf("%f",&fangle);
		printf("輸入角度中止值(°)>=起始值:endangle=");
		scanf("%f",&endangle);
		if(fangle>endangle)printf("輸入錯誤,請重新輸入！");
	}while(fangle>endangle);
	do{
		 printf("輸入角度增量(>=0):dangle=");
		 scanf("%f",&dangle);
		 if(dangle<0)printf("輸入錯誤，請重新輸入!");
	}while(dangle<0);
	printf("*******參數輸入完畢*******\n");
	printf("進行射線追蹤，請稍等...\n");


	outfpxz=fopen("xz.dat","w");
	outfpp=fopen("pshu.dat","w");

    
    v= (float*) malloc(xn*zn*sizeof(float));
 
   if(vtype == 1)
   {
    for (ii = 0 ; ii < zn; ii++)
	{

		for (ij = 0;ij < xn;ij++)
		{
			fscanf(infp,"%f",&v[ii*xn+ij]);
		
		}
       
	}
   }else
   {
	   fread(v,sizeof(float),zn*xn,infp);
   }
  
   for(is = 0;is< nsx ; is++)
   {
      sx = fsx + is*dsx;n=0;
	  
	  while ((fangle + n*dangle) <= endangle)
	  {
		  x = sx;z = 0;
		  alpha = (float)PAI*(fangle + n*dangle)/180;
		  printf("射線追蹤：炮點ｓｘ＝%f，入射角度ａｌｐｈａ＝%f，",sx,(fangle+n*dangle));
		  nx = -(float)sin(alpha);
	      nz = (float)cos(alpha)  ; 
	      t  = 0;                      
	      p[0].x = x;    
	      p[0].z = 0;
         
	    do
		{
          vflag = deltv(x,z,nx,nz,&dvxa,&dvza,&va,v,fx,dx,xn,fz,dz,zn);
		  if(vflag == 0) break;
          sign=raytrace(&x,&z,&nx,&nz,&t,va,dvxa,dvza,dx,dz,fx,fz,p,&k,xend,zend);
		  if(sign == 0 || sign == -1) break;
		}while(1);
	
        
		if(sign != 0 )
        {
          for (ii = 0;ii < k ;ii++)
		  {
		   fprintf(outfpxz,"%.8f  ",p[ii].x);
	       fprintf(outfpxz,"%.8f\n ",p[ii].z);
		  }
		  fprintf(outfpp,"%d\n",k);
		  printf("射線交點數ｋ＝%d\n",k);
		}else{ printf("追蹤失敗!\n");}
		 k = 1;
		 n = n + 1;
	  }
    printf("完成炮點ｓｘ＝%f的射線追蹤\n",sx);
   }


    printf("完成射線追蹤！");
    fclose(infp);
	fclose(outfpxz);
	fclose(outfpp);

}
int deltv(float x,float z,float nx,float nz,float *DVX,float *DVZ,float *vp,float *v,float fx,float dx,int xn,
		  float fz ,float dz,int zn)
{
	float l,dvx,dvz;
	float v1,v2,v3,v4;
	float dv1,dv2,dv1x,dv1z,dv2x,dv2z;
	int i,j;

	
	j = (int)floor((x-fx)/dx);
	i = (int)floor((z-fz)/dz);



	/*求速度梯度*/
	if(j>=0 && j<= xn && i>=0 && i<=zn)
	{
	
          if(nx<0 &&nz<0)
		  {
			  if(i-1 <0 || j-1 <0) return 0;
			  v1 = v[(i-1)*xn+j-1];v2 = v[(i-1)*xn+j];
			  v3 = v[i*xn+j];      v4 = v[i*xn+j-1];
		  }
		  if(nx>0 && nz<0)
		  {
			  if(i-1 <0 || j+1 >xn) return 0;
			  v1 = v[(i-1)*xn + j] ; v2 = v[(i-1)*xn+j+1];
			  v3 = v[i*xn +j+1] ;    v4 = v[i*xn +j];
		  }
		  if(nx >0 && nz>0)
		  {
			  if(i+1 > zn || j+1 > xn) return 0;
			  v1 = v[i*xn+j]; v2 = v[i*xn +j+1];
			  v3 = v[(i+1)*xn+j+1]; v4 = v[(i+1)*xn +j];
		  }
          if(nx <0 && nz>0)
		  {
			  if(i+1 >zn || j-1<0) return 0;
			  v1 = v[i*xn +j-1]; v2 = v[i*xn +j];
			  v3 = v[(i+1)*xn+j];v4 = v[(i+1)*xn+j-1];
		  }
		  
	      l = (float)sqrt(dx*dx+dz*dz);
		  dv1 = (v3-v1)/l;dv1x = dv1*dx/l;dv1z = dv1*dz/l;
		  dv2 = (v4-v2)/l;dv2x =-dv2*dx/l;dv2z = dv2*dz/l;
		  *DVX = (float)0.5*(dv1x+dv2x);
		  *DVZ = (float)0.5*(dv1z+dv2z);
		 
	}else return 0;
    dvx = *DVX;dvz =*DVZ;
	if(j>= 0 && j<= xn && i >=0 && i <= zn)
	{
		if(fmod((x-fx),dx) == 0 && fmod((z-fz),dz) == 0)
			*vp = v[i*xn+j];
		else
			*vp = v[i*xn+j]+(float)fmod((x-fx),dx)*dvx+(float)fmod((z-fz),dz)*dvz;
        return 1;
	}else return 0;

}



float raytrace (float *X,float *Z,float *NX,float *NZ,float *T,float v,float dvx,float dvz,float dx,float dz,float fx,float fz,
			   struct dotinform*p,int *k,float xend,float zend)
{
	/*變量定義*/
   int i,j,flag;
   int  jj,ks;
   float x,z,nx,nz,t;
   jie so,s[8],ends[8];
   float deltx,zbefore,zafter;
   float xa,xb,xc1,xc2,delt,zs,za,zb,zc1,zc2,delt1,delt2,xs;
   float deltz,xbefore,xafter;
   float deltk;
  
   x = *X;
   z = *Z;
   nx = *NX;
   nz = *NZ;
   t  = *T;
  
   for(i=0;i<8;i++)
   {
	   s[i].ch= 'a';s[i].s = -1;s[i].val = 0;
	   ends[i].ch= 'a';ends[i].s= -1;ends[i].val=0;
   }
 
   so.ch = 'a' ; so.s = -1;so.val= 0;
  
  
   if(fmod(x,dx) == 0) {flag = 1;}   /*flag=1,射線從Z軸入射*/
   if(fmod(z,dz) == 0) {flag = 2;}
   if(flag != 1 && flag != 2) return 0; 
  
   if(v <= 0) return 0;
   i = 0;
   j = 0;
   ks= 0;   
   switch (flag){
   case 1:
	   if(nx >0)
	   {
		   deltx = dx;
	   }else{
		   deltx = -dx;
	   }
    
	   xa = (dvx*nx*nx+dvz*nx*nz-dvx)/(2*v);
	   xb = nx;
	   xc1= 0;
	   xc2= -deltx;
	   delt = xb*xb - 4*xa*xc2;
	   if(xa != 0)
	   {
	     s[i].s = -(xb/xa);s[i].ch = 'x';s[i].val = x;
	     i++;
	     if (delt >= 0)
		 {
		   s[i].s = (-xb+(float)sqrt(delt))/(2*xa);s[i].ch = 'x';s[i].val = (x+deltx);
		   i++;
           s[i].s = (-xb-(float)sqrt(delt))/(2*xa);s[i].ch = 'x';s[i].val = (x+deltx);
		   i++;
		 }
	   }else{
         s[i].s = -(xc2/xb);s[i].ch = 'x';s[i].val = (x+deltx); i++;
        
	   }
      
	   if(nz > 0)
	   {
		   zbefore = z - (float)fmod(z,dz);
	       zafter  = zbefore + dz;
	   }else{ 
		      if ( fmod(z,dz) == 0.0)
			  {
				  zafter = z ;
	              zbefore  = zafter - dz;
			  }else{
                 zbefore = z - (float)fmod(z,dz);
	             zafter  = zbefore + dz;
			  }
	   }
	   for(j = 0; j<i ;j++)
	   {
		   if (s[j].s >0) 
		   {
			   zs = s[j].s*s[j].s*(dvz*nz*nz+dvx*nx*nz-dvz)/(2*v)+nz*s[j].s +z;
		       if (zbefore <= zs && zs <= zafter) 
			   {
				   ends[ks].s = s[j].s;ends[ks].ch = s[j].ch;ends[ks].val = s[j].val;ks++;
			   }
		   }
	   }
	   

	   za = (dvz*nz*nz+dvx*nx*nz-dvz)/(2*v);
	   zb = nz;
	   zc1= z-zbefore;
	   zc2= z-zafter;
	   delt1 = zb*zb-4*za*zc1;                   
	   delt2 = zb*zb-4*za*zc2;
	   if(za != 0)
	   {
	     if (delt1 >= 0)
		 {
		   s[i].s = (-zb+(float)sqrt(delt1))/(2*za);s[i].ch = 'z';s[i].val = zbefore;i++;
           s[i].s = (-zb-(float)sqrt(delt1))/(2*za);s[i].ch = 'z';s[i].val = zbefore;i++;
		 }
	     if (delt2 >= 0)
		 {
           s[i].s = (-zb+(float)sqrt(delt2))/(2*za);s[i].ch = 'z';s[i].val = zafter;i++;
           s[i].s = (-zb-(float)sqrt(delt2))/(2*za);s[i].ch = 'z';s[i].val = zafter;i++;       
		 }
	   }else{
              s[i].s = -(zc1/nz);s[i].ch = 'z';s[i].val = zbefore;i++;
              s[i].s = -(zc2/nz);s[i].ch = 'z';s[i].val = zafter;i++;
	   }

	   for(jj = j ;jj <i ;jj++)
	   {
		   if (s[jj].s >0) 
		   {
			   xs = s[jj].s*s[jj].s*(dvx*nx*nx+dvz*nx*nz-dvx)/(2*v)+nx*s[jj].s+x;
		       if ((x<=xs && xs<=x+deltx)||(x+deltx <= xs && xs <= x)) 
			   {
				   ends[ks].s = s[jj].s;ends[ks].ch = s[jj].ch;ends[ks].val = s[jj].val;ks++;
			   }
		   }
	   }
	   so.s = 999999;
	   for(jj = 0 ;jj<ks ; jj++)
	   {
		   if(ends[jj].s<=so.s)
		   {
			   so.s = ends[jj].s;so.ch = ends[jj].ch;so.val = ends[jj].val;
		   }
	   }
	   break;

   case 2:
       if (nz >0) 
	   {
		   deltz = dz;
	   }else{
		   deltz = -dz;                       
	   }
	  
	   za = (dvz*nz*nz+dvx*nx*nz-dvz)/(2*v);
	   zb = nz;
	   zc1= 0;
	   zc2= -deltz;
	   delt= zb*zb - 4*za*zc2;
	   if(za != 0)
	   {
	     s[i].s=-(zb/za);s[i].ch = 'z';s[i].val = z;
	     i++;
	
	     if (delt >= 0)
		 {
		   s[i].s = (-zb+(float)sqrt(delt))/(2*za);s[i].ch = 'z';s[i].val = (z+deltz);i++;
           s[i].s = (-zb-(float)sqrt(delt))/(2*za);s[i].ch = 'z';s[i].val = (z+deltz);i++;
		 }
	   }else{
              s[i].s = -(zc2/zb);s[i].ch = 'z';s[i].val = (z+deltz);
	          i++;
	   }
	   if(x>0)
	   {
	     if(nx >=0)
		 {
	       xbefore = x - (float)fmod(x,dx);
	       xafter  = xbefore + dx;
		 }else{
		      if (fmod(x,dx) == 0.0)
			  { 
				  xafter = x;
				  xbefore = xafter -dx;
			  }else{ 
				   xbefore = x - (float)fmod(x,dx);
	               xafter  = xbefore + dx;
			  }
		 }
	   }else{
           if(nx <=0)
		   {
	         xafter = x - (float)fmod(x,dx);
	         xbefore  = xafter - dx;
		   }else{
		      if (fmod(x,dx) == 0.0)
			  { 
				  xbefore = x;
				  xafter = xbefore +dx;
			  }else{ 
				   xafter = x - (float)fmod(x,dx);
	               xbefore  = xafter - dx;
			  }
		   }
	   }

	   for (j=0;j<i;j++)
	   {
		   if(s[j].s>0) 
		   {
			   xs = s[j].s*s[j].s*(dvx*nx*nx+dvz*nx*nz-dvx)/(2*v)+nx*s[j].s+x;
		       if(xbefore<= xs&& xs<= xafter) 
			   {
			     ends[ks].s = s[j].s;ends[ks].ch = s[j].ch;ends[ks].val = s[j].val;ks++;
			   }
		   }
	   }
	  

	   xa = (dvx*nx*nx+dvz*nx*nz-dvx)/(2*v);
	   xb = nx;
	   xc1=x-xbefore;
	   xc2=x-xafter;
	   delt1 = xb*xb-4*xa*xc1;
	   delt2 = xb*xb-4*xa*xc2;
	   
	   if(xa != 0)
	   {
	     if(delt1>=0)
		 {
		   s[i].s = (-xb+(float)sqrt(delt1))/(2*xa);s[i].ch = 'x';s[i].val = xbefore;i++;
           s[i].s = (-xb-(float)sqrt(delt1))/(2*xa);s[i].ch = 'x';s[i].val = xbefore;i++;
		 }

         if(delt2>=0)
		 {
		   s[i].s = (-xb+(float)sqrt(delt2))/(2*xa);s[i].ch = 'x';s[i].val = xafter;i++;
           s[i].s = (-xb-(float)sqrt(delt2))/(2*xa);s[i].ch = 'x';s[i].val = xafter;i++;
		 }
	   }else{
           s[i].s = -(xc1/nx);s[i].ch = 'x';s[i].val = xbefore;i++;
           s[i].s = -(xc2/nx);s[i].ch = 'x';s[i].val = xafter;i++;
	   }
 
	   for(jj = j;jj<i;jj++)
	   {
         
		   if(s[jj].s>0) 
		   {
			   zs = s[jj].s*s[jj].s*(dvz*nz*nz+dvx*nx*nz-dvz)/(2*v)+nz*s[jj].s+z;
		       if((z<= zs &&zs <= z+deltz)||(z+deltz <= zs && zs<=z))
			   {
			     ends[ks].s = s[jj].s;ends[ks].ch = s[jj].ch;ends[ks].val = s[jj].val;ks++;
			   
			   }
		   }
	   }
      
	   
	   so.s=999999;
	   for(jj=0;jj<ks;jj++)
	   {
		   if(ends[jj].s<=so.s) 
		   {  
			   so.s = ends[jj].s;so.ch = ends[jj].ch;so.val = ends[jj].val;
		   }
	   }
	   break;
	}

	if (so.s != 999999)
	{
	   
		if(so.ch == 'x')
		{ 
			*X = so.val;
            *Z = so.s*so.s*(dvz*nz*nz+dvx*nx*nz-dvz)/(2*v)+nz*so.s+z;
		}
		if(so.ch == 'z')
		{
            *X = so.s*so.s*(dvx*nx*nx+dvz*nx*nz-dvx)/(2*v)+nx*so.s+x;
		    *Z = so.val;
		}
		*T = so.s*(1-(dvx*nx+dvz*nz)*so.s/(2*v))/v+t;
        *NX = nx+(dvx*nx*nx+dvz*nz*nx-dvx)*so.s/v;
        *NZ = nz+(dvx*nx*nz+dvz*nz*nz-dvz)*so.s/v;
     	deltk = ((*NX)*(*NX)+(*NZ)*(*NZ)-1)/2;
		*NX = (*NX)/(1+deltk);
		*NZ = (*NZ)/(1+deltk);
	
/*		if(fabs((*NX)*(*NX)) >=0.8 && fabs((*NX)*(*NX)) <1)
		{
            if(*Z >z) *NZ = (float)sqrt(1-(*NX)*(*NX));
			else *NZ = (float)-sqrt(1-(*NX)*(*NX));
		}
        if(fabs((*NZ)*(*NZ)) >=0.8 && fabs((*NZ)*(*NZ)) <1)
		{
           	if(*X >x) *NX = (float)sqrt(1-(*NZ)*(*NZ));
			else *NX = (float)-sqrt(1-(*NZ)*(*NZ));
		}
*/
	    if(fabs(*NX)>1)
		{
			if(*X >x) *NX = (float)sqrt(1-(*NZ)*(*NZ));
			*NX = (float)-sqrt(1-(*NZ)*(*NZ));
		}
		if(fabs(*NZ)>1)
		{
            if(*Z >z) *NZ = (float)sqrt(1-(*NX)*(*NX));
			*NZ = (float)-sqrt(1-(*NX)*(*NX));
		} 
	
		if(*X == x)  return 0;
		if ((fx<=*X)&&(*X<= xend)&&(fz<=*Z)&&(*Z<=zend))
		{
		   p[*k].x= *X;
		   p[*k].z= *Z;
		   p[*k].nx=*NX;
		   p[*k].nz=*NZ;
		   p[*k].tsum=*T;
		   (*k)++;
		}
	   if(*X <= fx || *X >= xend || *Z<= fz || *Z>= zend) return -1;
	   return 1 ;
	}
	else {
		
		return 0;
		
	}
   
}