sutetraray.c

su 的源代码库

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

/* SUTETRARAY: $Revision: 1.3 $ ; $Date: 2006/11/07 23:11:28 $	*/

#include "su.h"
#include "segy.h"
#include "Tetra/tetra.h"

/*********************** self documentation *****************************/
char *sdoc[] = {
"									",
" SUTETRARAY - 3-D TETRAhedral WaveFront construction RAYtracing	",
"									",
"   sutetray < tetrafile [parameters] > ttfile	 			",
" 									",
" 									",
" Optional Parameters:							",	
" nxgd=				number of x samples (fastest dimension)	",
" nygd=				number of y samples (second dimension) 	",
" fxgd=xmin(from tetrafile)	first x sample			  	",
" fygd=ymin(from tetrafile)	first y sample				",
" shootupward=1			=1 shooting upward, -1 downward		",
" dxgd=(xmax-xmin)/(nxgd-1)	x sampling interval		 	",
" dygd=(ymax-ymin)/(nygd-1)	y sampling interval			",
" takeoff=30.0			max takeoff angle for shooting aperture	",
" ntakeoff=5			number of samples in takeoff angle	",
" nazimuth=15			number of samples in azimuth		",
" maxntris=1000			maximum number of triangles allowed	",
" maxnrays=600			maximum number of rays allowed		",
" maxnsegs=100			maximum number of ray segments allowed 	",
" tmax=10			maximum traveltime traced		",
" dtwf=0.05			wavefront step size for raytracing	",
" ntwf=(int)(tmax/dtwf)+1	maximum number of wavefronts	  	",
" ntwf=MAX(ntwf,ifwf2dump+nwf2dump)				     	",
" edgemax=4			maximum triangle length (in dxgd) not to",
"				be split				",
" fsx=fxgd+nxgd/2*dxgd		first source in x			",
" dsx=dxgd			x increment in source			",
" nsx=1				number of source in x			",
" fsy=fygd+nygd/2*dygd		first source in y			",
" dsy=dygd			y increment in source		   	",
" nsy=1				number of source in y		   	",
" fsz=zmin+(zmax-zmin)/2	first source in z			",
" dsy=(zmax-zmin)/2		y increment in source		   	",
" nsy=1				number of source in y			",
" nxt=nxgd*2/2+1	number of x-samples in ttable for sukdmig3d	",	
" nyt=nygd*2/2+1	number of y-samples in ttable for sukdmig3d	",
" irefseq=nhz-1,...     reflector sequence				",
" crfile=NULL	   file for cosines and ray paths r (for sukdmig3d)	",
" sttfile=NULL	  file for surface traveltimes (for visualization)	",
" verbose=0		=1 print some useful information		",
"									",
"	The following two files are output for viewer3 to display   	",
" rayfile=NULL		ray path file		    			",
" wffile=NULL		wavefront file (then must give ifwf2dump)	",
" ifwf2dump=20		the first wavefront to dump  			",
" nwf2dump=1		the number of wavefronts to dump		",
" jpfile=stderr		job print (default: stderr)			",
"									",
" 									",
" 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 DEBUG*/
#define OUT -99999		/* x and y value is the ray is out-of-domain */
#define EPS 0.00001		/* small number to avoid being divided by 0 */
#define TraveltimeInfinity 9.0  /* traveltime infinity */
#define DistanceInfinity 9000	/* distance infinity */
#define SmallReal 0.1e-5 	/* must be smaller than any slowness */
#define SmallSigma 0.0001	/* smallest sigma value for comparison */
#define SigmaInfinity 1.0e+10
#define MaxOnce 1000
#define Margin 0.1              /* some margin for which_tetra */

segy tro;

struct WF
{
      float x[3];   /* intersection */
      float p[3];   /* slowness vector at intersection */	
};

struct TRI
{
      int i1;	        /* index for first vertex */
      int i2;           /* index for second vertex */
      int i3;           /* index for third vertex */
      struct TRI *next;	/* pointer to next tri */
};

enum ICODE {NORMAL,SPOILED,HIT_SURFACE};

struct RAY
{
      enum ICODE icode;
      int itetra;     /* tetra in which this ray now travels */
      int lastfacet;  /* last facet that this ray penetrates */
      int ionce;      /* number of adjustments applied when cornered */
      int iref;       /* index for reflections */
      int shootupward;/* 1 upward, -1 downward */
      int nseg;       /* number of segments */
      float t;        /* current traveltime within dt */
      float x[3];     /* cuttent ray position */
      float ttotal;   /* total traveltime from source */
      float r;        /* total distance traveled from source */
      float v;        /* velocity at current point */
      float **xseg;   /* segment positions */
      float *tseg;
      float p[3];
};

enum RESULT {HIT_FACET_FIRST,
      HIT_WF_FIRST,
      HIT_SURFACE_FIRST,
      OUT_OF_MODEL};

enum CHECK_TRIS {PLEASE,PLEASE_DO_NOT};

enum DONE {ALREADY,NOT_YET};

enum IN_TRI{WF_YES,WF_NO};

enum ENOUGH_NTWF{WF_TRUE,WF_FALSE};
/* Prototypes of functions used internally */
void crossprod(float u1, float u2, float u3, float v1, float v2,
		 float v3, float *w1, float *w2, float *w3);
float tetrasolver(float ct[10], float p[3], float gs[3],
			float xst[3],float sg0); 
int in_which_tetra(float xmin, float ymin, float xmax,
			float ymax, float zmax, float x[3],
			 struct POINT *point, struct TETRA *tetra,
			 int ntetra);
float tri_intp(float x, float y, float x0, float y0, float t0,
		float x1, float y1, float t1, float x2, float y2, float t2); 
void tri_multi_intp( struct RAY *ray0, struct RAY *ray1, struct RAY *ray2,
		float fxgd, float fygd, float dxgd, float dygd,
		int nxgd, int nygd, float **ttable, float **ctable,
		 float **rtable);
enum IN_TRI in_tri(float x0, float y0, float x1, float y1, float x2,
		float y2, float x, float y, float areat);
float quadsolver(float a,float b,float c);
float cubicsolver(float d,float a,float b,float c);

enum RESULT trace1tetra( int maxnsegs, float xmin, float ymin, float xmax,
      float ymax, float zmax, struct RAY *ray, int npoint,
      struct POINT *point, int nfacet, struct FACET *facet,
      int ntetra, struct TETRA *tetra, int *irefseq,
      int nrefseq, FILE *jpfp, float dt);  

int    
oneshot(int shootupward, int maxnsegs, int npoint, struct POINT *point,
        int nfacet, struct FACET *facet, int ntetra, struct TETRA *tetra, 
        int *irefseq, int nrefseq, FILE *jpfp, float takeoff, float dtwf,
        int ntwf, int maxntris, int maxnrays, float edgemax, 
        int nxt, int nyt, int ixtf, int iytf, float source[3], int ntakeoff,
        int nazimuth, float fxgd, float dxgd, int nxgd, float fygd, 
        float dygd, int nygd, float xmin, float ymin, float xmax, float ymax,
        float zmax, FILE *rayfp, FILE *wffp, FILE *crfp, FILE *sttfp,
        int ifwf2dump, int nwf2dump);
/*******************************************************************
Inputs:
int shootupward flag: shootupward if =1 else downward
int maxnsegs    max number of ray segments
int npoint      number of point
struct *point   control points
int nfacet      number of facet
struct *facet   facet information
int ntetra      number of tetra
struct *tetra   tetra information
int *irefseq    reflector sequence
int nrefseq     number of reflectors
FILE *jpfp	file pointer to the job print
float takeoff   maximum takeoff angle
float dtwf        time step for wavefront updating
int ntwf          total number of time steps allowed
int maxntris    maximum number of triangles allowed
int maxnrays    maximum number of rays allowed
float edgemax   maximum edge length in dxgd
int nxt         number of x-samples in ttable for sukdmig3d
int nyt         number of y-samples in ttable for sukdmig3d
int nxtl        number of samples in ttable on left of source
int nytf        number of samples in ttable in front of source
float source[3] source position
int ntakeoff    number of samples in azimuth
int nazimuth    number of take-off angles
float fxgd      first sample in x
float dxgd      sample interval in x
int nxgd        number of sampels in x
float fygd      first sample in y
float dygd      sample interval in y
int nygd        number of sampels in y
FILE *rayfp     file pointer to rayfile
FILE *wffp      file pointer to wavefront file
int ifwf2dump   the first wavefront to dump
int nwf2dump    the number of wavefronts to dump

Output:         none
******************************************************************/

float f(float sigma, float s, float p, float a);
float zhorz( float xx, float yy, float n10, float n20, float n30,
      float x00, float y00, float z00);
void WF_normalize(float *n1,float *n2,float *n3);
void read_point( struct POINT *point, int npoint);
void read_facet( struct FACET *facet, int nfacet);
void read_tetra( struct TETRA *tetra, int ntetra);
float s_intp( struct POINT *point, struct TETRA *tetra, float x[3]);
float area( float x0, float y0, float x1, float y1, float x2, float y2);
void cross_slim_tetra( struct RAY *ray, struct POINT *point,
		 struct FACET *facet, struct TETRA *tetra, int *irefseq,
		int nrefseq, float sincid, FILE *jpfp); 
float facet_fit( struct POINT *point, float x[3], struct FACET *facet);
enum RESULT RTlaw( struct RAY *ray, struct POINT *point, struct FACET *facet,
      struct TETRA *tetra, int *irefseq, int nrefseq, int ifhit, float sincid,
      int itetranew, FILE *jpfp);
void s_normalize( float p[3], float s);

int 
main(int argc, char **argv)
{
      int shootupward;      /* flag: 1 upward, -1 downward */
      int npoint;           /* number of control points */
      struct POINT *point=NULL;/* control points */
      int nfacet;           /* number of facets */
      struct FACET *facet=NULL;/* facet information */
      int ntetra;           /* number of tetra */
      struct TETRA *tetra=NULL;/* tetra information */

      int nygd;    	    /* number of samples in y direction */
      int nxgd; 	    /* number of samples in x direction */

      /**********************************************************
      These four parameters are used for defining the ttable 
      size for sukdmig3d
      **********************************************************/
      int nxt;	      /* number of x-samples in ttable for sukdmig3d */
      int nyt;	      /* number of y-samples in ttable for sukdmig3d */
      int ixtf;       /* number of samples in ttable on left of source */
      int iytf;       /* number of samples in ttable in front of source */
      int ntakeoff;   /* number of samples in take-off */
      int nazimuth;   /* number of samples in azimuth */

      float fygd;     /* first sample in y */
      float fxgd;     /* first sample in x */

      float fsx;      /* first source in x */
      float fsy;      /* first source in y */
      float fsz;      /* first source in z */

      float dsx;      /* x increment of sources */
      float dsy;      /* y increment of sources */
      float dsz;      /* z increment of sources */

      int nsx;	      /* number of sources in x direction */
      int nsy;	      /* number of sources in y direction */ 
      int nsz;        /* number of sources in z direction */

      int *irefseq;   /* reflector sequence */

      float dygd;     /* y sample interval */
      float dxgd;     /* x sample interval */

      float dtwf;	      /* time step in raytracing */
      int ntwf;	      /* number of time steps in raytracing */
      float tmax;     /* maximum time traced */

      char *rayfile="";/* file for ray information */
      char *wffile="";
      char *crfile="";
      char *sttfile="";
      char *jpfile="";/* file for job print */

      FILE *rayfp=NULL;
      FILE *rayfp0=NULL;
      FILE *crfp=NULL;
      FILE *sttfp=NULL;
      FILE *wffp=NULL;
      FILE *jpfp=stderr;

      int ifwf2dump;
      int nwf2dump;

      int nrefseq;
      int ns;         /* number of sources */

      int maxntris;   /* maximum number of triangles allowed */
      int maxnrays;   /* maximum number of rays allowed */
      float edgemax;  /* maximum edge length in dxgd */
	
      int verbose;    /* if =1 print some useful information */

      float source[3];/* source position */

      int isx;	      /* source index for nsx */
      int isy;	      /* source index for nsy */
      int isz;        /* source index for nsz */

      int ixsf;       /* first source in gd */
      int ixsr;       /* source spacing and gd spacing ratio */
      int iysf;       /* first source in gd */
      int iysr;       /* source spacing and gd spacing ratio */

      float xmin,xmax,ymin,ymax,zmin,zmax;

      float takeoff;  /* max takeoff angle */

      int iref,maxnsegs;

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

      /* get parameters */
      if (!getparint("verbose",&verbose)) verbose=0;

      if (!getparint("shootupward",&shootupward)) shootupward=1;

      if (getparstring("jpfile",&jpfile))
	    if ((jpfp=fopen(jpfile,"w"))==NULL)
		 jpfp=stderr;
		 
      if (!getparfloat("takeoff",&takeoff)) takeoff=30.0;
      takeoff=takeoff/180.0*PI;

      nrefseq=countparval("irefseq");
      irefseq=alloc1int(nrefseq);
      getparint("irefseq",irefseq);

      if (verbose)
            for (iref=0;iref<nrefseq;iref++)
                  fprintf(jpfp,"irefseq[%d]=%d\n",iref,irefseq[iref]);

      fscanf(stdin,"%f=xmin\n",&xmin);
      fscanf(stdin,"%f=xmax\n",&xmax);
      fscanf(stdin,"%f=ymin\n",&ymin);
      fscanf(stdin,"%f=ymax\n",&ymax);
      fscanf(stdin,"%f=zmin\n",&zmin);
      fscanf(stdin,"%f=zmax\n",&zmax);

      if (verbose)
	    fprintf(jpfp,"xmin=%f\nxmax=%f\nymin=%f\nymax=%f\nzmin=%f\nzmax=%f\n",
		  xmin,xmax,ymin,ymax,zmin,zmax);

      if (!getparint("nxgd",&nxgd)) nxgd=21;
      if (nxgd<2) err("nxgd too small");
      nxgd=(nxgd/2)*2+1;

      if (!getparfloat("fxgd",&fxgd)) fxgd=xmin;
      if (fxgd<xmin) err("fxgd too small");
      if (!getparfloat("dxgd",&dxgd)) dxgd=(xmax-xmin)/(nxgd-1);
      if (fxgd+(nxgd-1)*dxgd>xmax*1.01) {
            fprintf(stderr,"fxgd=%f nxgd=%d dxgd=%f xmax=%f\n",
            fxgd,nxgd,dxgd,xmax);
            err("fxgd+(nxgd-1)*dxgd too large");
      }

      if (!getparint("nygd",&nygd)) nygd=21;