tetramod.c

su 的源代码库

		   (ix,iy+1) (ix+1,iy+1)	
			________
			|\     |
			| \  1 |
			|  \   |
			|   \  | 
			| 0  \ |
			|     \|
			--------  
		    (ix,iy) (ix+1,iy) 

	    **********************************************************/
	    for (ixhz=0;ixhz<nxhz-1;ixhz++) {
		   for (iyhz=0;iyhz<nyhz-1;iyhz++) {
		 	  tm_normal(xhz[ihz][iyhz][ixhz],xhz[ihz][iyhz][ixhz+1],
			        xhz[ihz][iyhz+1][ixhz],
			        yhz[ihz][iyhz][ixhz],yhz[ihz][iyhz][ixhz+1],
       		                yhz[ihz][iyhz+1][ixhz],
		                zhz[ihz][iyhz][ixhz],zhz[ihz][iyhz][ixhz+1],
       		          	zhz[ihz][iyhz+1][ixhz],n1ll[iyhz]+ixhz,
				n2ll[iyhz]+ixhz,n3ll[iyhz]+ixhz);

			  tm_normal(xhz[ihz][iyhz][ixhz+1],xhz[ihz][iyhz+1][ixhz+1],
				xhz[ihz][iyhz+1][ixhz],
			       	yhz[ihz][iyhz][ixhz+1],yhz[ihz][iyhz+1][ixhz+1],
				yhz[ihz][iyhz+1][ixhz],
			        zhz[ihz][iyhz][ixhz+1],zhz[ihz][iyhz+1][ixhz+1],
				zhz[ihz][iyhz+1][ixhz],n1ur[iyhz]+ixhz,
                                n2ur[iyhz]+ixhz,n3ur[iyhz]+ixhz);
	       	   }
       	    }

	    /****************************************
	    Using linear interpolate the normals:
	    Four corners first; four edges second;
	    interior last
	    ****************************************/
	    n1hz[ihz][0][0]=n1ll[0][0];
	    n2hz[ihz][0][0]=n2ll[0][0];
	    n3hz[ihz][0][0]=n3ll[0][0];
	    n1hz[ihz][0][nxhz-1]=n1ur[0][nxhz-2];
	    n2hz[ihz][0][nxhz-1]=n2ur[0][nxhz-2];
	    n3hz[ihz][0][nxhz-1]=n3ur[0][nxhz-2];
	    n1hz[ihz][nyhz-1][0]=n1ll[nyhz-2][0];
	    n2hz[ihz][nyhz-1][0]=n2ll[nyhz-2][0];
	    n3hz[ihz][nyhz-1][0]=n3ll[nyhz-2][0];
	    n1hz[ihz][nyhz-1][nxhz-1]=n1ur[nyhz-2][nxhz-2];
	    n2hz[ihz][nyhz-1][nxhz-1]=n2ur[nyhz-2][nxhz-2];
	    n3hz[ihz][nyhz-1][nxhz-1]=n3ur[nyhz-2][nxhz-2];
       
            for (ixhz=1;ixhz<nxhz-1;ixhz++) {
	       	  n1hz[ihz][0][ixhz]=n1ll[0][ixhz-1]
	      		+n1ur[0][ixhz-1]+n1ll[0][ixhz];
	       	  n2hz[ihz][0][ixhz]=n2ll[0][ixhz-1]
		        +n2ur[0][ixhz-1]+n2ll[0][ixhz];       
		  n3hz[ihz][0][ixhz]=n3ll[0][ixhz-1]
		        +n3ur[0][ixhz-1]+n1ll[0][ixhz];       

		  n1hz[ihz][nyhz-1][ixhz]=n1ur[nyhz-2][ixhz-1]
                        +n1ll[nyhz-2][ixhz]+n1ur[nyhz-2][ixhz]; 
		  n2hz[ihz][nyhz-1][ixhz]=n2ur[nyhz-2][ixhz-1]
                        +n2ll[nyhz-2][ixhz]+n2ur[nyhz-2][ixhz];  
		  n3hz[ihz][nyhz-1][ixhz]=n3ur[nyhz-2][ixhz-1]
                        +n3ll[nyhz-2][ixhz]+n3ur[nyhz-2][ixhz];  
	    }

	    for (iyhz=1;iyhz<nyhz-1;iyhz++) {
	       	  n1hz[ihz][iyhz][0]=n1ll[iyhz-1][0]
                        +n1ur[iyhz-1][0]+n1ll[iyhz][0];
	          n2hz[ihz][iyhz][0]=n2ll[iyhz-1][0]
                        +n2ur[iyhz-1][0]+n2ll[iyhz][0];
		  n3hz[ihz][iyhz][0]=n3ll[iyhz-1][0]
                        +n3ur[iyhz-1][0]+n3ll[iyhz][0];
                  n1hz[ihz][iyhz][nxhz-1]=n1ur[iyhz-1][nxhz-2]
                        +n1ll[iyhz][nxhz-2]+n1ur[iyhz][nxhz-2];
		  n2hz[ihz][iyhz][nxhz-1]=n2ur[iyhz-1][nxhz-2]
                        +n2ll[iyhz][nxhz-2]+n2ur[iyhz][nxhz-2];
	          n3hz[ihz][iyhz][nxhz-1]=n3ur[iyhz-1][nxhz-2]
                        +n3ll[iyhz][nxhz-2]+n3ur[iyhz][nxhz-2];
            }
		
            for (iyhz=1;iyhz<nyhz-1;iyhz++) {
	     	  for (ixhz=1;ixhz<nxhz-1;ixhz++) {
	       	        n1hz[ihz][iyhz][ixhz]
		              =n1ll[iyhz][ixhz-1]
			      +n1ur[iyhz][ixhz-1]
		     	      +n1ll[iyhz][ixhz]
	       		      +n1ur[iyhz-1][ixhz]
	       		      +n1ll[iyhz-1][ixhz]
	       		      +n1ur[iyhz-1][ixhz-1];
	       	        n2hz[ihz][iyhz][ixhz]
                              =n2ll[iyhz][ixhz-1]
                              +n2ur[iyhz][ixhz-1]
                              +n2ll[iyhz][ixhz]
                              +n2ur[iyhz-1][ixhz]
                              +n2ll[iyhz-1][ixhz]
                              +n2ur[iyhz-1][ixhz-1];
		        n3hz[ihz][iyhz][ixhz]
                              =n3ll[iyhz][ixhz-1]
                              +n3ur[iyhz][ixhz-1]
                              +n3ll[iyhz][ixhz]
                              +n3ur[iyhz-1][ixhz]
                              +n3ll[iyhz-1][ixhz]
                              +n3ur[iyhz-1][ixhz-1];
		 }
	  }

	  for (iyhz=0;iyhz<nyhz;iyhz++) 
		 for (ixhz=0;ixhz<nxhz;ixhz++) 
	            	tm_normalize(&n1hz[ihz][iyhz][ixhz],
		              &n2hz[ihz][iyhz][ixhz],
		       	      &n3hz[ihz][iyhz][ixhz]);
      }

      /******************************************
      Write model parameters to hzfile
      ******************************************/
      fwrite(&nhz,sizeof(int),1,hzfp);
      fwrite(&nxhz,sizeof(int),1,hzfp);
      fwrite(&nyhz,sizeof(int),1,hzfp);
      fwrite(&xmin,sizeof(float),1,hzfp);
      fwrite(&xmax,sizeof(float),1,hzfp);
      fwrite(&ymin,sizeof(float),1,hzfp);
      fwrite(&ymax,sizeof(float),1,hzfp);
      fwrite(&zmin,sizeof(float),1,hzfp);
      fwrite(&zmax,sizeof(float),1,hzfp);
      
      /**********************************************************
      Output the horizon information to file hzfile:
      **********************************************************/	
      for (ihz=0;ihz<nhz;ihz++) {
       	    if (fwrite(xhz[ihz][0],sizeof(float),nxhz*nyhz,hzfp)!=nxhz*nyhz)
		  err("Can not write xhz to hzfp");	

	    if (fwrite(yhz[ihz][0],sizeof(float),nxhz*nyhz,hzfp)!=nxhz*nyhz)
	          err("Can not write yhz to hzfp");

       	    if (fwrite(zhz[ihz][0],sizeof(float),nxhz*nyhz,hzfp)!=nxhz*nyhz)
		  err("Can not write zhz to hzfp");

	    if (fwrite(v0hz[ihz][0],sizeof(float),nxhz*nyhz,hzfp)!=
	          nxhz*nyhz)
	          err("Can not write v0hz to hzfp");

	    if (fwrite(v1hz[ihz][0],sizeof(float),nxhz*nyhz,hzfp)!=
	          nxhz*nyhz)
	          err("Can not write v1hz to hzfp");		
      }

      /**********************************************
      For even iyhz: Lower left   Upper right
                       2             1__2
                       |\            \  |
                       | \            \ |
                       |__\            \|
                       0   1            0

      For odd iyhz: Lower left    Upper right
                       1  2             2
                       |--/            /|
                       | /            / |
                       |/            /__|
                       0             0   1
      ***********************************************/
 
      ntrap=(nxhz-1)*(nyhz-1);         /*number of trapzoids per hz*/
      nhztrap=nhz*ntrap;               /*total number of trapzoids*/
      nhztrap1=(nhz+1)*ntrap;
      nfacetx=(nyhz-1)*nxhz;
      nhzfacetx=nhz*nfacetx;
      nfacety=nyhz*(nxhz-1);
      nhzfacety=nhz*nfacety;
      nhzpoint=nhz*nxhz*nyhz;

      if (verbose) {
	    warn("ntrap=%d nhztrap=%d nhztrap1=%d",ntrap,nhztrap,nhztrap1);
	    warn("nfacetx=%d nhzfacetx=%d",nfacetx,nhzfacetx);
	    warn("nfacety=%d nfacety=%d",nfacety,nhzfacety);
      }

      npoint=2*nhzpoint;
      point=(struct POINT *)malloc(npoint*sizeof(struct POINT));
      if (verbose) fprintf(stderr,"npoint=%d\n",npoint);

      nfacet=2*nhztrap1+             /*facet on horizons: if1 and if2*/
            6*nhztrap+               /*slant: if3 if4 if5 if6 if7 and if8*/
            2*nhzfacetx+             /*normal to x*/
            2*nhzfacety;             /*normal to y*/

      facet=(struct FACET *)malloc(nfacet*sizeof(struct FACET));
      if (verbose) fprintf(stderr,"nfacet=%d\n",nfacet);

      ntetra=6*nhztrap;
      tetra=(struct TETRA *)malloc(ntetra*sizeof(struct TETRA));
      if (verbose) fprintf(stderr,"ntetra=%d\n",ntetra);

      /*Assign the control point information*/
      for (ihz=0;ihz<nhz;ihz++) {
	    fprintf(stderr,"processing: ihz=%d\n",ihz);
            ihz0=MIN(nhz-1,ihz+1);
            for (iyhz=0;iyhz<nyhz;iyhz++) {
                  for (ixhz=0;ixhz<nxhz;ixhz++) {
                        ip=ihz*nxhz*nyhz+iyhz*nxhz+ixhz;
                        point[ip].x[0]=xhz[ihz][iyhz][ixhz];
		        point[ip].x[1]=yhz[ihz][iyhz][ixhz];
		        point[ip].x[2]=zhz[ihz][iyhz][ixhz];
		        point[ip].s=1.0/v0hz[ihz][iyhz][ixhz]/v0hz[ihz][iyhz][ixhz];
                        point[ip].n[0]=n1hz[ihz][iyhz][ixhz];
                        point[ip].n[1]=n2hz[ihz][iyhz][ixhz];
			point[ip].n[2]=n3hz[ihz][iyhz][ixhz];

                        ip+=nhzpoint;
                        point[ip].x[0]=xhz[ihz0][iyhz][ixhz];
		        point[ip].x[1]=yhz[ihz0][iyhz][ixhz];
		        point[ip].x[2]=zhz[ihz+1][iyhz][ixhz];
		        point[ip].s=1.0/v1hz[ihz][iyhz][ixhz]/v1hz[ihz][iyhz][ixhz];
                        if (ihz==nhz-1) {
			      point[ip].n[0]=0.0;
                              point[ip].n[1]=0.0;
			      point[ip].n[2]=1.0;
                        } else {
			      point[ip].n[0]=n1hz[ihz+1][iyhz][ixhz];
                              point[ip].n[1]=n2hz[ihz+1][iyhz][ixhz];
			      point[ip].n[2]=n3hz[ihz+1][iyhz][ixhz];
                        }
                  }
            }
      }

      for (ihz=0;ihz<nhz;ihz++) {
            ihz0=MIN(nhz-1,ihz+1);
            /*For each trapzoid*/
            for (iyhz=0;iyhz<nyhz-1;iyhz++) {
	          for (ixhz=0;ixhz<nxhz-1;ixhz++) {

                        if1=ihz*ntrap+iyhz*(nxhz-1)+ixhz; /*left facet on hz*/
                        iwh=nhztrap1;
                        if2=iwh+if1;                      /*right facet on hz*/
                        iwh=nhztrap1*2;
                        if3=iwh+if1;                      /*upper diagonal facet*/
                        iwh=nhztrap1*2+nhztrap;
                        if4=iwh+if1;                      /*lower diagonal facet*/
                        iwh=nhztrap1*2+nhztrap*2;
                        if5=iwh+if1;                      /*left slant1 facet*/
		        iwh=nhztrap1*2+nhztrap*3;
		        if6=iwh+if1;                      /*right slant1 facet*/
		        iwh=nhztrap1*2+nhztrap*4;
                        if7=iwh+if1;                      /*left slant2 facet*/
		        iwh=nhztrap1*2+nhztrap*5;
		        if8=iwh+if1;                      /*right slant2 facet*/
		        iwh=nhztrap1*2+nhztrap*6;
                        if9=iwh+ihz*nfacetx+iyhz*nxhz+ixhz;/*upper facet facing x*/
                        if10=nhzfacetx+if9;               /*lower facet facing x*/
		        iwh=nhztrap1*2+nhztrap*6+nhzfacetx*2;
                        if11=iwh+ihz*nfacety+iyhz*(nxhz-1)+ixhz;/*upper facet facing y*/
                        if12=nhzfacety+if11;              /*lower facet facing y*/

                        it1=ihz*ntrap+iyhz*(nxhz-1)+ixhz; /*left top tetra*/
                        it2=nhztrap+it1;                  /*left middle tetra*/
                        it3=nhztrap*2+it1;                /*left bottom tetra*/
                        it4=nhztrap*3+it1;                /*right top tetra*/
                        it5=nhztrap*4+it1;                /*right middle tetra*/
                        it6=nhztrap*5+it1;                /*right bottom tetra*/

                        i000=ihz*nxhz*nyhz+iyhz*nxhz+ixhz;/*upper 4 control point*/
                        i010=i000+nxhz;
                        i100=i000+1;
                        i110=i000+nxhz+1;

                        i001=nhzpoint+i000;    /*lower 4 control point*/
                        i011=nhzpoint+i010;
                        i101=nhzpoint+i100;
                        i111=nhzpoint+i110;
                           
		        tetra[it1].ireg=ihz;
		        tetra[it2].ireg=ihz;
		        tetra[it3].ireg=ihz;
		        tetra[it4].ireg=ihz;
		        tetra[it5].ireg=ihz;
		        tetra[it6].ireg=ihz;

                        /*normals to other facet will be calculated later*/
                        if (iyhz%2==0) {
                              tetra[it1].ip[0]=i000;     /*left top tetra*/
                              tetra[it1].ip[1]=i010;
                              tetra[it1].ip[2]=i100;
                              tetra[it1].ip[3]=i101;

                              tetra[it2].ip[0]=i000;     /*left middle tetra*/
                              tetra[it2].ip[1]=i010;
                              tetra[it2].ip[2]=i001;
                              tetra[it2].ip[3]=i101;

                              tetra[it3].ip[0]=i001;     /*left bottom tetra*/
                              tetra[it3].ip[1]=i011;
                              tetra[it3].ip[2]=i101;
                              tetra[it3].ip[3]=i010;

                              tetra[it4].ip[0]=i010;     /*right top tetra*/
                              tetra[it4].ip[1]=i110;
                              tetra[it4].ip[2]=i100;
                              tetra[it4].ip[3]=i101;

                              tetra[it5].ip[0]=i010;     /*right middle tetra*/
                              tetra[it5].ip[1]=i011;
                              tetra[it5].ip[2]=i101;
                              tetra[it5].ip[3]=i110;

                              tetra[it6].ip[0]=i011;     /*right bottom tetra*/
                              tetra[it6].ip[1]=i101;
                              tetra[it6].ip[2]=i111;
                              tetra[it6].ip[3]=i110;

                              tetra[it1].ifacet[0]=if1;  /*upper left tetra*/
                              tetra[it1].ifacet[1]=if3;
                              tetra[it1].ifacet[2]=if11;        
                              tetra[it1].ifacet[3]=if5;

                              tetra[it2].ifacet[0]=if5;  /*middle left tetra*/
                              tetra[it2].ifacet[1]=if9;
                              tetra[it2].ifacet[2]=if12;        
                              tetra[it2].ifacet[3]=if7;

                              tetra[it3].ifacet[0]=if7;  /*lower left tetra*/
                              tetra[it3].ifacet[1]=if10;
                              tetra[it3].ifacet[2]=if4;        
                              tetra[it3].ifacet[3]=if1+ntrap;

                              tetra[it4].ifacet[0]=if2;  /*upper right tetra*/
                              tetra[it4].ifacet[1]=if3;
                              tetra[it4].ifacet[2]=if9+1;        
                              tetra[it4].ifacet[3]=if6;

                              tetra[it5].ifacet[0]=if6;  /*middle right tetra*/