tetramod.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* TETRAMOD: $Revision: 1.1 $ ; $Date: 1997/10/22 18:22:32 $       */

#include "tetra.h"

/*********************** self documentation *****************************/
char *sdoc[] = {
"						 			",
" TETRAMOD - TETRAhedron MODel builder. In each layer, velocity gradient",
"       is constant or a 2-D grid; horizons could be a uniform grid and/or",
"       added by a 2-D grid specified.   				",
"									",
" tetramod [parameters] > tetrafile		 			",
" 									",
" Required parameters:							",
"                                                                       ",
" nxhz=		number of samples (2nd dimension) for horizons		",
" nyhz=		number of samples (1st dimension) for horizons		",
" hzfile=	output xhz,yhz,zhz,v0hz,v1hz for viewer3                ",
"                                                                       ",
" Optional parameters:							",
"                                                                       ",
" xmin=0       	x of the lower left point in the model	        	",
" ymin=0        y of the lower left point in the model                  ",
" xmax=2        x of the upper right point in the model                 ",
" ymax=2        y of the upper right point in the model                 ",
" zmax=(max z)  max z in the model                                      ",
" blt=1.0       bottom layer thickness                                  ",
" nhz=1		number of layers in the model (except the model base)	",
" ficth=-1      ficticious horizons (no velocity interpolation based on them)",
" The following four numbers define the four corners for 3-D model;	",
" z00=0,0.6,1.2,... 	z at (xmin,ymin) on each horizon		",
" z01=0,0.6,1.2,... 	z at (xmin,ymax) on each horizon		",
" z10=0,0.6,1.2,... 	z at (xmax,ymin) on each horizon                ",
" z11=0,0.6,1.2,... 	z at (xmax,ymax) on each horizon                ",
" v00=1.0,2.0,3.0,...	v at (xmin,ymin) on each horizon		",
" v01=1.0,2.0,3.0,...	v at (xmin,ymax) on each horizon		",
" v10=1.0,2.0,3.0,...	v at (xmax,xmin) on each horizon                ",
" v11=1.0,2.0,3.0,...	v at (xmax,xmax) on each horizon                ",
" dvdz00=0,0,0,...   	dvdz at (xmin,ymin) on each horizon		",
" dvdz01=0,0,0,...   	dvdz at (xmin,ymax) on each horizon		",
" dvdz10=0,0,0,...   	dvdz at (xmax,xmin) on each horizon             ",
" dvdz11=0,0,0,...   	dvdz at (xmax,xmax) on each horizon             ",
" x0file=   	x grid for horizon 0 					",
" y0file=   	y grid for horizon 0 					",
" z0file=	z grid for horizon 0 added to 4-z interpolation   	",
" v0file=	v grid for horizon 0 added to 4-v interpolation   	",
" dvdz0file=	dvdz grid for horizon 00 added to 4-dvdz interpolation	",
"									",
" x1file=	x grid for horizon 1					",
" y1file=	y grid for horizon 1					",
" z1file=	z grid for horizon 1 added to 4-z interpolation  	",
" v1file=	v grid for horizon 1 added to 4-v interpolation  	",
" dvdz1file=	dvdz grid for horizon 1 added to 4-dvdz interpolation	",
" ...		for horizon  #2, #3 #4, etc... 				",
" verbose=0	=1 print some useful iinformation			",
"									",
" Remarks:								",
" TETRAMOD defines its own grammar to describe a 3-D model (including a ",
" tetra file for ray tracer SUTETRARAY and a horizon file for VIEWER3   ",
"									",
" Disclaimer:                                                           ",
" This is a research code that will take considerable work to get into  ",
" the form of a a production level 3D migration code. The code is       ",
" offered as is, along with tetramod and sukdmig3d, to provide a starting",
" point for researchers who wish to write their own 3D migration codes.",
"                                                                       ",
NULL};

/*
 * Credits:
 *  	CWP: Zhaobo Meng, 1996
 *
 *  Reference:
 *  Zhaobo Meng and Norman Bleistein, Wavefront Construction (WF) Ray
 * "Tracing in Tetrahedral Models -- Application to 3-D traveltime and
 *  ray path computations, CWP report 251, 1997
 *
 */
/**************** end self doc *******************************************/

