sutetraray.c

su 的源代码库

                        facet,   /* facet information */
                        ntetra,  /* number of tetra */
                        tetra,   /* tetra information */
                        irefseq, /* reflector sequence */
                        nrefseq, /* number of reflectors */
		       	jpfp,
		       	dtwf)))	/* step size for t */
        	  {};

                  if (result==OUT_OF_MODEL) 
                        ray[iray].icode=SPOILED;

		  if (itwf<ntwf-1) {	
		       	wf[inext+iray].x[0]=ray[iray].x[0];
		       	wf[inext+iray].x[1]=ray[iray].x[1];
		       	wf[inext+iray].x[2]=ray[iray].x[2];
		       	wf[inext+iray].p[0]=ray[iray].p[0];
		       	wf[inext+iray].p[1]=ray[iray].p[1];
		       	wf[inext+iray].p[2]=ray[iray].p[2];
		  } else {
		        fprintf(jpfp,"ntwf is not enough\n");
		        enough_ntwf=WF_FALSE;
		  }

                  /* plot the rays wherever they go */
                  if (result==HIT_SURFACE_FIRST) {
                        check_tris=PLEASE_DO_NOT;
                        ray[iray].icode=HIT_SURFACE;

                        if (fabs(ray[iray].x[2])>0.0001) 
                              fprintf(jpfp,"x[2]!=0.0, =%f\n",ray[iray].x[2]);
                  }
	    }

            #ifdef DEBUG
            fprintf(jpfp,"itwf=%d done\n",itwf);
            #endif
      }/* next it */

      /******************************************************************
      Now, let's output the ray information for viewer3 to check if all
      the rays are correct. 
      ******************************************************************/
      if (rayfp!=NULL) {
            fprintf(rayfp,
                  "%d =Maximum number of segments\n",maxnsegs);

            ngoodrays=0;
            for (iray=0;iray<nrays;iray++)
                  if (ray[iray].icode!=SPOILED) ngoodrays++;

            fprintf(rayfp,
                  "%d =Number of rays\n",ngoodrays);

	    for (iray=0;iray<nrays;iray++) {
                  if (ray[iray].icode==SPOILED) continue;
                  fprintf(rayfp,"%d=nseg %f=ttotal\n",
                        ray[iray].nseg+1,ray[iray].ttotal);

                  tsum=0.0;
                  for (iseg=0;iseg<=ray[iray].nseg;iseg++) {
		        tsum+=ray[iray].tseg[iseg];
                        fprintf(rayfp,"%f %f %f %f %f\n",
                              ray[iray].xseg[iseg][0],
			      ray[iray].xseg[iseg][1],
			      ray[iray].xseg[iseg][2],
			      tsum,
			      ray[iray].tseg[iseg]);
                  }
            }
      }

      if (enough_ntwf==WF_FALSE) return 1;

      /* out the surface traveltime, only for visualization */
      if (sttfp!=NULL) {
	    for (trip=headp,ntr=0;trip!=(struct TRI*)NULL;
		  trip=trip->next) { 

	          if (ray[trip->i1].icode==SPOILED  ||
	    	      ray[trip->i2].icode==SPOILED  ||
	       	      ray[trip->i3].icode==SPOILED  ||
                      fabs(ray[trip->i1].x[2])>0.01 ||
		      fabs(ray[trip->i2].x[2])>0.01 ||
		      fabs(ray[trip->i3].x[2])>0.01) continue;

		  ntr++;
            }

            fprintf(stderr,"ntr=%d\n",ntr);

            fprintf(sttfp,"%d = ntris\n",ntr);

            if (ntr>0) {
                  for (trip=headp;trip!=(struct TRI*)NULL;
		       trip=trip->next) {

	                if (ray[trip->i1].icode==SPOILED  ||
	    	            ray[trip->i2].icode==SPOILED  ||
	       	            ray[trip->i3].icode==SPOILED  ||
                            fabs(ray[trip->i1].x[2])>0.01 ||
		            fabs(ray[trip->i2].x[2])>0.01 ||
		            fabs(ray[trip->i3].x[2])>0.01) continue;

	                fprintf(sttfp,"%f %f %f %f %f %f %f %f %f %f %f %f\n",
                              ray[trip->i1].x[0],
	                      ray[trip->i1].x[1],
			      ray[trip->i1].x[2],
			      ray[trip->i1].ttotal,
			      ray[trip->i2].x[0],
			      ray[trip->i2].x[1],
			      ray[trip->i2].x[2],
			      ray[trip->i2].ttotal,
			      ray[trip->i3].x[0],
			      ray[trip->i3].x[1],
			      ray[trip->i3].x[2],
			      ray[trip->i3].ttotal);
                  }
            }
            fclose(sttfp);
            sttfp=NULL;
      }

      /**************************************************
      To do the evaluation on the surface, for each 
      triangle, receivers inside this triangle will be
      evaluated. If a receiver is evaluated more than
      once, only the one with shortest traveltime
      is kept.
      *************************************************/
      for (trip=headp;trip!=NULL;trip=trip->next) {

	    if (ray[trip->i1].icode==SPOILED ||
	    	ray[trip->i2].icode==SPOILED ||
	       	ray[trip->i3].icode==SPOILED ) continue;

            tri_multi_intp(
                  &ray[trip->i1],
                  &ray[trip->i2],
                  &ray[trip->i3],
                  fxgd,
                  fygd,
                  dxgd,
                  dygd,
                  nxgd,
                  nygd,
		  ttable,
		  ctable,
		  rtable);
      }

      /******************************************************
      Now let us output the traveltime table to stdout, this
      will be used for sukdmig3d. Remember sx and sy are in 
      dxgd and dygd already.
      ******************************************************/
      for (iygd=iytf;iygd<iytf+nyt;iygd++) {
	    for (ixgd=ixtf;ixgd<ixtf+nxt;ixgd++) 
       		  tro.data[ixgd-ixtf]=MAX(0.0,ttable[iygd][ixgd]*1000.0);
            
       	    /* source left/front range */
       	    tro.gx=ixtf;
       	    tro.gy=iytf;

       	    /* source positions */
       	    tro.sx=(int)(source[0]*1000.0);
       	    tro.sy=(int)(source[1]*1000.0);
            tro.unscale=source[2];  /* this is a float */

       	    fputtr(stdout,&tro);

            if (crfp==NULL) continue;

            if (fwrite(&ctable[iygd][ixtf],sizeof(float),nxt,crfp)!=nxt)
                  err("Can not write to crfile");

            if (fwrite(&rtable[iygd][ixtf],sizeof(float),nxt,crfp)!=nxt)
                  err("Can not write to crfile");

      }

      fflush(stdout);
      fprintf(jpfp,"traveltimes for this shot written.\n");
      free2float(ttable);
	
      return 0;
}

