triseis.c

该程序主要用于三角网

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

/* TRISEIS: $Test Release: 1.1 $ ; $Date: 1997/10/29 17:23:49 $	*/

#include "par.h"
#include "Triangles/tri.h"
#include "Triangles/sloth.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" TRISEIS - Gaussian beam synthetic seismograms for a sloth model	",
"									",
"  triseis <modelfile >seisfile xs= zs= xg= zg= [optional parameters]	",
"									",
" Required Parameters:							",
" xs=            x coordinates of source surface			",
" zs=            z coordinates of source surface			",
" xg=            x coordinates of receiver surface			",
" zg=            z coordinates of receiver surface			",
"									",
" Optional Parameters:							",
" ns=1           number of sources uniformly distributed along s surface",
" ds=            increment between source locations (see notes below)	",
" fs=0.0         first source location (relative to start of s surface)	",
" ng=101         number of receivers uniformly distributed along g surface",
" dg=            increment between receiver locations (see notes below)	",
" fg=0.0         first receiver location (relative to start of g surface)",
" dgds=0.0       change in receiver location with respect to source location",
" krecord=1      integer index of receiver surface (see notes below)	",
" kreflect=-1    integer index of reflecting surface (see notes below)	",
" prim           =1, only single-reflected rays are considered 		",     
"                =0, only direct hits are considered  			",
" bw=0           beamwidth at peak frequency				",
" nt=251         number of time samples					",
" dt=0.004       time sampling interval					",
" ft=0.0         first time sample					",
" nangle=101     number of ray takeoff angles				",
" fangle=-45     first ray takeoff angle (in degrees)			",
" langle=45      last ray takeoff angle (in degrees)			",
" reftrans=0     =1 complex refl/transm. coefficients considered 	",
" atten=0        =1 add noncausal attenuation 				",
"                =2 add causal attenuation 				",
" lscale=        if defined restricts range of extrapolation		",
" fpeak=0.1/dt   peak frequency of ricker wavelet 			",
" aperture=      maximum angle of receiver aperture 			",
"									",
" NOTES:								",
" Only rays that terminate with index krecord will contribute to the	",
" synthetic seismograms at the receiver (xg,zg) locations.  The		",
" source and receiver locations are determined by cubic spline		",
" interpolation of the specified (xs,zs) and (xg,zg) coordinates.	",
" The default source location increment (ds) is determined to span	",
" the source surface defined by (xs,zs).  Likewise for dg.		",
"									",
NULL};

/*
 * AUTHOR:  Dave Hale, Colorado School of Mines, 02/09/91
 * MODIFIED:  Andreas Rueger, Colorado School of Mines, 08/18/93
 *	Modifications include: 2.5-D amplitudes, correction for ref/transm,
 *			timewindow, lscale, aperture, beam width, etc.
 */
/**************** end self doc ***********************************/

#define TAPER 0.97
#define C 0.99
#define EPS 0.0001

/* parameters for one step of dynamic ray tracing */
typedef struct RSStruct {
	float sigma,x,z,px,pz,t;
	float q1,p1,q2,p2;
        float atten;
        float ampli,ampliphase;
	int kmah,nref;
	EdgeUse *eu;
	Face *f;
} RayStep;

/* prototypes for functions defined and used internally */
static void makexzs (int nu, float *xu, float *zu, 
	int ns, float ds, float fs, float *xs, float *zs);
static void makesyn (float fpeak,int prim,int nre, RayEnd *re, float lscale,
	int krecord, int ng, float *xg, float *zg,
	float bw, int nt, float dt, float ft,
        int natten, int reftrans, float aperture);
static complex cricker (float w, float wpeak, float delay);
static int insideModel (Model *m, float x, float z);
static void shootRays (Model *m, float xs, float zs,
	int nangle, float dangle, float fangle,
	int kstop, int kreflect,
	RayEnd re[]);
static void traceRayInTri (RayStep *rs, RayStep *rsnew);
static int traceRayAcrossEdge (RayStep *rs, RayStep *rsnew);
static void reflectRayFromEdge (RayStep *rs, RayStep *rsnew);
static void evaluateDynamic (float sigma, 
	float px, float pz, float dsdx, float dsdz,
	float *q1, float *p1, float *q2, float *p2, int *kmah);
int checkIfSourceOnEdge (Face *tris, float zs, float xs);