#define EPS 0.0001
#define SmallReal 0.1e-10
#define SWAP(a,b) {float temp=(a);(a)=(b);(b)=temp;}
#define np 9
#define maxiterations 5

/* structures used internally */
void write_point(
/************************************************************
write_point - write control points to a file
************************************************************
Function Prototype:
void write_point(struct POINT *point,int npoint);
************************************************************
Input:
int npoint number of points
Ouput:
struct POINT *point control points
************************************************************/
      struct POINT *point,
      int npoint);

void write_facet(
/************************************************************
Input:
int nfacet number of facets
Ouput:
struct FACET *facet facet information
************************************************************/
      struct FACET *facet,
      int nfacet);

void write_tetra(
/************************************************************
Input:
int ntetra number of tetra
Ouput:
struct TETRA *tetra tetra information
************************************************************/
      struct TETRA *tetra,
      int ntetra);

void tm_normal(
/************************************************************
Inputs:
float x0,y0,z0 first point
float x1,y1,z1 second point
float x2,y2,z2 third point
outputs:
float *n1,*n2,*n3 the normal
************************************************************/
      float x0,float x1,float x2,
      float y0,float y1,float y2,
      float z0,float z1,float z2,
      float *n1,float *n2,float *n3);

void tm_normalize(
/************************************************************
Inputs/outputs:
float *n1,*n2,*n3 normalize this vector
************************************************************/
float *n1,float *n2,float *n3);

void gradient(
/************************************************************
Inputs:
float x0[3],x1[3],x2[3],x3[3] 4 vertices
float s0,s1,s2,s3 sloth at the 4 vertices
Outputs:
float gs[3] the gradient of sloth
************************************************************/
      float x0[3],  /*point0*/
      float s0,     /*sloth at point0*/
      float x1[3],  /*point1*/
      float s1,     /*sloth at point1*/
      float x2[3],  /*point2*/
      float s2,     /*sloth at point2*/
      float x3[3],  /*point3*/
      float s3,     /*sloth at point3*/
      float gs[3]);

int gaussj(
/************************************************************
Inputs:
float a[3][3] 
Input/Output:
float b[3] the soluation of a x = b
************************************************************/
      float a[3][3],
      float b[3]);

float norm2(float xy[],int n);

void ax(
/************************************************************
Inputs:
float a[np][np]
float x[np] 
Output:
float y[3] y=ax
************************************************************/
float a[np][np],
float x[np],
float y[np]);

void tetracgr(
/************************************************************
Inputs:
float a[np][np]
float rhs[np] 
Input/Output:
float ct[np+1] the soluation of a x = rhs, x->ct
************************************************************/
float a[np][np],
float ct[np+1],
float rhs[np]);

