sutetraray.c

su 的源代码库

			************************************************/
			
                        if (check_tris==PLEASE) {
			      dist[0]=(x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)
                                     +(z0-z1)*(z0-z1);
		      	      dist[1]=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)
                                     +(z2-z1)*(z2-z1);
			      dist[2]=(x2-x0)*(x2-x0)+(y2-y0)*(y2-y0)
                                     +(z2-z0)*(z2-z0);

		       	      distmax=dist[0];
			      iedgemax=0;
			      for (iedge=1;iedge<3;iedge++) {
			       	    if (dist[iedge]>distmax) {
				          iedgemax=iedge;
				      	  distmax=dist[iedge];
				    }
			      }

			      /* this triangle is ok. all edges within limits,
			      nothing needs to be done with this triangle */
			      if (distmax<edgemax*edgemax) 
			       	    continue; 

			      /* this triangle is too big. needs to be splitted 
			      into two triangles */
			      if (iedgemax==0) {	/* vertex 1 & 2 */
				    xvertex0=x0;
				    xvertex1=x1;
				    yvertex0=y0;
				    yvertex1=y1;
				    zvertex0=z0;
				    zvertex1=z1;
				    tripi0=trip->i1;
				    tripi1=trip->i2;
			      } else if (iedgemax==1) {/* vertex 2 & 3 */
			     	    xvertex0=x1;
				    xvertex1=x2;
				    yvertex0=y1;
				    yvertex1=y2;
				    zvertex0=z1;
				    zvertex1=z2;
				    tripi0=trip->i2;
				    tripi1=trip->i3;
			      } else {		/* vertex 1 & 3 */
			       	    xvertex0=x0;
				    xvertex1=x2;
				    yvertex0=y0;
				    yvertex1=y2;
				    zvertex0=z0;
				    zvertex1=z2;
				    tripi0=trip->i1;
				    tripi1=trip->i3;
			      }

			      if (ray[tripi0].shootupward!=ray[tripi1].shootupward)
                                    continue; /* better not add one here */

			      /***********************************************
			      Add one more ray: nrays: Interpolation of ray 
                              data, whatever you need, could be traveltime; 
                              amplitude; phase, and even more.
	      		      ***********************************************/
			      wf[ithis+nrays].x[0]=(xvertex0+xvertex1)/2.0;
			      wf[ithis+nrays].x[1]=(yvertex0+yvertex1)/2.0;
			      wf[ithis+nrays].x[2]=(zvertex0+zvertex1)/2.0;

			      /* Needs to determine which tetra it is in */
                              fprintf(jpfp,"Add one more ray\n");
		      	      ray[nrays].itetra=in_which_tetra(
                                    xmin,
                                    ymin,
                                    xmax,
                                    ymax,
                                    zmax,
                                    wf[ithis+nrays].x, 
				    point,
                                    tetra,
				    ntetra);

                              if (ray[nrays].itetra<0) {
				    ray[nrays].icode=SPOILED;
                                    fprintf(jpfp,"attempt to add this ray but failed\n");
			      }

			      ray[nrays].shootupward=ray[tripi0].shootupward;

			      #ifdef DEBUG
			      fprintf(jpfp,"add a ray (nray=%d) in tetra=%d\n",
			            nrays,ray[nrays].itetra);
			      fprintf(jpfp,"at x=%f %f %f\n",wf[ithis+nrays].x[0],
			            wf[ithis+nrays].x[1],wf[ithis+nrays].x[2]);

			      fprintf(jpfp,"this tetra: ip=%d %d %d %d\n",
			            tetra[ray[nrays].itetra].ip[0],
			            tetra[ray[nrays].itetra].ip[1],
			            tetra[ray[nrays].itetra].ip[2],
			            tetra[ray[nrays].itetra].ip[3]);
			      fprintf(jpfp,"ray[nrays].shootupward=%d\n",
				    ray[nrays].shootupward);
                              #endif

                              ray[nrays].r=(ray[tripi0].r+ray[tripi1].r)/2.0;
                              ray[nrays].tseg[0]=itwf*dtwf;
                              ray[nrays].ttotal=itwf*dtwf;

                              #ifdef DEBUG
			      fprintf(jpfp,"itwf=%d ray[%d].xseg=%f %f %f\n",
			            itwf,nrays,ray[nrays].xseg[0][0],
			            ray[nrays].xseg[0][1],ray[nrays].xseg[0][2]);
			      #endif

			      wf[ithis+nrays].p[0]=(wf[ithis+tripi0].p[0]+
			       	    wf[ithis+tripi1].p[0])/2.0;
			      wf[ithis+nrays].p[1]=(wf[ithis+tripi0].p[1]+
			       	    wf[ithis+tripi1].p[1])/2.0;
			      wf[ithis+nrays].p[2]=(wf[ithis+tripi0].p[2]+
			       	    wf[ithis+tripi1].p[2])/2.0;

			      s=s_intp(
				    point,              /* control points */
				    &tetra[initetra],   /* this tetra */
				    wf[ithis+nrays].x); /* the current position */

                              s_normalize(wf[ithis+nrays].p,s);
			      
                              r=(ray[tripi0].r+ray[tripi1].r)/2.0;
                              d2=(xvertex0-wf[ithis+nrays].x[0])*
                                 (xvertex0-wf[ithis+nrays].x[0])+
				 (yvertex0-wf[ithis+nrays].x[1])*
                                 (yvertex0-wf[ithis+nrays].x[1])+
				 (zvertex0-wf[ithis+nrays].x[2])*
                                 (zvertex0-wf[ithis+nrays].x[2]);
                              d=sqrt(d2);
                              sintheta=d/MAX(0.2,r);
                              sintheta=MIN(1.0,sintheta);
                              costheta=sqrt(1.0-sintheta*sintheta);
                              adjust=r-r*costheta;

			      #ifdef DEBUG
			      fprintf(jpfp,"adjust=%f,r=%f\n",adjust,r);
			      fprintf(jpfp,"p1,p2,p3=%f %f %f\n",
			      wf[ithis+nrays].p[0],
			      wf[ithis+nrays].p[1],
			      wf[ithis+nrays].p[2]);
			      #endif

                              wf[ithis+nrays].x[0]+=adjust*wf[ithis+nrays].p[0]/sqrt(s);
                              wf[ithis+nrays].x[1]+=adjust*wf[ithis+nrays].p[1]/sqrt(s);
			      wf[ithis+nrays].x[2]+=adjust*wf[ithis+nrays].p[2]/sqrt(s);

                              ray[nrays].xseg[0][0]=wf[ithis+nrays].x[0];
			      ray[nrays].xseg[0][1]=wf[ithis+nrays].x[1];
			      ray[nrays].xseg[0][2]=wf[ithis+nrays].x[2];

			      /* Needs to determine which tetra it is in */
		      	      ray[nrays].itetra=in_which_tetra(
                                    xmin,
                                    ymin,
                                    xmax,
                                    ymax,
                                    zmax,
                                    wf[ithis+nrays].x, 
				    point,
                                    tetra,
				    ntetra);

                              if (ray[nrays].itetra<0) {
				    ray[nrays].icode=SPOILED;
                                    fprintf(jpfp,"tried to adjust this ray, failed\n");
			      }

                              /********************************************
			      add 2 triangles.  Since any one edge is shared
                              by 2 triangles, we need to add one triangle on
                              one side, and the other one on the other side.
                              i.e. add one triangle to trip, and one to
                              trip0
                              ********************************************/

			      i1=trip->i1;
			      i2=trip->i2;
			      i3=trip->i3;

                              ii1=ii2=ii3=0;
                              iedgemax0=999;
                              /* find the next neighboring triangle */
		              for (trip0=headp;trip0!=(struct TRI*)NULL;
		     	            trip0=trip0->next) {
                                    if (trip0==trip) continue;
                                    ii1=trip0->i1;
                                    ii2=trip0->i2;
                                    ii3=trip0->i3;

                                    if ((ii1==i1 || ii1==i2 || ii1==i3) &&
                                        (ii2==i1 || ii2==i2 || ii2==i3)) {
                                          iedgemax0=0;
                                          break;
                                    } else if 
                                        ((ii2==i1 || ii2==i2 || ii2==i3) &&
                                         (ii3==i1 || ii3==i2 || ii3==i3)) {
                                          iedgemax0=1;
                                          break;
				    } else if 
                                        ((ii1==i1 || ii1==i2 || ii1==i3) &&
                                         (ii3==i1 || ii3==i2 || ii3==i3)) {
                                          iedgemax0=2;
                                          break;
                                    }
                              }

			      #ifdef DEBUG
                              fprintf(jpfp,
                              "ii1=%d ii2=%d ii3=%d i1=%d i2=%d i3=%d iedgemax0=%d\n",
                              ii1,ii2,ii3,i1,i2,i3,iedgemax0);
			      #endif

			      temtrip=trip->next;
			      if (iedgemax==0) {	/* vertex 1 & 2 */
			       	    trip->i2=nrays;
			            trip->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			            trip=trip->next;
				    trip->i1=nrays;
			            trip->i2=i2;
			      	    trip->i3=i3;
			      } else if (iedgemax==1) {/* vertex 2 & 3 */
			       	    trip->i3=nrays;
			            trip->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			            trip=trip->next;
			      	    trip->i1=i1;
			            trip->i2=nrays;
			       	    trip->i3=i3;
			      } else {		/* vertex 1 & 3 */
			       	    trip->i3=nrays;
			            trip->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			       	    trip=trip->next;
			            trip->i1=nrays;
			      	    trip->i2=i2;
			            trip->i3=i3;
			      }
	
			      trip->next=temtrip;
			      ntris++;

                              if (trip0==(struct TRI*)NULL || iedgemax0==999 ||
                                  ii1==ii2 || ii1==ii3 || ii2==ii3) {
                                    nrays++;
				    continue;
			      }

			      temtrip=trip0->next;
			      if (iedgemax0==0) {	/* vertex 1 & 2 */
			       	    trip0->i2=nrays;
			            trip0->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			            trip0=trip0->next;
				    trip0->i1=nrays;
			            trip0->i2=ii2;
			      	    trip0->i3=ii3;
			      } else if (iedgemax0==1) {/* vertex 2 & 3 */
			       	    trip0->i3=nrays;
			            trip0->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			            trip0=trip0->next;
			      	    trip0->i1=ii1;
			            trip0->i2=nrays;
			       	    trip0->i3=ii3;
			      } else {		/* vertex 1 & 3 */
			       	    trip0->i3=nrays;
			            trip0->next=(struct TRI *)
			       		  malloc(sizeof(struct TRI));
			       	    trip0=trip0->next;
			            trip0->i1=nrays;
			      	    trip0->i2=ii2;
			            trip0->i3=ii3;
			      }
	
			      trip0->next=temtrip;
			      ntris++;
                              nrays++;
			}
		  } /* end if check_tris */

		  /************************************************
	          Write the current wavefront to the wffile. This 
	      	  can be visualized by viewer3
		  ************************************************/
		  if (wffp!=NULL && itwf>=ifwf2dump && itwf<ifwf2dump+nwf2dump) {

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

			      if (wf[ithis+trip->i1].x[0]>0.0 &&
			          wf[ithis+trip->i1].x[1]>0.0 &&
			          wf[ithis+trip->i1].x[2]>0.0 &&
			          wf[ithis+trip->i2].x[0]>0.0 &&
			          wf[ithis+trip->i2].x[1]>0.0 &&
			          wf[ithis+trip->i2].x[2]>0.0 &&
			          wf[ithis+trip->i3].x[0]>0.0 &&
			          wf[ithis+trip->i3].x[1]>0.0 &&
			          wf[ithis+trip->i3].x[2]) ntr++;
                              else {
                                    ntr=0;
                                    break;
                              }
			}

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

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

			            fprintf(wffp,"%f %f %f %f %f %f %f %f %f\n",
                                          wf[ithis+trip->i1].x[0],
				          wf[ithis+trip->i1].x[1],
				          wf[ithis+trip->i1].x[2],
				          wf[ithis+trip->i2].x[0],
				          wf[ithis+trip->i2].x[1],
				          wf[ithis+trip->i2].x[2],
				          wf[ithis+trip->i3].x[0],
				          wf[ithis+trip->i3].x[1],
				          wf[ithis+trip->i3].x[2]);
                        }
                  }

		  if (ntris<2 || nrays<4) break;
            } /* end if it!=0 */

            if (itwf>=ntwf-1) break;

	    #ifdef DEBUG
            fprintf(jpfp,"Check done\n");
            #endif

            done=ALREADY;

	    for (iray=0;iray<nrays;iray++) {

	       	  if (ray[iray].icode==SPOILED || 
                        ray[iray].icode==HIT_SURFACE) continue;

                  done=NOT_YET;

      		  ray[iray].x[0]=wf[ithis+iray].x[0];
		  ray[iray].x[1]=wf[ithis+iray].x[1];
		  ray[iray].x[2]=wf[ithis+iray].x[2];

		  ray[iray].t=0.0;
		  ray[iray].p[0]=wf[ithis+iray].p[0];
		  ray[iray].p[1]=wf[ithis+iray].p[1];
		  ray[iray].p[2]=wf[ithis+iray].p[2];

		  /******************************************
		  Check if the current ray point source is out
		  of bounds
		  ******************************************/	

		  #ifdef DEBUG
                  fprintf(jpfp,"\nicode(iray=%d)=%d\n",
                        iray,ray[iray].icode);
		  #endif

		  while (HIT_FACET_FIRST==
                        (result=trace1tetra(
			maxnsegs,
		        xmin,ymin,xmax,ymax,zmax,
			ray+iray,/* this ray */
                        npoint,  /* number of points */
                        point,   /* control points */
                        nfacet,  /* number of facets */