rayt3d.c

seismic software,very useful

#include "ttc2.h"
#include "usgrid.h"
#include "grid.h"
#include "su.h"
#include "header.h"
#include <assert.h>
#include <pthread.h>

#ifndef WORDS_BIGENDIAN
#include "ttc_byteswap.h"
#endif


char *sdoc[] = {
"rayt3d - traveltime tables calulated by paraxial ray tracing ",
" ",
"rayt3d vfile= tfile= [optional parameters] ",
" ",
"Required Parameters:",
"vfile=stdin            file containing velocitiy v(nz,nx,ny)    	",
"tfile=stdout           file containing traveltime tables  ",
"			t(nzo,nxo,nyo,nxs,nys)	    ",
" ",
"Optional Parameters:",
"dt=0.008  		time sampling interval in ray tracing	",
"nt=401  		number of time sampling in ray tracing	",
"                       (one-way travel time calculation)	",
" ",
"the following nine parameters are from velocity grid header ONLY",
"fz=from-vfile	        first depth sample in velocity 	",
"nz=from-vfile       	Number of depth samples in velocity 	",
"dz=from-vfile       	depth interval in velocity		",
"fx=from-vfile       	first inline coordinate in velocity 	", 
"nx=from-vfile		Number of inline coordinates in velocity	",
"dx=from-vfile		inline interval in velocity			",
"fy=from-vfile       	first crossline coordinate in velocity 	" ,
"ny=from-vfile       	Number of crossline coordinates in velocity	",
"dy=from-vfile        	crossline interval in velocity			",
" ",
"fxo=fx                 first x coordinate in traveltime table",
"fyo=fy                 first y coordinate in traveltime table",
"fzo=fz                 first z coordinate in traveltime table",
"dxo=dx                 x interval in traveltime table",
"dyo=dy                 y interval in traveltime table",
"dzo=dz                 z interval in traveltime output",
"nxo=nx                 number of x samples in traveltime table",
"nyo=ny                 number of y samples in traveltime table",
"nzo=nz                 number of z samples in traveltime table",
"ntline=0		number of target lines (inlinel, crossline or	",
"			slant line). Traveltime outputs on each target 	",
"			line instead of whole volume			",
"when ntline>0, specify 			",
"	v0=5000		reference velocity at the surface		",
"	dvz=0		reference velocity vertical gradient		",
"and for each target line specify					",
"	xbeg		x coordinate at which the line begins		",
"	xend		x coordinate at which the line ends		",
"	ybeg		y coordinate at which the line begins		",
"	yend		y coordinate at which the line ends		",
"	nout		number of output points in the line		",
"	tofile		filename containing traveltime output  		",
"                       tfile will be ignored				",
"                       the above 6 parameters must be given ntline times ",
" ",
"nxs=1                  number of x samples for source locations	",
"nys=1                  number of y samples for source locations	",
"fxs=fx                 x coordinate of first source ",
"fys=fy                 y coordinate of first source ",
"dxs=2*dxo              x coordinate increment of sources ",
"dys=2*dyo              y coordinate increment of sources ",
"aperx=0.5*nx*dx        ray tracing aperature in x-direction ",
"apery=0.5*ny*dy        ray tracing aperature in y-direction ",
" ",
"fa=0            	first take-off polar angle ",
"da=1            	increment of take-off polar angle ",
"na=61            	number of polar angles",
"amin=0            	minimum emergence polar angle  ",
"                       measured from vertical --- moving upwards	",
"amax=90            	maximum emergence polar angle	",
"azhmin=0            	minimum take-off azimuth angle ",
"                       measured from y axis (clockwise)	",
"azhmax=360            	maximum take-off azimuth angle",
" ",
"fac=0.01            	factor to determine radius for extrapolation	",
"ek=1              	=1 to implement eikonal equation in shadow zones",
"ms=1			print verbal information at every ms finished sources",
"ncpu=1			number of cpu to use in the computation ",
"jpfile=stderr          name of job print file; default to standard error",
"restart=n		job is restarted (=y yes; =n no)		",
"isres=                 restart source index (sequential) ",
"                       (default will be determined from the travel ",
"                        file size) ",
"diskvxx=       disk file for velocity derivative	 		",
"diskvxy=       disk file for velocity derivative			",
"diskvxz=       disk file for velocity derivative			",
"diskvyy=       disk file for velocity derivative			",
"diskvyz=       disk file for velocity derivative			",
"diskvzz=       disk file for velocity derivative			",
"               (if the above 6 parameters are not given, ",
"               the memory is going to be used)				",
"idiskc=0       the above 6 disk files already computed (0=no 1=yes)", 
"compress=0     travel time compress ratio (0=no 5=ratio of 5) ",
" ",
"Note:	",
"1. Each traveltime table is calculated by paraxial ray tracing; then 	",
"   traveltimes in shadow zones are computed by solving eikonal equation.",
"2. A smoothed velocity is prefered. All sampling information of the	",
"   velocity is included in the file grid header.			",
"3. Traveltime table and source locations must be within velocity model. ",
"4. Ray tracing apertures can be choosen as sum of migration apertures ",
"   plus half of maximum offset. ",
"5. When ntline=0, the program outputs the whole traveltime tables.",
"   Otherwise, the outputs are given at the single target lines.	",
"6. When a single target line is slant, the specified parameters for  	",
"   output must be consistent with those used in migration. 	",	
"7. memory requirement for this program is about,		",
"   (7*nz*(nx*ny+ncpu*nxyt)+(1+2*cpu)*nxo*nyo*nzo+28*128*1001)*4 bytes	",
"   where 1001 is the maximum number of nt, 128 is maximum number of	",
"   rays to compute at vector unit, 28 is the number of auxiliary arrays",
"   used in ray tracing, nxyt=(2*aperx/dx)*(2*apery/dy)			",
" ",
" ",
NULL};