/* the main program */
int main (int argc, char **argv)
{
	int nxs,nzs,nxzs,nxg,nzg,nxzg,
		ns,ng,krecord,kreflect,nt,nangle,it,is,ig;
	int natten,reftrans,prim;
	float ds,fs,dg,fg,dgds,s,bw,dt,ft,cpuray,cpuseis,
		fangle,langle,dangle,lscale,fpeak,aperture,
		*xsu,*zsu,*xgu,*zgu,
		*xs,*zs,*xg,*zg,*zeros;
	Model *m;	
	RayEnd *re;

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

	/* read sloth model */
	m = readModel(stdin);

	/* source location parameters */
	nxs = countparval("xs");
	nzs = countparval("zs");
	if (nxs==0 || nzs==0) err("must specify both xs and zs");
	if (nxs!=nzs) err("number of xs must equal number of zs");
	nxzs = nxs;
	xsu = ealloc1float(nxzs);
	zsu = ealloc1float(nxzs);
	getparfloat("xs",xsu);
	getparfloat("zs",zsu);
	if (!getparint("ns",&ns)) ns = 1;
	if (!getparfloat("ds",&ds)) ds = 0.0;
	if (!getparfloat("fs",&fs)) fs = 0.0;

	/* receiver location parameters */
	nxg = countparval("xg");
	nzg = countparval("zg");
	if (nxg==0 || nzg==0) err("must specify both xg and zg");
	if (nxg!=nzg) err("number of xg must equal number of zg");
	nxzg = nxg;
	xgu = ealloc1float(nxzg);
	zgu = ealloc1float(nxzg);
	getparfloat("xg",xgu);
	getparfloat("zg",zgu);
	if (!getparint("ng",&ng)) ng = 101;
	if (!getparfloat("dg",&dg)) dg = 0.0;
	if (!getparfloat("fg",&fg)) fg = 0.0;
	if (!getparfloat("dgds",&dgds)) dgds = 0.0;
	if (!getparint("krecord",&krecord)) krecord = 1;
	if (!getparint("prim",&prim)) prim = INT_MAX;

	/* other parameters */
	if (!getparint("kreflect",&kreflect)) kreflect = -1;
	if (!getparfloat("bw",&bw)) bw = 0.0;
	if (!getparint("nt",&nt)) nt = 251;
	if (!getparfloat("dt",&dt)) dt = 0.004;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("fpeak",&fpeak)) fpeak = 0.1/dt;

	if (0.5/dt<2.0*fpeak) {
		fprintf(stderr,"WARNING: ALIASING POSSIBLE\n");
		fprintf(stderr,"decrease dt or reduce fpeak\n");
	}

	if (!getparint("nangle",&nangle)) nangle = 101;
	if (!getparfloat("fangle",&fangle)) fangle = -45.0;
	if (!getparfloat("langle",&langle)) langle = 45.0;
	if (!getparint("reftrans",&reftrans)) reftrans = 0;
	if (!getparint("atten",&natten)) natten = 0;
	if (!getparfloat("aperture",&aperture)) aperture = FLT_MAX;
	if (!getparfloat("lscale",&lscale)) lscale = FLT_MAX;
	if (aperture != FLT_MAX) aperture = asin(aperture*PI/180.0);

	/* convert angles to radians and determine angle increment */
	fangle *= PI/180.0;
	langle *= PI/180.0;
	dangle = (nangle>1) ? (langle-fangle)/(nangle-1) : 1.0;

	/* allocate space */
	xs = ealloc1float(ns);
	zs = ealloc1float(ns);
	xg = ealloc1float(ng);
	zg = ealloc1float(ng);
	zeros = ealloc1float(nt);
	re = ealloc1(nangle,sizeof(RayEnd));

	/* determine source locations */
	makexzs(nxzs,xsu,zsu,ns,ds,fs,xs,zs);

	/* make trace with all zeros */
	for (it=0; it<nt; ++it)
		zeros[it] = 0.0;

	/* loop over source locations */
	for (is=0,s=fs,cpuray=cpuseis=0.0; is<ns; ++is,s+=ds) {

		/* uncomment for diagnostic print */
		fprintf(stderr,"shot %3d of %3d at x=%g\n",is,ns,xs[is]);
/*		fprintf(stderr,"shot %d at x=%.25g,z=%.25g\n",
			is,xs[is],zs[is]);
		fprintf(stderr,"rec %d at x=%.25g\n",is,fg+dgds*(s-fs));

		if (is<32 || is>33) continue;
*/
		/* if source is outside of model boundaries */
		if (!insideModel(m,xs[is],zs[is])) {

			/* write zeros */
			for (ig=0; ig<ng; ++ig)
				if (efwrite(zeros,sizeof(float),nt,stdout)!=nt)
					err("error writing output traces!\n");

			/* continue with next source location */
			continue;
		}

		/* shoot rays */
		/* cputemp = cpusec();*/
		shootRays(m,xs[is],zs[is],nangle,dangle,
				fangle,krecord,kreflect,re);
		/* cpuray += cpusec()-cputemp;*/

		/* compute receiver locations */
		makexzs(nxzg,xgu,zgu,ng,dg,fg+dgds*(s-fs),xg,zg);

		/* compute and write synthetic seismograms */
		/* cputemp = cpusec(); */
		makesyn(fpeak,prim,nangle,re,lscale,krecord,						
			ng,xg,zg,bw,nt,dt,ft,natten,reftrans,aperture);
	/*	cpuseis += cpusec()-cputemp; */
	}
	cpuray += 0.0;
	cpuseis += 0.0;

	/* cpu timings */
	/*
         * fprintf(stderr,"\ntriseis:  total cpu time = %g s\n",cpusec());
         * fprintf(stderr,"triseis:  ray tracing cpu time = %g s\n",cpuray);
         * fprintf(stderr,"triseis:  seismogram cpu time = %g s\n",cpuseis);
	 */

	return EXIT_SUCCESS;
}

