tetramod.c

su 的源代码库

      if (!getparstring("hzfile",&hzfile)) 	
       	    err("Must specify hzfile");
 
      if (!getparint("verbose",&verbose)) 	
       	    verbose=0;

      nficth=countparval("ficth");
      ficth=alloc1int(nficth);
      getparint("ficth",ficth);
      if (nficth==1 && ficth[0]==-1) nficth=0;

      if (verbose) {
	    fprintf(stderr,"nficth=%d\n",nficth);
            for (ifict=0;ifict<nficth;ifict++)
                  fprintf(stderr,"ficth[%d]=%d",ifict,ficth[ifict]);
      }

      if (!getparfloat("blt",&blt)) blt=1.0;

      if (verbose)
	    warn("hzfile=%s\n",hzfile);

      if (!getparfloat("xmin",&xmin)) xmin=0.0;
      if (!getparfloat("xmax",&xmax)) xmax=2.0;
      if (!getparfloat("ymin",&ymin)) ymin=0.0;
      if (!getparfloat("ymax",&ymax)) ymax=2.0;
      if (!getparfloat("zmax",&zmax)) zmax=0.0;

      if ((hzfp=fopen(hzfile,"w"))==NULL)
       	    err("Can not open hzfile");

      if (verbose) {
       	    warn("nhz=%d",nhz);
            warn("xmin=%g\nxmax=%g\nymin=%g\nymax=%g\n",xmin,xmax,ymin,ymax);
      }
	
      z00   =alloc1float(nhz);
      z01   =alloc1float(nhz);
      z10   =alloc1float(nhz);
      z11   =alloc1float(nhz);
      v00   =alloc1float(nhz);	
      v01   =alloc1float(nhz);
      v10   =alloc1float(nhz);
      v11   =alloc1float(nhz);
      dvdz00=alloc1float(nhz);
      dvdz01=alloc1float(nhz);
      dvdz10=alloc1float(nhz);
      dvdz11=alloc1float(nhz);

      iihz=ealloc1int(nhz);

      if (!getparfloat("z00",z00))
	    for(ihz=0;ihz<nhz;ihz++) 
	       	  z00[ihz] = ihz*0.6;
      if (!getparfloat("z01",z01))
	    for(ihz=0;ihz<nhz;ihz++) 
		  z01[ihz] = ihz*0.6;
      if (!getparfloat("z10",z10))
	    for(ihz=0;ihz<nhz;ihz++) 
		  z10[ihz] = ihz*0.6;
      if (!getparfloat("z11",z11))
	    for(ihz=0;ihz<nhz;ihz++) 
	       	  z11[ihz] = ihz*0.6;

      if (verbose) {	
       	    for (ihz=0;ihz<nhz;ihz++) 
	       	  warn("layer %d z: %f %f %f %f",ihz,
			z00[ihz],z01[ihz],z10[ihz],z11[ihz]);
      }

      /************************************************
      The default velocity is 2.0,3.0,...,in km/s
      ************************************************/
      if (!getparfloat("v00",v00))
	    for(ihz=0;ihz<nhz;ihz++) v00[ihz] = ihz+1;
      if (!getparfloat("v01",v01))
	    for(ihz=0;ihz<nhz;ihz++) v01[ihz] = ihz+1;
      if (!getparfloat("v10",v10))
	    for(ihz=0;ihz<nhz;ihz++) v10[ihz] = ihz+1;
      if (!getparfloat("v11",v11))
	    for(ihz=0;ihz<nhz;ihz++) v11[ihz] = ihz+1;

      if (verbose) 
	    for (ihz=0;ihz<nhz;ihz++) 
		  warn("layer %d v: %f %f %f %f",ihz,
		       v00[ihz],v01[ihz],v10[ihz],v11[ihz]);

      if (!getparfloat("dvdz00",dvdz00))
	    for(ihz=0;ihz<nhz;ihz++) dvdz00[ihz] = 0.0;
      if (!getparfloat("dvdz01",dvdz01))
	    for(ihz=0;ihz<nhz;ihz++) dvdz01[ihz] = 0.0;
      if (!getparfloat("dvdz10",dvdz10))
	    for(ihz=0;ihz<nhz;ihz++) dvdz10[ihz] = 0.0;
      if (!getparfloat("dvdz11",dvdz11))
	    for(ihz=0;ihz<nhz;ihz++) dvdz11[ihz] = 0.0;

      if (verbose) 
	    for (ihz=0;ihz<nhz;ihz++) 
		  warn("layer %d dvdz: %f %f %f %f",ihz,
		       dvdz00[ihz],dvdz01[ihz],dvdz10[ihz],dvdz11[ihz]);

      /*************************************************
      Get the size of the horizon (in triangles), which is
      not necessary to be equal to nx and ny, and these
      grids are not uniform 
      *************************************************/
      warn("nxhz=%d",nxhz);
      warn("nyhz=%d",nyhz);
      if (nxhz<=1) err("nxhz must >1");
      if (nyhz<=2) err("nyhz must >2");
      if (nyhz/2*2==nyhz) err("nyhz must be an odd number\n");

      n1hz=ealloc3float(nxhz,nyhz,nhz);
      n2hz=ealloc3float(nxhz,nyhz,nhz);
      n3hz=ealloc3float(nxhz,nyhz,nhz);

      n1ll=ealloc2float(nxhz-1,nyhz-1);
      n1ur=ealloc2float(nxhz-1,nyhz-1);
      n2ll=ealloc2float(nxhz-1,nyhz-1);
      n2ur=ealloc2float(nxhz-1,nyhz-1);
      n3ll=ealloc2float(nxhz-1,nyhz-1);
      n3ur=ealloc2float(nxhz-1,nyhz-1);

      xhz=ealloc3float(nxhz,nyhz,nhz);
      yhz=ealloc3float(nxhz,nyhz,nhz);
      zhz=ealloc3float(nxhz,nyhz,nhz+1);
      v0hz=ealloc3float(nxhz,nyhz,nhz);
      v0hztem=ealloc3float(nxhz,nyhz,nhz);
      v1hz=ealloc3float(nxhz,nyhz,nhz);
      dvdzhz=ealloc3float(nxhz,nyhz,nhz);
      dvdzhztem=ealloc3float(nxhz,nyhz,nhz);

      zflat=ealloc3float(nxhz,nyhz,nhz);

      /*these dxhz and dyhz will only be used if there is no xhz and yhz 
      provided*/

      dxhz=(xmax-xmin)/(nxhz-1);
      dyhz=(ymax-ymin)/(nyhz-1);

      for (ihz=0;ihz<nhz;ihz++) {

	    sprintf(xfile,"x%dfile",ihz);
	    sprintf(yfile,"y%dfile",ihz);
	    sprintf(zfile,"z%dfile",ihz);
	    sprintf(vfile,"v%dfile",ihz);
	    sprintf(dvdzfile,"dvdz%dfile",ihz);

	    if (verbose==1) 
		  warn("%s %s %s %s %s",xfile,
		       	yfile,zfile,vfile,dvdzfile);
		
	    if (getparstring(xfile,&xfile)) {
		  if (fread(xhz[ihz][0],sizeof(float),nxhz*nyhz,
		       	fopen(xfile,"r"))!=nxhz*nyhz)
		       	err("Can not read in %s\n",xfile);
		  warn("read in %s",xfile);
		  for (iyhz=0;iyhz<nyhz;iyhz++) {
		        if (fabs(xhz[ihz][iyhz][0]-xmin)>EPS) 
		              err("xhz[%d][%d][0] != xmin, not allowed",ihz,iyhz);
		        else
			      xhz[ihz][iyhz][0]=xmin;
			if (fabs(xhz[ihz][iyhz][nxhz-1]-xmax)>EPS)
			      err("xhz[%d][%d][%d]=%e != xmax=%e,not allowed",
			      ihz,iyhz,nxhz-1,xhz[ihz][iyhz][nxhz-1],xmax);
			else
			      xhz[ihz][iyhz][nxhz-1]=xmax;
		  }

	    } else { 
		  warn("Using uniform Cartesian grids for x-horizon%d",ihz);
		  for (ixhz=0;ixhz<nxhz;ixhz++)
		       	for (iyhz=0;iyhz<nyhz;iyhz++)
		       		xhz[ihz][iyhz][ixhz]=xmin+dxhz*ixhz; 
	    }

	    if (getparstring(yfile,&yfile)) {
	       	  if (fread(yhz[ihz][0],sizeof(float),nxhz*nyhz,
		       	fopen(yfile,"r"))!=nyhz*nxhz)
		       	      err("Can not read in %s\n",yfile);
		  warn("read in %s",yfile);
		  for (ixhz=0;ixhz<nxhz;ixhz++) {
		        if (fabs(yhz[ihz][0][ixhz]-ymin)>EPS) 
		      	      err("yhz[%d][0][%d]=%e != ymin=%e,not allowed",ihz,ixhz,
			      yhz[ihz][0][ixhz],ymin);
			else
			      yhz[ihz][0][ixhz]=ymin;
			if (fabs(yhz[ihz][nyhz-1][ixhz]-ymax)>EPS)
			      err("yhz[%d][%d][%d]=%e!=ymax=%e,not allowed",
			           ihz,nyhz-1,ixhz,yhz[ihz][nyhz-1][ixhz],ymax);
			else
			      yhz[ihz][nyhz-1][ixhz]=ymax;
		  }
	    } else {
		  warn("Using uniform Cartesian grids for y-horizon%d",ihz);
		  for (ixhz=0;ixhz<nxhz;ixhz++)
		        for (iyhz=0;iyhz<nyhz;iyhz++)
		      	      yhz[ihz][iyhz][ixhz]=ymin+iyhz*dyhz;
	    }

	    if (getparstring(zfile,&zfile)) {
		  if (fread(zhz[ihz][0],sizeof(float),nyhz*nxhz,
		        fopen(zfile,"r"))!=nyhz*nxhz)
		        err("Can not read in %s\n",zfile);
		  warn("read in %s",zfile);
	    } else {
	          warn("Using uniform cartesian grids for z-horizon%d",ihz);
	          for (ixhz=0;ixhz<nxhz;ixhz++)
	       	        for (iyhz=0;iyhz<nyhz;iyhz++)
			      zhz[ihz][iyhz][ixhz]=0;
	    }
	
	    if (getparstring(vfile,&vfile)) {       
	          if (fread(v0hz[ihz][0],sizeof(float),
	       	        nxhz*nyhz,fopen(vfile,"r"))!=
		       	      nxhz*nyhz)       
			      err("Can not read in %s\n",vfile);  
		  warn("read in %s",vfile);
	    } else {
		  warn("Using default for v-horizon %d",ihz);
	          for (ixhz=0;ixhz<nxhz;ixhz++)
		        for (iyhz=0;iyhz<nyhz;iyhz++)
		       	      v0hz[ihz][iyhz][ixhz]=0;
	    }

	    if (getparstring(dvdzfile,&dvdzfile)) {
	          if (fread(dvdzhz[ihz][0],sizeof(float),nyhz*nxhz,
			fopen(dvdzfile,"r"))!=nyhz*nxhz)       
		       	      err("Can not read in %s\n",dvdzfile);  
		  warn("read in %s",dvdzfile);
	    } else {
		  warn("Using default for dvdz-horizon %d",ihz);
		  for (ixhz=0;ixhz<nxhz;ixhz++)
			for (iyhz=0;iyhz<nyhz;iyhz++)
			      dvdzhz[ihz][iyhz][ixhz]=0;
	    }

	    /*************************************************************
	    zhz is added to the plane determined by the four corners
	    *************************************************************/
	    for (ixhz=0;ixhz<nxhz;ixhz++){
	          for (iyhz=0;iyhz<nyhz;iyhz++) {
		        alpha=(xhz[ihz][iyhz][ixhz]-xmin)/(xmax-xmin);
		        fz0=alpha*z01[ihz]+(1-alpha)*z00[ihz];
		        fz1=alpha*z11[ihz]+(1-alpha)*z10[ihz];
		        beta =(yhz[ihz][iyhz][ixhz]-ymin)/(ymax-ymin);
		        zflat[ihz][iyhz][ixhz]=beta*fz1+(1-beta)*fz0;
			zhz[ihz][iyhz][ixhz]+=zflat[ihz][iyhz][ixhz];
		  }
	    }
      }

      for (iyhz=0;iyhz<nyhz;iyhz++)
            for (ixhz=0;ixhz<nxhz;ixhz++)
	          zmax=MAX(zmax,zhz[nhz-1][iyhz][ixhz]);

      zmax+=blt;
      for (iyhz=0;iyhz<nyhz;iyhz++)
            for (ixhz=0;ixhz<nxhz;ixhz++)
	          zhz[nhz][iyhz][ixhz]=zmax;

      for (ihz=0;ihz<nhz;ihz++) {
  	    if (nficth!=0 && ihz!=0) {
	          for (iiihz=ihz;iiihz>=0;iiihz--) {
		        for (ifict=0;ifict<nficth;ifict++) 
                              if (iiihz==ficth[ifict]) break;
                        if (ifict==nficth) break;
                  }
                  iihz[ihz]=iiihz;
                  if (iiihz==nhz) iihz[ihz]=ihz;
            } else
                  iihz[ihz]=ihz;

            /* iihz is the previous non-fictitious horizon */
            fprintf(stderr,"iihz[%d]=%d\n",ihz,iihz[ihz]);
      }

      for (ihz=0;ihz<nhz;ihz++)
	    for (ixhz=0;ixhz<nxhz;ixhz++)
	          for (iyhz=0;iyhz<nyhz;iyhz++) {
                        v0hztem[ihz][iyhz][ixhz]=v0hz[ihz][iyhz][ixhz];
			dvdzhztem[ihz][iyhz][ixhz]=dvdzhz[ihz][iyhz][ixhz];
	          }
      for (ihz=0;ihz<nhz;ihz++){
	    /*************************************************************
	    v0hz and dvdzhz are added to the velocity plane determined by 
	    velocity at the four corners
	    *************************************************************/
	    for (ixhz=0;ixhz<nxhz;ixhz++){
	          for (iyhz=0;iyhz<nyhz;iyhz++){

                        iiihz=iihz[ihz];
                        if (ihz<nhz-1) {
			      if (zhz[ihz][iyhz][ixhz]>zhz[ihz+1][iyhz][ixhz]) { 
			            SWAP(zhz[ihz][iyhz][ixhz],zhz[ihz+1][iyhz][ixhz]);
                                    iiihz=iihz[ihz+1];
			      }
                        }

	     	      	alpha=(xhz[ihz][iyhz][ixhz]-xmin)/(xmax-xmin);
	       	        fv0=alpha*v01[iiihz]+(1-alpha)*v00[iiihz];
	       	        fv1=alpha*v11[iiihz]+(1-alpha)*v10[iiihz];
	       	        fdvdz0=alpha*dvdz01[iiihz]+(1-alpha)*dvdz00[iiihz];
	       	        fdvdz1=alpha*dvdz11[iiihz]+(1-alpha)*dvdz10[iiihz];
	       	        beta=(yhz[ihz][iyhz][ixhz]-ymin)/(ymax-ymin);

			vtemp=v0hztem[iiihz][iyhz][ixhz]+
	       	              beta*fv1+(1-beta)*fv0;
	       	        dvdzhz[ihz][iyhz][ixhz]=dvdzhztem[iiihz][iyhz][ixhz]+
	       	              beta*fdvdz1+(1-beta)*fdvdz0;

			v0hz[ihz][iyhz][ixhz]=vtemp+
			      dvdzhz[ihz][iyhz][ixhz]*(zhz[ihz][iyhz][ixhz]-
			      zflat[iiihz][iyhz][ixhz]);

                        v1hz[ihz][iyhz][ixhz]=vtemp+
                              dvdzhz[ihz][iyhz][ixhz]*(
                              zhz[ihz+1][iyhz][ixhz]-
                              zflat[iiihz][iyhz][ixhz]);

                        #ifdef DEBUG
                        if (ixhz==0 && iyhz==0) {
                              fprintf(stderr,"iiihz=%d v0hz=%f v1hz=%f \n",
                              iiihz,v0hz[ihz][iyhz][ixhz],v1hz[ihz][iyhz][ixhz]);
                              fprintf(stderr,"vtemp=%f\n",vtemp);
                              fprintf(stderr,"zflat=%f zhz=%f\n",
                              zflat[iiihz][iyhz][ixhz],zhz[ihz][iyhz][ixhz]);
                        }
                        #endif
		  }
	    }

	    /**********************************************************
	    Output the normals for suwfrayt3d
	    All triangles look like: