susynvxzcs.c

su 的源代码库

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

/* SUSYNVXZCS: $Revision: 1.13 $ ; $Date: 2006/11/07 22:58:42 $	*/

#include "su.h" 
#include "segy.h" 

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUSYNVXZCS - SYNthetic seismograms of common shot in V(X,Z) media via	",
" 		Kirchhoff-style modeling				",
" 									",
" susynvxzcs<vfile >outfile  nx= nz= [optional parameters]		",
" 									",
" Required Parameters:							",
" <vfile        file containing velocities v[nx][nz]			",
" >outfile      file containing seismograms of common ofset		",
" nx=           number of x samples (2nd dimension) in velocity ",
" nz=           number of z samples (1st dimension) in velocity ",
" 									",
" Optional Parameters:							",
" nt=501        	number of time samples				",
" dt=0.004      	time sampling interval (sec)			",
" ft=0.0        	first time (sec)				",
" fpeak=0.2/dt		peak frequency of symmetric Ricker wavelet (Hz)	",
" nxg=			number of receivers of input traces		",
" dxg=15		receiver sampling interval (m)			",
" fxg=0.0		first receiver (m)				",
" nxd=5         	skipped number of receivers			",
" nxs=1			number of offsets				",
" dxs=50		shot sampling interval (m)			",
" fxs=0.0		first shot (m)				",
" dx=50         	x sampling interval (m)				",
" fx=0.         	first x sample (m)				",
" dz=50         	z sampling interval (m)				",
" nxb=nx/2    	band width centered at midpoint (see note)	",
" nxc=0         hozizontal range in which velocity is changed	",
" nzc=0         vertical range in which velocity is changed	",
" pert=0        =1 calculate time correction from v_p[nx][nz]	",
" vpfile        file containing slowness perturbation array v_p[nx][nz]	",
" ref=\"1:1,2;4,2\"	reflector(s):  \"amplitude:x1,z1;x2,z2;x3,z3;...\"",
" smooth=0		=1 for smooth (piecewise cubic spline) reflectors",
" ls=0			=1 for line source; =0 for point source		",
" tmin=10.0*dt		minimum time of interest (sec)			",
" ndpfz=5		number of diffractors per Fresnel zone		",
" verbose=0		=1 to print some useful information		",
" 									",
" Notes:								",
" This algorithm is based on formula (58) in Geo. Pros. 34, 686-703,	",
" by N. Bleistein.							",
"									",
" Traveltime and amplitude are calculated by finite difference which	",
" is done only in one of every NXD receivers; in skipped receivers, 	",
" interpolation is used to calculate traveltime and amplitude.		", 
" For each receiver, traveltime and amplitude are calculated in the 	",
" horizontal range of (xg-nxb*dx, xg+nxb*dx). Velocity is changed by 	",
" constant extropolation in two upper trianglar corners whose width is 	",
" nxc*dx and height is nzc*dz.						",
"									",
" Eikonal equation will fail to solve if there is a polar turned ray.	",
" In this case, the program shows the related geometric information. 	",
" There are three ways to remove the turned rays: smoothing velocity, 	",
" reducing nxb, and increaing nxc and nzc (if the turned ray occurs  	",
" in shallow areas). To prevent traveltime distortion from an over-	",
" smoothed velocity, traveltime is corrected based on the slowness 	",
" perturbation.								",
"									",
" More than one ref (reflector) may be specified.			",
" Note that reflectors are encoded as quoted strings, with an optional	",
" reflector amplitude: preceding the x,z coordinates of each reflector.	",
" Default amplitude is 1.0 if amplitude: part of the string is omitted.	",
"									",
NULL};

/*
 *	Author: Zhenyue Liu, 07/20/92, Center for Wave Phenomena
 *		Many subroutines borrowed from Dave Hale's program: SUSYNLV
 *
 *		Trino Salinas, 07/30/96, fixed a bug in the geometry
 *		setting to allow the spread move with the shots.
 *
 * Trace header fields set: trid, counit, ns, dt, delrt,
 *				tracl. tracr, fldr, tracf,
 *				sx, gx
 */
/**************** end self doc ***********************************/

/* global varibles for detecting error in solving eikonal equation */
int ierr=0;
float x_err, z_err, r_err;
int ia_err;


/* parameters for half-derivative filter */
#define LHD 20
#define NHD 1+2*LHD

static void makeone (float **ts, float **as, float **sgs, float **tg, 
	float **ag, float **sgg, float ex, float ez, float dx,  
	float dz, float fx, float vs0, float vg0, int ls, Wavelet *w,
	int nr, Reflector *r, int nt, float dt, float ft, float *trace);