static void makexzs (int nu, float *xu, float *zu, 
	int ns, float ds, float fs, float *xs, float *zs)
/* interpolate (x,z) coordinates uniformly sampled in distance s */
{
	int iu,nuu,iuu,is,js;
	float x,z,xlast,zlast,dx,dz,duu,uu,temp,
		*u,*s,(*xud)[4],(*zud)[4],*us;

	xud = (float(*)[4])ealloc1float(4*nu);
	zud = (float(*)[4])ealloc1float(4*nu);
	u = ealloc1float(nu);
	for (iu=0; iu<nu; ++iu)
		u[iu] = iu;
	csplin(nu,u,xu,xud);
	csplin(nu,u,zu,zud);
	nuu = 10*nu;
	duu = (u[nu-1]-u[0])/(nuu-1);
	s = ealloc1float(nuu);
	s[0] = 0.0;
	xlast = xu[0];
	zlast = zu[0];
	for (iuu=0,uu=0.0,s[0]=0.0; iuu<nuu; ++iuu,uu+=duu) {
		intcub(0,nu,u,xud,1,&uu,&x);
		intcub(0,nu,u,zud,1,&uu,&z);
		dx = x-xlast;
		dz = z-zlast;
		s[iuu] = s[iuu-1]+sqrt(dx*dx+dz*dz);
		xlast = x;
		zlast = z;
	}			
	us = ealloc1float(ns);
	if (ds==0.0) ds = (s[nuu-1]-s[0])/((ns>1)?ns-1:1);
	if (ds<0.0) {
		yxtoxy(nuu,duu,0.0,s,ns,-ds,fs+(ns-1)*ds,0.0,(float)(nu-1),us);
		for (is=0,js=ns-1; is<js; ++is,--js) {
			temp = us[is];
			us[is] = us[js];
			us[js] = temp;
		}
	} else {
		yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu-1),us);
	}
	intcub(0,nu,u,xud,ns,us,xs);
	intcub(0,nu,u,zud,ns,us,zs);
	free1float(us);
	free1float(s);
	free1float(u);
	free1float((float*)xud);
	free1float((float*)zud);
}

static void makesyn (float fpeak, int prim, int nre, RayEnd *re, float lscale,
	int krecord, int ng, float *xg, float *zg, float bw,
	int nt, float dt, float ft, int natten, int reftrans, float aperture)