/* This structure is used to pass arguments to compress_ttc through the
pthreads call. */
typedef struct
{
  float *data;
  int *num_points;
  float ratio;
  float null_value;
  int thread_number;
  unsigned char *work_space;
  long long  work_bytes;
  unsigned char *buffer;
} compress_ttc_call_struct;

/* This function is used to call compress_ttc from a thread.  */
void compress_ttc_pt(void *call_struct)
{
 compress_ttc_call_struct *tmp = (compress_ttc_call_struct*)call_struct;

 compress_ttc(tmp->buffer, tmp->data, tmp->num_points, tmp->ratio, tmp->null_value, tmp->thread_number, tmp->work_space, tmp->work_bytes);
}

/* multiple-cpu ray tracing subroutine */
void rt3dp_(int *np,int *nzyx,int *nzyxo,int *nz,int *ny,int *nx,
  int *nzo,int *nyo,int *nxo,int *ek,int *nt, int *na,
  float *fa,float *da,float *amin,float *amax,float *azhmin,float *azhmax,
  float *dt,float *tmax,
  float *fxo,float *fyo,float *fzo,float *dxo,float *dyo,float *dzo,
  float *fac,float *fx,float *fy,float *fz,
  float *dx,float *dy,float *dz,float *aperx,float *apery,
  float *v,float *vxx,float *vxy,float *vxz,float *vyy,float *vyz,
  float *vzz,float *ov2,
  float *vt,float *vxxt,float *vxyt,float *vxzt,float *vyyt,
  float *vyzt,float *vzzt,
  float *tt1,float *tt2,float *t,float *s,
  float *xsp,float *ysp,float *azhnp,float *azhxp,float *fxtp,float *fytp,
  int *nxtp,int *nytp,int *nzyxt,
  float *xp,float *yp,float *zp,float *pxp,float *pyp,float *pzp,
  float *e1xp,float *e1yp,float *e1zp,float *e2xp,float *e2yp,float *e2zp,
  float *q111p,float *q112p,float *q121p,float *q122p,float *p211p,
  float *p212p,float *p221p,float *p222p,float *q211p,float *q212p,
  float *q221p,float *q222p,float *vp,float *dvdxp,float *dvdyp,
  float *dvdzp,int *nrsp,float *a0p,float *azh0p,int *n1,int *n2,float *p2p,
  float *q2p,float *hp,float *gradvp,float *d2tp,int *i2,int *i3,int *i6, 
  int *map,float *vs,float *dvdxs,float *dvdys,float *dvdzs,
  float *uxx,float *uxy,float *uxz,float *uyy,float *uyz,float *uzz,float *tzt,
  float *xx,float *yy,float *zz,float *pxs,float *pys,float *pzs,
  float *e1xs,float *e1ys,float *e1zs,float *e2xs,float *e2ys,float *e2zs,
  float *p111s,float *p112s,float *p121s,float *p122s,float *q111s,
  float *q112s,float *q121s,float *q122s,float *p211s,float *p212s,
  float *p221s,float *p222s,float *q211s,float *q212s,float *q221s,float *q222s,
  float *dxx,float *dyy,float *dzz,float *dpxs,float *dpys,float *dpzs,
  float *de1x,float *de1y,float *de1z,float *de2x,float *de2y,float *de2z,
  float *dp111,float *dp112,float *dp121,float *dp122,float *dq111,
  float *dq112,float *dq121,float *dq122,float *dp211,float *dp212,
  float *dp221,float *dp222,float *dq211,float *dq212,float *dq221,float *dq222,
  float *xt,float *yt,float *zt,float *pxt,float *pyt,float *pzt,
  float *e1xt,float *e1yt,float *e1zt,float *e2xt,float *e2yt,float *e2zt,
  float *p111t,float *p112t,float *p121t,float *p122t,float *q111t,
  float *q112t,float *q121t,float *q122t,float *p211t,float *p212t,
  float *p221t,float *p222t,float *q211t,float *q212t,float *q221t,
  float *q222t,float *dxt,float *dyt,float *dzt,float *dpxt,float *dpyt,
  float *dpzt,float *de1xt,float *de1yt,float *de1zt,float *de2xt,
  float *de2yt,float *de2zt,float *dp111t,float *dp112t,float *dp121t,
  float *dp122t,float *dq111t,float *dq112t,float *dq121t,float *dq122t,
  float *dp211t,float *dp212t,float *dp221t,float *dp222t,float *dq211t,
  float *dq212t,float *dq221t,float *dq222t,
  float *vv, int *idisk,
  char *filevxx, char *filevxy, char *filevxz, 
  char *filevyy, char *filevyz, char *filevzz,
  int *lenvxx, int *lenvxy, int *lenvxz,
  int *lenvyy, int *lenvyz, int *lenvzz);
				   