/* additional prototypes for eikonal equation functions */
void delta_t (int na, float da, float r, float dr, 
	float uc[], float wc[], float sc[], float bc[],
	float un[], float wn[], float sn[], float bn[]);
void eiktam_pert (float xs, float zs, int nz, float dz, float fz, int nx, 
	float dx, float fx, float **vel, float **vel_p, int pert, 
	float **time, float **angle, float **sig, float **bet);

/* segy trace */
segy tr;

int
main (int argc, char **argv)
{
	int 	nr,ir,ls,smooth,ndpfz,ns,ixs,ixg,nxs,nxg,nt,
		nx,nz,nxb,nxd,ixd,i,ix,iz,nx1,nx0,nxd1,nxc,nzc,
		verbose,tracl,pert,
		*nxz;
	float   tmin,temp,temp1,
		dsmax,fpeak,dx,dz,fx,ex,fx1,ex1,
		dxg,dxs,dt,fxg,fxs,ft,xs,xg,vs0,vg0,ez,
		*ar,**xr,**zr,
		**vold,**vpold=NULL,**told,**aold,**sgold,
		**ts,**as,**sgs,**tg,**ag,**sgg,**bas,**bag,
		**v,**vp=NULL,**tg1=NULL,**ag1=NULL,**sgg1=NULL;
	FILE *vfp=stdin;
	char *vpfile="";
	Reflector *r;
	Wavelet *w;

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);

	/* get required parameters */
	if (!getparint("nx",&nx)) err("must specify nx!\n");
	if (!getparint("nz",&nz)) err("must specify nz!\n");
	
	/* get optional parameters */
	if (!getparint("nt",&nt)) nt = 501; CHECK_NT("nt",nt);
	if (!getparfloat("dt",&dt)) dt = 0.004;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("fpeak",&fpeak)) fpeak = 0.2/dt;
	if (!getparint("nxg",&nxg)) nxg = 101;
	if (!getparfloat("dxg",&dxg)) dxg = 15;
	if (!getparfloat("fxg",&fxg)) fxg = 0.0;
	if (!getparint("nxd",&nxd)) nxd = 5;
	if (!getparint("nxs",&nxs)) nxs = 1;
	if (!getparfloat("dxs",&dxs)) dxs = 50;
	if (!getparfloat("fxs",&fxs)) fxs = 0.0;
	if (!getparfloat("dx",&dx)) dx = 50;
	if (!getparfloat("fx",&fx)) fx = 0.0;
	if (!getparfloat("dz",&dz)) dz = 50;
	if (!getparint("nxb",&nxb)) nxb = nx/2;
	if (!getparint("nxc",&nxc)) nxc = 0;
	nxc = MIN(nxc,nxb);
	if (!getparint("nzc",&nzc)) nzc = 0;
	nzc = MIN(nzc,nz);
	if (!getparint("ls",&ls)) ls = 0;
	if (!getparint("pert",&pert)) pert = 0;
	if (!getparfloat("tmin",&tmin)) tmin = 10.0*dt;
	if (!getparint("ndpfz",&ndpfz)) ndpfz = 5;
	if (!getparint("smooth",&smooth)) smooth = 0;
	if (!getparint("verbose",&verbose)) verbose = 0;
	
	/* check the ranges of shots and receivers */
	ex = fx+(nx-1)*dx;
	ez = (nz-1)*dz;
	for (ixs=0; ixs<nxs; ++ixs) {
		/* compute shot coordinates */
		xs = fxs+ixs*dxs;
		if (fx>xs || ex<xs) 
		err("shot lies outside of specified (x,z) grid\n");
	}

	for (ixs=0; ixs<2; ++ixs) {
		for (ixg=0; ixg<nxg; ++ixg) {
			/* compute receiver coordinates */
			xg = fxg+ixg*dxg+ixs*(nxs-1)*dxs;
			if (fx>xg || ex<xg)
			err("receiver lies outside of specified (x,z) grid\n");
		}
	}


	decodeReflectors(&nr,&ar,&nxz,&xr,&zr);
	if (!smooth) breakReflectors(&nr,&ar,&nxz,&xr,&zr);

	/* allocate space */
	vold = ealloc2float(nz,nx);
	aold = ealloc2float(nz,nx);
	told = ealloc2float(nz,nx);
	bas = ealloc2float(nz,nx);
	sgold = ealloc2float(nz,nx);
	if(pert) vpold = ealloc2float(nz,nx);
	/* read velocities and slowness perturbation*/
	if(fread(vold[0],sizeof(float),nx*nz,vfp)!=nx*nz)
		err("cannot read %d velocities from file %s",nx*nz,vfp);
	if(pert) {
		/* read slowness perturbation*/
		if (!getparstring("vpfile",&vpfile)) 
			err("must specify vpfile!\n");
		if(fread(vpold[0],sizeof(float),nx*nz,fopen(vpfile,"r"))
			!=nx*nz)
		err("cannot read %d slowness perturbation from file %s"
			,nx*nz,vpfile);
		/* calculate 1/2 perturbation of slowness squares*/
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				vpold[ix][iz] = vpold[ix][iz]/vold[ix][iz];
	}
	/* determine maximum reflector segment length */
	tmin = MAX(tmin,MAX(ft,dt));
	dsmax = vold[0][0]/(2*ndpfz)*sqrt(tmin/fpeak);
 	
	/* make reflectors */
	makeref(dsmax,nr,ar,nxz,xr,zr,&r);

	/* count reflector segments */
	for (ir=0,ns=0; ir<nr; ++ir)
		ns += r[ir].ns;

	/* make wavelet */
	makericker(fpeak,dt,&w);
	
	/* if requested, print information */
	if (verbose) {
		warn("\nSUSYNVXZCS:");
		warn("Total number of small reflecting segments is %d.\n",ns);
	}
	
	/* set constant segy trace header parameters */
	memset((void *) &tr, 0, sizeof(segy));
	tr.trid = 1;
	tr.counit = 1;
	tr.ns = nt;
	tr.dt = 1.0e6*dt;
	tr.delrt = 1.0e3*ft;
	tr.d2 = dxg;
	
	/* allocate space */
	nx1 = 1+2*nxb;
	ts = ealloc2float(nz,nx1);
	as = ealloc2float(nz,nx1);
	sgs = ealloc2float(nz,nx1);
	tg = ealloc2float(nz,nx1);
	ag = ealloc2float(nz,nx1);
	sgg = ealloc2float(nz,nx1);
	v = ealloc2float(nz,nx1);
	bag = ealloc2float(nz,nx1);
	if(pert) vp = ealloc2float(nz,nx1);
	if(nxd>1) {
		/* allocate space for interpolation */
		tg1 = ealloc2float(nz,nx1);
		ag1 = ealloc2float(nz,nx1);
		sgg1 = ealloc2float(nz,nx1);
  	}

	/* change velocity values within two trianglar upper corners to
 		avoid possible polar-turned rays	*/
	for(iz=0; iz<nzc; ++iz)
		for(ix=iz*nxc/nzc; ix<nxc; ++ix){
			vold[nxc-ix-1][iz] = vold[nxc-iz*nxc/nzc][iz];
			vold[nx1-nxc+ix][iz] = vold[nx1-nxc+iz*nxc/nzc-1][iz];
		}

		
	/* loop over shots */
	for (ixs=0, tracl=0; ixs<nxs; ++ixs){
	    xs = fxs+ixs*dxs;
	    /* calculate traveltimes from shot */
	    eiktam_pert(xs,0,nz,dz,0,nx,dx,fx,vold,vpold,pert,told,aold,sgold,bas);
	    if(ierr){
		/* Ekonal equation suffers from turned rays*/
		err("\tEikonal equation fails to solve at "
		    "x=%g, z=%g\n\tfrom shot xs=%g .\n",x_err+xs,z_err,xs);
		}
	    /* calculate velocity at shot 	*/
	    temp = (xs-fx)/dx;
	    ixd = temp;
	    temp = temp-ixd;
	    vs0 = (1-temp)*vold[ixd][0]+temp*vold[ixd+1][0];
	
	    /* save the skipping number */
	    nxd1 = nxd;
	    
	    /* loop over receivers */
	    for (ixg=0; ixg<nxg; ixg +=nxd1){
		xg = fxg+ixs*dxs+ixg*dxg;

		/* set range for traveltimes' calculation */
		fx1 = xg-nxb*dx;
		nx0 = (fx1-fx)/dx;
		temp = (fx1-fx)/dx-nx0;
	       	for(ix=0;ix<nx1;++ix){
		    if(ix<-nx0) 
			for(iz=0;iz<nz;++iz)
				v[ix][iz] = vold[0][iz];
		    else if(ix+nx0>nx-2)
			for(iz=0;iz<nz;++iz)
				v[ix][iz]=vold[nx-1][iz];
		    else	
			for(iz=0;iz<nz;++iz){
				v[ix][iz] = vold[ix+nx0][iz]*(1.0-temp)
					+temp*vold[ix+nx0+1][iz];
				ts[ix][iz] = told[ix+nx0][iz]*(1.0-temp)
					+temp*told[ix+nx0+1][iz];
				as[ix][iz] = aold[ix+nx0][iz]*(1.0-temp)
					+temp*aold[ix+nx0+1][iz];
				sgs[ix][iz] = sgold[ix+nx0][iz]*(1.0-temp)
					+temp*sgold[ix+nx0+1][iz];
				}
		}
	if(pert){
	    for(ix=0;ix<nx1;++ix)
		for(iz=0;iz<nz;++iz){
		    if(ix<-nx0) 
			vp[ix][iz] = vpold[0][iz];
		    else if(ix+nx0>nx-2) 
			vp[ix][iz]=vpold[nx-1][iz];
		    else
			vp[ix][iz] = vpold[ix+nx0][iz]*(1.0-temp)
				+temp*vpold[ix+nx0+1][iz];
		    }
	}
			
		if(ixg==0 || nxd1==1){
		/* No interpolation */
	
			/* compute traveltimes, propagation angles, sigmas 
	  		 from receiver */
		eiktam_pert(xg,0.,nz,dz,0.,nx1,dx,fx1,v,vp,pert,tg,ag,sgg,bag);
		if(ierr){
			/* Ekonal equation fails to solve */
			warn("\tEikonal equation fails to solve at "
			     "x=%g, z=%g\n\tfrom receiver xg=%g .",
			     x_err+xg,z_err,xg);
			break;
		}
			/* calculate velocity at receiver 	*/
			vg0 = (1-temp)*vold[nx0+nxb][0]
				+temp*vold[nx0+nxb+1][0];

			/* make one trace */
			ex1 = MIN(ex,xg+nxb*dx);
			makeone(ts,as,sgs,tg,ag,sgg,ex1,ez,dx,dz,fx1,vs0,vg0,
				ls,w,nr,r,nt,dt,ft,tr.data);
			/* set segy trace header parameters */
			tr.tracl = tr.tracr = ++tracl;
			tr.sx = NINT(xs);
			tr.gx = NINT(xg);
			/* write trace */
			puttr(&tr);
		}
		else {
			/* Linear interpolation */
			
		eiktam_pert(xg,0,nz,dz,0,nx1,dx,fx1,v,vp,pert,tg1,ag1,sgg1,bag);
		if(ierr){
			/* Ekonal equation fails to solve */
			warn("\tEikonal equation fails to solve at "
			     "x=%g, z=%g\n\tfrom receiver xg=%g .",
			     x_err+xg,z_err,xg);
			break;
		}

		/* calculate velocity at receiver 	*/
		vg0 = (1-temp)*vold[nx0+nxb][0]+temp*vold[nx0+nxb+1][0];
			
			/* interpolate quantities between two midpoints	*/
		    	xg -= nxd1*dxg;
		    for(i=1; i<=nxd1; ++i) {
		    	xg += dxg;
			fx1 = xg-nxb*dx;	
			ex1 = xg+nxb*dx;
			temp = nxd1-i;
			temp1 = 1.0/(nxd1-i+1);
			for(ix=0;ix<nx1;++ix)
			    for(iz=0;iz<nz;++iz){
			    if(i==nxd1){
			   	tg[ix][iz] = tg1[ix][iz];
			   	ag[ix][iz] = ag1[ix][iz];
			   	sgg[ix][iz] = sgg1[ix][iz];
				}
			    else{
			   	tg[ix][iz] = (temp*tg[ix][iz]
					+tg1[ix][iz])*temp1;
			   	ag[ix][iz] = (temp*ag[ix][iz]
					+ag1[ix][iz])*temp1;
			   	sgg[ix][iz] = (temp*sgg[ix][iz]
					+sgg1[ix][iz])*temp1;
				}
			}

			/* calculate quantities at shot 	*/
	    		nx0 = (fx1-fx)/dx;
	    		temp = (fx1-fx)/dx-nx0;
	   		for(ix=MAX(0,-nx0); ix<MIN(nx1,nx-1-nx0); ++ix){
			    for(iz=0;iz<nz;++iz){
			   	ts[ix][iz] = told[ix+nx0][iz]*(1.0-temp)
					+temp*told[ix+nx0+1][iz];
			    	sgs[ix][iz] = sgold[ix+nx0][iz]*(1.0-temp)
					+temp*sgold[ix+nx0+1][iz];
			        as[ix][iz] = aold[ix+nx0][iz]*(1.0-temp)
					+temp*aold[ix+nx0+1][iz];
		   		 }
			}
				
			/* make one trace */
			ex1 = MIN(ex,xg+nxb*dx);
			makeone(ts,as,sgs,tg,ag,sgg,ex1,ez,dx,dz,fx1,vs0,vg0,
				ls,w,nr,r,nt,dt,ft,tr.data);
			/* set segy trace header parameters */
			tr.tracl = tr.tracr = ++tracl;
			tr.sx = NINT(xs);
			tr.gx = NINT(xg);
			/* write trace */
			puttr(&tr);
		    }