sutetraray.c

su 的源代码库

      struct TRI head;	        /* first triangle */
      struct TRI *headp=NULL;	/* pointer to the first usable triangle */
      struct TRI *trip;	        /* pointer to tris */
      struct TRI *trip0;        /* runing pointer */
      struct TRI *itrip;	/* index for trip */
      struct TRI *temtrip;	/* temporary pointer to triangles */

      int nrays=0;	        /* number of rays */
      int ngoodrays;
      int ntris=0;	        /* number of triangles */

      float dist[3];

      int ixgd,iygd;
      float x0,x1,x2;
      float y0,y1,y2;
      float z0,z1,z2;

      float xvertex0,xvertex1;
      float yvertex0,yvertex1;
      float zvertex0,zvertex1;

      int initreg;
      int ntr;
      int initetra;
      int tripi0,tripi1;

      float **ttable;        /* traveltime table */
      float **ctable;
      float **rtable;
      enum CHECK_TRIS check_tris=PLEASE;
      enum RESULT result;
      enum DONE done;
      struct RAY *ray;

      float temp,sqrts,tsum;

      int iseg;
      float adjust;     /* when interpolating, ajust this amount */
      float r,d,d2,sintheta,costheta;

      ray=(struct RAY *)
	    malloc(maxnrays*sizeof(struct RAY));

      fprintf(jpfp,"source at %f %f %f\n",source[0],source[1],source[2]);

      ttable=ealloc2float(nxgd,nygd);
      ctable=ealloc2float(nxgd,nygd);
      rtable=ealloc2float(nxgd,nygd);

      for (ixgd=0;ixgd<nxgd;ixgd++)
            for (iygd=0;iygd<nygd;iygd++)
                   ttable[iygd][ixgd]=TraveltimeInfinity;

      memset((void *)ctable[0],0,
            nxgd*nygd*sizeof(float));
      memset((void *)rtable[0],0,
            nxgd*nygd*sizeof(float));

      /* the ray information is bookkept on horizons */
      wf=(struct WF *) 
	    malloc(sizeof(struct WF)*maxnrays*ntwf);

      initetra=in_which_tetra(
            xmin,
            ymin,
            xmax,
            ymax,
            zmax,
	    source,    /* source position */
	    point,     /* control points */
            tetra,     /* tetra array */
            ntetra);   /* number of tetra */

      if (initetra<0) err("Source out of model");
      fprintf(jpfp,"source in tetra %d\n",initetra);

      initreg=tetra[initetra].ireg;
      fprintf(jpfp,"initreg=%d\n",initreg);
	
      /*********************************************
      Initialize the ray position
      to be the source position.
      ********************************************* */
      for (iray=0;iray<maxnrays;iray++) {
      	    wf[iray].x[0]=source[0];
	    wf[iray].x[1]=source[1];
	    wf[iray].x[2]=source[2];
            ray[iray].xseg=ealloc2float(3,maxnsegs);
            ray[iray].tseg=ealloc1float(maxnsegs);
            for (iseg=0;iseg<maxnsegs;iseg++)
                  ray[iray].tseg[iseg]=0.0;
	    ray[iray].icode=NORMAL;
            ray[iray].shootupward=shootupward;
	    ray[iray].itetra=initetra;
            ray[iray].lastfacet=-1;
            ray[iray].ionce=0;
            ray[iray].iref=0;
	    ray[iray].nseg=0;
            ray[iray].r=0.0;
            ray[iray].ttotal=0.0;
      }

      for (iray=0;iray<nazimuth*ntakeoff+1;iray++) {
            wf[iray].x[0]=source[0];
	    wf[iray].x[1]=source[1];
	    wf[iray].x[2]=source[2];
	    ray[iray].xseg[0][0]=source[0];
	    ray[iray].xseg[0][1]=source[1];
            ray[iray].xseg[0][2]=source[2];
      }

      for (itwf=1;itwf<ntwf-1;itwf++) {
	     for (iray=0;iray<maxnrays;iray++) {
		   ithis=itwf*maxnrays+iray;
		   wf[ithis].x[0]=OUT;
		   wf[ithis].x[1]=OUT;
		   wf[ithis].x[2]=OUT;
	     }
      }

      azimuth=2.0*3.1415926;
      dazimuth=azimuth/nazimuth;
      dtakeoff=takeoff/ntakeoff;

      /*********************************************************
      The maximum number of rays must not be more than maxnrays
      *********************************************************/
      if (ntakeoff*nazimuth>=maxnrays) 
      	    err("too many rays: (ntakeoff*nazimuth=)%d>=(maxnrays=%d",
		  ntakeoff*nazimuth,maxnrays);

      done=NOT_YET;

      s=s_intp(
	    point,             /* control points */
            &tetra[initetra],  /* this tetra */
            source);           /* the current position */

      #ifdef DUBUG
      fprintf(jpfp,"sloth at the source=%f\n",s);
      #endif

      sqrts=sqrt(s);

      for (itwf=0;itwf<ntwf;itwf++) {
	    ithis=itwf*maxnrays;
	    inext=(itwf+1)*maxnrays;

            if (done==ALREADY) break;

	    if (itwf==0) {
		  /*****************************************
		  Pointing to the first triangle
		  **************************************** */
		  headp=&head;
		  /********************************************
		  There are nazimuth*ntakeoff+1 rays starting for each shot
		  ********************************************/
		  for (iray=0;iray<=nazimuth*ntakeoff;iray++) {
			/************************************
	    	      	Ray 0 is the center; takeoff is the 
		       	take-off angle; azimuth is the azimuth
		       	************************************/
		       	if (iray==0) {
                              itakeoff=0;
			      sintakeoff=0.0;
                              costakeoff=1.0;
			} else {
                              itakeoff=1+(iray-1)/nazimuth;
                              temp=dtakeoff*itakeoff;
	       		      sintakeoff=sin(temp);
	       	              costakeoff=cos(temp);
                        }

		       	/***************************************
		       	Index for azimuth
		       	****************************************/
		       	if (iray==0) {
                              iazimuth=0;
			      sinazimuth=0.0;
                              cosazimuth=1.0;
		      	} else {	
                              iazimuth=(iray-1)%nazimuth;
	  		      temp=(iazimuth+0.5*(itakeoff-1))*dazimuth+0.1;
		              sinazimuth=sin(temp);
		       	      cosazimuth=cos(temp);
                        }

		       	/***************************************
		       	the equation of the slowness vector is
		       	px^2 + py^2 + pz^2 = s,
		       	there are two freedoms.
		       	***************************************/
		       	wf[iray].p[0]=sintakeoff*cosazimuth*sqrts;
		       	wf[iray].p[1]=sintakeoff*sinazimuth*sqrts;
                        /* If the rays are shot upward, put -1 here; 
                        otherwise, put +1 here */
			if (shootupward==1)
                              wf[iray].p[2]=-costakeoff*sqrts;
                        else
		      	      wf[iray].p[2]=costakeoff*sqrts;
	       	  }

		  /******************************************
		  At first step, we have nazimuth*ntakeoff+1 rays,
		  which form (2*ntakeoff-1)*nazimuth triangles.
			
		  		4   3   2	
				 \ 3|2 /
				  \ | /
				 4 \|/ 1
			      5-----0-----1
				 5 /|\ nazimuth-1
				  / | \ 
				 / 6|  \
				6    nazimuth
		  ******************************************/
		  head.i1=0;
		  head.i2=1;
		  head.i3=2;

		  trip=headp;
		  for (iray=2;iray<=nazimuth;iray++) {
		       	trip->next=(struct TRI*) 
		      	      malloc(sizeof(struct TRI));
			trip=trip->next;
		      	trip->i1=0;
		       	trip->i2=iray;
		       	trip->i3=iray%nazimuth+1;
		  }

		  for (iazimuth=0;iazimuth<ntakeoff-1;iazimuth++) {
		       	for (itakeoff=0;itakeoff<nazimuth;itakeoff++) {
		       	      trip->next=(struct TRI*) 
			      malloc(sizeof(struct TRI));
			      trip=trip->next;
			      trip->i1=itakeoff+1+iazimuth*nazimuth;
			      trip->i2=(itakeoff+1)%nazimuth+1
                                    +iazimuth*nazimuth;
			      trip->i3=trip->i1+nazimuth;
			      trip->next=(struct TRI*) 
			       	    malloc(sizeof(struct TRI));
			      trip=trip->next;
				    trip->i1=(itakeoff+1)%nazimuth+1
                                           +iazimuth*nazimuth;
			      trip->i2=itakeoff+1+(iazimuth+1)*nazimuth;
			      trip->i3=trip->i1+nazimuth;
		        }
		  }

		  trip->next=(struct TRI*)NULL;
		  nrays=nazimuth*ntakeoff+1;
		  ntris=(2*ntakeoff-1)*nazimuth;
		  fprintf(jpfp,"Starting: ntrips=%d nrays=%d\n",
		       	ntris,nrays);
	     } else {
	          /****************************************************
		  Test if the head is spoiled. If the first ray is spoiled,
                  the next ray will take the place; All triangles connected
                  with these rays need to be removed;
		  ****************************************************/
		  while ((ray[headp->i1].icode==SPOILED ||
			  ray[headp->i2].icode==SPOILED ||
			  ray[headp->i3].icode==SPOILED) && 
			       ntris>1) {
			  headp=headp->next;
			  ntris--;
			  fprintf(jpfp,
			       	"head spoiled, cut off: ntris=%d\n",ntris);
		  }

		  if (ntris<2) break;
			

                  /****************************************************
                  Need to test if all triangles are still good. Loop over
                  all triangle, to see is all three rays of which the
                  triangle is composed to be NORMAL.
	      	  ****************************************************/
		  for (trip=headp;;) { 

			/* make sure that next triangle is not the end*/
			if (trip->next==(struct TRI*)NULL)
		       	      break;

                        /* if next triangle is good, then continue to
                        test the next */
			if (  ray[trip->next->i1].icode!=SPOILED && 
		       	      ray[trip->next->i2].icode!=SPOILED && 
			      ray[trip->next->i3].icode!=SPOILED) {
				    trip=trip->next;
				    continue; 
			 }

                         /************************************************
                         if the program execuates here, it means that
                         triangle trip is good; but trip->next is bad.
                         We need to cut off this triangle. What could 
                         happen to the one or more rays of this this triangle
                         for this misfortune to happen?
                         ************************************************/
		      	 if (ntris<1) break;

			 ntris--;

                         /***********************************************
                         trip is good; but trip->next is bad. We need to see
                         if trip0->next->next is also bad, or even further.
                         In these case, we need to cut off all the consecutive
                         bad triangles.
                         ************************************************/
		       	 for (itrip=trip->next;;itrip=itrip->next) {

		               /* bad to the very end */
		       	       if (itrip->next==(struct TRI*)NULL)
			       	     break;	
			       /* good. the bad chain has stopped */
		               if (ray[itrip->next->i1].icode!=SPOILED &&
				   ray[itrip->next->i2].icode!=SPOILED &&
				   ray[itrip->next->i3].icode!=SPOILED) 
                                         break;

			       ntris--;
                               free(itrip);
		       	 }

                         /* now: both trip and trip->next are good ones */
		       	 trip->next=itrip->next;
		   }
                   /* from headp to the end, there are ntris triangles
                   they are all good ones. */	

                  /* to dump the wavefront */
		  if (wffp!=NULL && itwf==ifwf2dump)
		        fprintf(wffp,"%d = nwf2dump\n",nwf2dump);

	       	  /****************************************************
		  Test if we need to add some triangles
		  ****************************************************/
		  for (trip=headp;trip!=(struct TRI*)NULL;
		     	trip=trip->next) {

		        if (nrays>=maxnrays) break;
		        if (ntris>=maxntris) {
		       	      fprintf(jpfp,
		       	"ntris >= maxntris: wavefront after this step is spoiled\n");
		              break;
		       	}

		       	if (ray[trip->i1].icode==SPOILED ||
		       	    ray[trip->i2].icode==SPOILED ||
		       	    ray[trip->i3].icode==SPOILED) 
			       err("Triangle wrong: Impossible! But that's life!"); 

                        /* ray one */
		        x0=wf[ithis+trip->i1].x[0];
		        y0=wf[ithis+trip->i1].x[1];
			z0=wf[ithis+trip->i1].x[2];
                        
			/* ray two */
			x1=wf[ithis+trip->i2].x[0];
			y1=wf[ithis+trip->i2].x[1];
			z1=wf[ithis+trip->i2].x[2];

			/* ray three */
			x2=wf[ithis+trip->i3].x[0];
			y2=wf[ithis+trip->i3].x[1];
		        z2=wf[ithis+trip->i3].x[2];

		        /************************************************
		        Write the current wavefront to the wffile. This 
			can be visualized by viewer3