void ov2int_(float *v,int *nx,int *ny,int *nz,float *fx,float *fy,float *fz,
  float *dx,float *dy,float *dz,float *ov2,int *nxo,int *nyo,int *nzo,
  float *fxo,float *fyo,float *fzo,float *dxo,float *dyo,float *dzo);
 
void dv2_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vxx,float *vxy,float *vxz,float *vyy,float *vyz,float *vzz);

void dv2xx_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vxx);
void dv2xy_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vxy);
void dv2xz_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vxz);
void dv2yy_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vyy);
void dv2yz_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vyz);
void dv2zz_(int *nx,int *ny,int *nz,float *dx,float *dy,float *dz,float *v,
  float *vzz);
 
void timeb_(int *nr,int *nz,float *dr,float *dz,float *fz,float *a,
  float *v0,float *t);

void resit_(float *t,int *nx,int *ny,int *nz,int *nr,float *dx,float *dy,
  float *dr,float *fx,float *fy,float *x0,float *y0,float *tb);

void interp_(int *nx,int *ny,int *nz,int *nout,float *fx,float *fy,
  float *dx,float *dy,float *x,float *y,float *t,float *tout);

void recot_(float *t,int *nout,int *nz,int *nr,float *dr,float *x,float *y,
  float *x0,float *y0,float *tb);