int 
main(int argc, char **argv)
{
      int nhz; 	      /*number of horizons except the surface*/
      int nxhz;       /*number of samples in x */
      int nyhz;       /*number of samples in y (1st dimension)*/

      float dxhz;    
      float dyhz;   

      float ***xhz;   /*x of upper horizon: fastest x; slowest ihz*/
      float ***yhz;   /*y of upper horizon: fastest x; slowest ihz*/
      float ***zhz;   /*z of upper horizon: fastest x; slowest ihz*/
      float ***v0hz;  /*v below upper horizon: fastest x; slowest y*/
      float ***v0hztem;/*an extra copy of v0hz*/
      float ***v1hz;  /*v at the bottom of this region*/
      float ***dvdzhz;/*v'z below upper horizon: fastest x; slowest y*/
      float ***dvdzhztem;

      float ***n1hz;  /*normals on horizon vertexes*/
      float ***n2hz;  /*normals on horizon vertexes*/
      float ***n3hz;  /*normals on horizon vertexes*/

      /*The following quatities are 1 smaller in size in each 
      dimension than n1hz, n2hz and n3hz. This is because the
      following quantities are the normals to the triangles*/
      float **n1ll;   /*n1ll[nyhz-1][nxhz-1]*/
      float **n1ur;   /*n1ur[nyhz-1][nxhz-1]*/
      float **n2ll;   /*n2ll[nyhz-1][nxhz-1]*/
      float **n2ur;   /*n2ur[nyhz-1][nxhz-1]*/
      float **n3ll;   /*n3ll[nyhz-1][nxhz-1]*/
      float **n3ur;   /*n3ur[nyhz-1][nxhz-1]*/

      int npoint;          /*number of control points*/
      struct POINT *point; /*control points of size npoints*/

      int nfacet;          /*number of facets*/
      struct FACET *facet; 

      int ntetra;          /*number of tetra*/
      struct TETRA *tetra; /*tetra of size: nhz*(nxhz-1)*(nyhz-1)*6*/

      /**********************************************************
      The users need to define horizons this way: first gives the
      z values at the four corners for each horizon (if not given,
      0 value will be assigned), the horizon will be linearly 
      interpolated using the four values; then the user could define a
      z grid, then this grid will be added to the interpolated 
      horizon (of course if not given, nothing will be added to the
      interpolated horizon.
      Likewise, v, dvdz are all processed this way.
      ***********************************************************/ 
      float *z00; 	/*z at (xmin,ymin) on each horizon*/
      float *z01; 	/*z at (xmin,ymax) on each horizon*/
      float *z10; 	/*z at (xmax,ymin) on each horizon*/
      float *z11; 	/*z at (xmax,ymax) on each horizon*/

      float *v00; 	/*v at (xmin,ymin) on each horizon*/
      float *v01; 	/*v at (xmin,ymax) on each horizon*/
      float *v10; 	/*v at (xmax,ymin) on each horizon*/
      float *v11; 	/*v at (xmax,ymax) on each horizon*/

      float *dvdz00; 	/*dvdz at (xmin,ymin) on each horizon*/
      float *dvdz01;	/*dvdz at (xmin,ymax) on each horizon*/
      float *dvdz10;    /*dvdz at (xmax,ymin) on each horizon*/
      float *dvdz11;    /*dvdz at (xmax,ymax) on each horizon*/

      char *xfile;      /*grid file for x0*/
      char *yfile;	/*grid file for y0*/
      char *zfile;	/*grid file for z0 (added to z)*/
      char *vfile;	/*grid file for v0 (added to v)*/
      char *dvdzfile;   /*grid file for dvdz (added to dvdz)*/

      /*********************************************************
      These are the output. Will be used by viewer3, suwfrayt3d
      and trip, etc.
      *********************************************************/
      char *hzfile; 	  /*horizon file for viewer3 and suwfrayt3d*/
      FILE *hzfp;

      float ***zflat;
      float vtemp;

      int *ficth,ifict;
      int nficth;

      /*define the four corners of the model*/
      float xmin,ymin;    /*first reference point*/
      float xmax,ymax;    /*second reference point*/
      float zmax=-999999.9;
      float zmin=0.0;

      float a[np][np],b[np];

      int ixhz,iyhz;
      int ihz; 	          /*index for nhz*/
      int *iihz,iiihz;           /*index for ihz*/
      int verbose;        /*=1 then print some useful iinformation*/
      float fz0,fz1,fv0,fv1,fdvdz0,fdvdz1;
      float beta,alpha;
      float dis0=0.0;
      float dis1=0.0;
      float dis2=0.0;
      int i000,i001,i010,i100,i110,i101,i011,i111;
      int if1,if2,if3,if4,if5,if6,if7,if8,if9,if10,if11,if12;
      int it1,it2,it3,it4,it5,it6;
      int ntrap,nhztrap,nhztrap1;
      int nhzfacetx,nfacetx;
      int nhzfacety,nfacety;
      int iwh,ip;
      int nhzpoint,ihz0=0;
      float blt;         /*bottom layer thickness*/
      int ip0,ip1,ip2,ip3;
      int itetra;
      int jumphz,if1p,if2p;
      int ifs;

      /* hook up getpar */
      initargs(argc,argv);
      requestdoc(1);

      xfile=(char *)malloc(80);
      yfile=(char *)malloc(80);
      zfile=(char *)malloc(80);
      vfile=(char *)malloc(80);
      dvdzfile=(char *)malloc(80);

      /* get parameters */
      if (!getparint("nhz",&nhz))	
       	    nhz=1;
      if (!getparint("nxhz",&nxhz))   
	    err("Must specify nxhz");
      if (!getparint("nyhz",&nyhz))   
	    err("Must specify nyhz");