float zhorz(float xx,float yy,float n10,float n20,
float n30,float x00,float y00,float z00)
/*****************************************************************
zhorz - linearly interpolation
*****************************************************************
Function prototype:
float zhorz(float xx,float yy,float n10,float n20,
float n30,float x00,float y00,float z00);
*****************************************************************
Inputs:
float xx 	x coordinate of the point where z coordinate will be returned
float yy        y coordinate of the point where z coordinate will be returned
float n10       x component of the normal to the plane
float n20       y component of the normal to the plane
float n30       z component of the normal to the plane
float x00       x coordinate of the point in the plane
float y00       y coordinate of the point in the plane
float z00       z coordinate of the point in the plane

Return:
the z coordinate of the point

Notes:
The equation of the plane is: 
	
(xx-x00)*n1ll+(yy-y00)*n2ll+(zz-z00)*n3ll=0

thus 

zz=z00-((xx-x00)*n1ll+(yy-y00)*n2ll)/n3ll;
*****************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
****************************************************************/
{
      if (fabs(n30)<0.001) return DistanceInfinity;
      return (z00-((xx-x00)*n10+(yy-y00)*n20)/n30);
}

float f(float sigma,float s,float p,float a)
/**********************************************************
f - from sigma to x
***********************************************************
Inputs:
float sigma	ray running parameter
float s         sloth
float p         ray parameter
float a         gradient of sloth

Return:
component of the ray position
Notices:
According to the formula (Cheveny 1987), the ray path satisfies:

x_i(sigma) = x_i(0) + p_i*sigma + 0.25*s_i*sigma^2

This function f will return the new cartesian coordinate:
***********************************************************
Author: CWP: Zhaobo Meng, Sept 1997
**********************************************************/
{
      return(s+p*sigma+0.25*a*sigma*sigma);
}