main(int argc, char **argv)
{
	int 	na,nt,nxs,nys,nxo,nyo,nzo,nx,ny,nz;
	int 	nxot,nyot,nzot,ixot,iyot,ixo,iyo,itemp,ixt,iyt;
	int 	ix,iy,nxtmax,nytmax,nzyxt;
	float   dt,xs,ys,fxs,fys,dxs,dys,exs,eys,
		fxo,fyo,fzo,dxo,dyo,dzo,exo,eyo,
		fa,amin,amax,da,azhmin,azhmax,
		fx,fy,fz,dx,dy,dz,ex,ey,ez,
		fac,tmax,aperx,apery,temp,
		*v,*vxx,*vxy,*vxz,*vyy,*vyz,*vzz,*vv,
		*vt,*vxxt,*vxyt,*vxzt,*vyyt,*vyzt,*vzzt,		
 		*t,*s,*ov2,*tt1,*tt2;
	float	fxot,fyot,fzot,dxot,dyot,dzot,xomin,xomax,yomin,yomax,
		xbmin,xbmax,xemin,xemax,yemin,yemax,ybmin,ybmax;
	float 	*ysp,*xsp,*azhnp,*azhxp,*fxtp,*fytp;
	int	*nxtp,*nytp,nzyx,nzyxo;
	int ixs,iys,ixs0,iys0,nsize,ek,ms,ncpu,np,is,is0,isnow,ip;
	int maxno,nev,iev,iout,*nout;
	float *xbeg,*xend,*ybeg,*yend,*dxout,*dyout,*xout,*yout,*tout;
	float v0,dvz,*tb,dr,rmax;
	int nr;
 	char *vfile, *tfile, *jpfile;
	string *tofile;
 	char *restart; 	
	int isres;
 	FILE *vfp, *tfp,**tofp,*jpfp;
	char *diskvxx,*diskvxy,*diskvxz,*diskvyy,*diskvyz,*diskvzz;
	char *filevxx,*filevxy,*filevxz,*filevyy,*filevyz,*filevzz;
	int lenvxx=0,lenvxy=0,lenvxz=0,lenvyy=0,lenvyz=0,lenvzz=0;
	FILE *diskfp;
	int idisk,idiskc;

	float gmin, gmax;
	int itmp;
  
        usghed ugh;
	ghed gh;
	int ierr; 
/*	int dtype, n1, n2, n3, n4, n5;
	float d1,d2,d3,d4,d5,o4,o5,dcdp2,dline3,ocdp2,oline3;*/

	long long isize;
	char *envs;
  
  float *xp,*yp,*zp,*pxp,*pyp,*pzp,*e1xp,*e1yp,*e1zp,*e2xp,*e2yp,*e2zp;
  float *q111p,*q112p,*q121p,*q122p,*p211p;
  float *p212p,*p221p,*p222p,*q211p,*q212p;
  float *q221p,*q222p,*vp,*dvdxp,*dvdyp,*dvdzp;
  int *nrsp;
  float *a0p,*azh0p;
  int n1,n2;
  float *p2p,*q2p,*hp,*gradvp,*d2tp;
  int i2,i3,i6,*map;
  float *vs,*dvdxs,*dvdys,*dvdzs,*uxx,*uxy,*uxz,*uyy,*uyz,*uzz,*tzt;
  float *xx,*yy,*zz,*pxs,*pys,*pzs;
  float *e1xs,*e1ys,*e1zs,*e2xs,*e2ys,*e2zs;
  float *p111s,*p112s,*p121s,*p122s,*q111s;
  float *q112s,*q121s,*q122s,*p211s,*p212s;
  float *p221s,*p222s,*q211s,*q212s,*q221s,*q222s;
  float *dxx,*dyy,*dzz,*dpxs,*dpys,*dpzs;
  float *de1x,*de1y,*de1z,*de2x,*de2y,*de2z;
  float *dp111,*dp112,*dp121,*dp122,*dq111;
  float *dq112,*dq121,*dq122,*dp211,*dp212;
  float *dp221,*dp222,*dq211,*dq212,*dq221,*dq222;
  float *xt,*yt,*zt,*pxt,*pyt,*pzt;
  float *e1xt,*e1yt,*e1zt,*e2xt,*e2yt,*e2zt;
  float *p111t,*p112t,*p121t,*p122t,*q111t;
  float *q112t,*q121t,*q122t,*p211t,*p212t;
  float *p221t,*p222t,*q211t,*q212t,*q221t;
  float *q222t,*dxt,*dyt,*dzt,*dpxt,*dpyt;
  float *dpzt,*de1xt,*de1yt,*de1zt,*de2xt;
  float *de2yt,*de2zt,*dp111t,*dp112t,*dp121t;
  float *dp122t,*dq111t,*dq112t,*dq121t,*dq122t;
  float *dp211t,*dp212t,*dp221t,*dp222t,*dq211t;
  float *dq212t,*dq221t,*dq222t;

  long long out_bytes;
  float compress, null=999999.;
  int num_points[3],nxyzot;
  long long work_bytes;

  float *in_data;

/* These variables are associated with setting up the threads. */
  int thread_count;
  pthread_t thread_id[MAXCPUS];
  compress_ttc_call_struct call_struct[MAXCPUS];
  pthread_attr_t attr;
	
/* hook up getpar to handle the parameters */
initargs(argc,argv);
requestdoc(1);
		
	/* get velocity information from header file */
	if( !getparstring("vfile",&vfile) ) {
		vfp = stdin;
	} else {
		vfp = efopen(vfile,"r"); 
	} 
	ierr = fgetusghdr(vfp,&ugh);
	if(ierr!=0) err("fgetusghdr error");
	nx = ugh.n3;
	ny = ugh.n2;
	nz = ugh.n1;
	fx = ugh.o3;
	fy = ugh.o2;
	fz = ugh.o1;
	dx = ugh.d3;
	dy = ugh.d2;
	dz = ugh.d1;
	if(nx<3 || ny<3 || nz<3 )
err("number of velocity samples in each direction must be not less than 3!\n");

	ex = fx+(nx-1)*dx;
	ey = fy+(ny-1)*dy;
	ez = fz+(nz-1)*dz;

	/* get optional parameters */
	if (!getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	}else {
		jpfp = fopen(jpfile,"w");
	}
	if (!getparint("nt",&nt)) nt = 401;
	if (!getparint("ncpu",&ncpu)) ncpu = 1;
	/* use ncpu=1 now */
	/*
	if(ncpu>1) { ncpu = 1; fprintf(jpfp,"ncpu reset to 1 \n"); }
	*/
	assert(ncpu<=MAXCPUS);
	envs = (char*) emalloc(80*sizeof(char));
        if(!getenv("PARALLEL")) {
                sprintf(envs,"%s=%d","PARALLEL",ncpu);
                putenv(envs);
		/* free(envs); */
        }
	if(nt>1001) err("nt cannot exceed 1001!\n");
	if (!getparfloat("dt",&dt)) dt = 0.008;
	tmax = (nt-1)*dt; 
	if (dt<0.000001) err("dt must be positive!\n");
 	if (!getparint("nyo",&nxot)) nxot = nx;
	if (!getparint("nxo",&nyot)) nyot = ny;
	if (!getparint("nzo",&nzot)) nzot = nz;
	if (!getparfloat("fyo",&fxot)) fxot = fx;
	if (!getparfloat("fxo",&fyot)) fyot = fy;
	if (!getparfloat("fzo",&fzot)) fzot = fz;
	if (!getparfloat("dyo",&dxot)) dxot = dx;
	if (!getparfloat("dxo",&dyot)) dyot = dy;
 	if (!getparfloat("dzo",&dzot)) dzot = dz;

	nzo = nzot;
	fzo = fzot;
	dzo = dzot;

	if (!getparint("ntline",&nev)) nev = 0;
	if(nev){
	    if (countparname("nout")!=nev)
		err("a nout array must be specified for each event");
	    if (countparname("xbeg")!=nev)
		err("a xbeg array must be specified for each event");
	    if (countparname("xend")!=nev)
		err("a xend array must be specified for each event");
	    if (countparname("ybeg")!=nev)
		err("a ybeg array must be specified for each event");
	    if (countparname("yend")!=nev)
		err("a yend array must be specified for each event");
	    if (countparname("tofile")!=nev)
		err("a tofile array must be specified for each event");

	    maxno = 0;