/* make and write synthetic seismograms */
{
	int ntfft,nw,iw,ig,ire,kend,kmah,nref;
	float dw,w,wpeak,delay,scale,wref,eps1,eps2,n,
		t,x,z,px,pz,q1,p1,q2,p2,v0,v,dsdx,dsdz,
		dx,dz,delt,dels,nn,phase,ta,cpoqr,cpoqi,disq,
		*syn,sigma,*lnvec,sqlscale,dangle;
        float ampli,ampliphase,atten,temp,c,eps,v2,v4,dvds,dvdn;
	double campr,campi,cosw,sinw,expw,cosd,sind,expd,tempd;
	complex ceps,cq,cp,cpoq,cqr,camp,
		*cwave,**csyn;
	
	/* constants */
	ntfft = npfaro(nt,2*nt);
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	wpeak = 2*PI*fpeak;
	wref = 1000.0;
	delay = 0.0;
	sqlscale = lscale*lscale;
	c = C;
	eps = EPS;

	/* allocate workspace */
	syn = ealloc1float(ntfft);
	cwave = ealloc1complex(nw);
	csyn = ealloc2complex(nw,ng);
	lnvec = ealloc1float(nw);

	/* initialize synthetics */
	for (ig=0; ig<ng; ++ig)
		for (iw=0; iw<nw; ++iw)
			csyn[ig][iw] = cmplx(0.0,0.0);

	/* precomputing causal attenuation term */
	if (natten==2) {
		temp = 1.0/PI;
		lnvec[0] = -1.0;
		for(iw=1,w=dw; iw<nw; ++iw,w+=dw)
			lnvec[iw] = log(w/wref)*temp;
	}

	/* loop over rays */
	for (ire=0; ire<nre; ++ire) {

		/* determine index of ray end */
		kend = re[ire].kend;

		/* if not the end index we want, skip this ray */
		if (kend!=krecord) continue;

		nref = re[ire].nref;

		/* if ray was not reflected */
		if (nref!=prim && prim!=INT_MAX) continue;

		px = re[ire].px;
		v = 1.0/sqrt(re[ire].se);

		/* if not within aperture, skip this ray */
		if (ABS(px*v)>aperture) continue;

		/* ray end parameters */
		sigma = 1/sqrt(re[ire].sigma);
		x = re[ire].x;
		z = re[ire].z;
		t = re[ire].t;
		pz = re[ire].pz;
		q1 = re[ire].q1;
		p1 = re[ire].p1;
		q2 = re[ire].q2;
		p2 = re[ire].p2;
		kmah = re[ire].kmah;
		atten = re[ire].atten;
		ampli = re[ire].ampli;
		ampliphase = re[ire].ampliphase;
		dangle = re[ire].dangle;
		dsdx = re[ire].dsdxe;
		dsdz = re[ire].dsdze;
		v0 = 1.0/sqrt(re[ire].sb);

		/* compute ray dependant constants */
		v2 = v*v;
		v4 = v2*v2;
		dvds = -0.5*v4*(dsdx*px+dsdz*pz);
		dvdn = -0.5*v4*(dsdx*pz-dsdz*px);

		scale = sigma*dangle*sqrt(v/(v0*2*PI));

		/* complex beam epsilon */
		if (bw!=0.0) {
			eps1 = 0.0;
			eps2 = -0.5*wpeak*bw*bw;
		} else {
			temp = c * ABS(p2/q2) + eps;
			eps1 = -q2/q1 * (p2*p1/(q2*q1) + temp*temp);
			eps1 = eps1/(p1*p1/(q1*q1) + temp*temp);
			eps2 = -temp/(p1*p1+q1*q1*temp*temp);		
		}
		ceps = cmplx(eps1,eps2);

		/* complex p, q, and p/q */
		cp = crmul(ceps,p1);  cp.r += p2;
		cq = crmul(ceps,q1);  cq.r += q2;
		cpoq = cdiv(cp,cq);
		cpoqr = cpoq.r;
		cpoqi = cpoq.i;

		/* if cpoqi negative, skip ray */
		/* (cpoqi<0 leads to exponential growth!) */
		if (cpoqi<0.0) {
			continue;
		}

		/* complex q, adjusted for reflections */
		cqr = (nref%2) ? cneg(cq) : cq;

		/* frequency- and ray-independent amplitude factor */
		phase = -PI*((kmah+1)/2);
		camp = csqrt(cdiv(ceps,cqr));

		/* correction for cmplex reflection/transmission */
		if (reftrans) {
			phase += ampliphase;
			camp = crmul(camp,ampli);
		}

                /* correction for 2.5-D Spreading */
        	phase += PI/4.0;

                camp = crmul(camp,scale);
		camp = cmul(camp,cmplx(cos(phase),sin(phase)));
		campr = camp.r;
		campi = camp.i;

		/* save computing time */
     		if (natten == 2)
			for (iw=0; iw<nw; ++iw)
				lnvec[iw] = atten*lnvec[iw];

		/* loop over receiver locations */
		for (ig=0; ig<ng; ++ig) {

			/* delta x and delta z */
			dx = xg[ig]-x;
			dz = zg[ig]-z;
			disq = dx*dx+dz*dz;

			/* delta t and delta s */
			n = v*(pz*dx-px*dz);
			nn = n*n;

			/* do not extrapolate too far */
			if (lscale!=FLT_MAX && disq>sqlscale) continue;

			dels = v*(px*dx+pz*dz)*(1-dvdn/v*n);
			delt = dels/v-0.5*dvds*dels*dels/v2;