float cubicsolver(float d,float a,float b,float c)
/*****************************************************************
cubicsolver - solver for cubic equation
*****************************************************************
Function prototype:
float cubicsolver(float d,float a,float b,float c);
*****************************************************************
Inputs:
float d,a,b,c  coefficients of the cubic equation
return the solution for sigma
*****************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
******************************************************************/
{
      float sigma,fderi,fsigma;
      int MaxIterations=20;
      int iter;
      
      sigma=quadsolver(a,b,c);

      if (d==0.0) return sigma;

      if (sigma>0.5*SigmaInfinity) {
            sigma=1.0;
            MaxIterations=30;
      }

      for (iter=0;iter<MaxIterations;iter++) {
  	    /* fderi=3dx^2+2ax+b */

	    fsigma=((d*sigma+a)*sigma+b)*sigma+c;

	    if (fabs(fsigma)<EPS) break;

	    fderi=(3.0*d*sigma+2.0*a)*sigma+b;

            if (fabs(fderi)<EPS) break;

	    sigma-=fsigma/fderi;
      }
      
      if (sigma<0.0) return 0.0;
      
      return sigma;
}

float quadsolver(float a,float b,float c)
/**************************************************************
quadsolver - solver for quadratic equation ax^2+bx+c=0
**************************************************************
Function prototype:
float quadsolver(float a,float b,float c);
**************************************************************
Inputs:
float a,b,c  coefficients of the quadratic equation
return the solution for sigma
**************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
*************************************************************/

{
      float qdelta;
      float sigma,sigma1,sigma2;

      if (c==0.0) {
            if (a==0.0) return SigmaInfinity;
            sigma=-b/a;
            if (sigma<0.0) return SigmaInfinity;
            return sigma;
      }

      if (a==0.0) {
            if (b==0.0) return SigmaInfinity; 
            sigma=-c/b;
            if (sigma<0.0) sigma=SigmaInfinity;
            return sigma;
      } else {
            qdelta=b*b-4*a*c;
            if (qdelta<-EPS) return SigmaInfinity; 
            qdelta=MAX(0.0,qdelta);
            qdelta=sqrt(qdelta);
            sigma1=(qdelta-b)/(2.0*a);
            sigma2=(-qdelta-b)/(2.0*a);
            if (sigma1<-EPS) sigma1=SigmaInfinity;
            if (sigma2<-EPS) sigma2=SigmaInfinity;
            sigma=MIN(sigma1,sigma2);
            return MAX(0.0,sigma);
      }
}

int in_which_tetra(float xmin,float ymin,float xmax,float ymax,
float zmax,float x[3],struct POINT *point,struct TETRA *tetra,int ntetra)
/*****************************************************************
in_which_tetra - x falls in which tetra?
******************************************************************
Function prototype:
int in_which_tetra(float xmin,float ymin,float xmax,float ymax,
float zmax,float x[3],struct POINT *point,struct TETRA *tetra,int ntetra);
******************************************************************
inputs:
float xmin              x min
float ymin              y min              
float xmax              x max
float ymax              y max
float zmax              z max
float x[3]              coordinate of the point
struct POINT *point     control points
struct TETRA *tetra     tetra array
int ntetra              number of tetra
******************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
******************************************************************/
{
      int it;
      int ip0,ip1,ip2,ip3;
      float vv,vv0,vv1,vv2,vv3;
      float xxmin[4],xxmax[4];

      float vd