modeling.c

su 的源代码库

		nuu = 20*nu[ir];
		duu = (u[nu[ir]-1]-u[0])/(nuu-1);
		s = ealloc1float(nuu);
		s[0] = 0.0;
		xlast = xu[ir][0];
		zlast = zu[ir][0];
		for (iuu=1,uu=duu; iuu<nuu; ++iuu,uu+=duu) {
			intcub(0,nu[ir],u,xud,1,&uu,&x);
			intcub(0,nu[ir],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;
		}

		/* compute u(s) from s(u) */
		ns = 1+s[nuu-1]/dsmax;
		ds = s[nuu-1]/ns;
		fs = 0.5*ds;
		us = ealloc1float(ns);
		yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu[ir]-1),us);

		/* compute reflector segments uniformly sampled in s */
		rs = ealloc1(ns,sizeof(ReflectorSegment));
		for (is=0; is<ns; ++is) {
			intcub(0,nu[ir],u,xud,1,&us[is],&rsx);
			intcub(0,nu[ir],u,zud,1,&us[is],&rsz);
			intcub(1,nu[ir],u,xud,1,&us[is],&rsxd);
			intcub(1,nu[ir],u,zud,1,&us[is],&rszd);
			rs[is].x = rsx;
			rs[is].z = rsz;
			rs[is].c = rsxd/sqrt(rsxd*rsxd+rszd*rszd);
			rs[is].s = -rszd/sqrt(rsxd*rsxd+rszd*rszd);
		}
		
		/* fill in reflector structure */
		rr[ir].ns = ns;
		rr[ir].ds = ds;
		rr[ir].a = ar[ir];
		rr[ir].rs = rs;

		/* free workspace */
		free1float(us);
		free1float(s);
		free1float(u);
		free1float((float*)xud);
		free1float((float*)zud);

		/* free space replaced by reflector segments */
		free1(xu[ir]);
		free1(zu[ir]);
	}

	/* free space replaced by reflector segments */
	free1(nu);
	free1(xu);
	free1(zu);
}

void raylv2 (float v00, float dvdx, float dvdz,
	float x0, float z0, float x, float z,
	float *c, float *s, float *t, float *q)
/*****************************************************************************
Trace ray between two points, for linear velocity v = v00+dvdx*x+dvdz*z
******************************************************************************
Input:
v00		velocity at (x=0,z=0)
dvdx		derivative dv/dx of velocity with respect to x
dvdz		derivative dv/dz of velocity with respect to z
x0		x coordinate at beginning of ray
z0		z coordinate at beginning of ray
x		x coordinate at end of ray
z		z coordinate at end of ray

Output:
c		cosine of propagation angle at end of ray
s		sine of propagation angle at end of ray
t		time along raypath
q		integral of velocity along raypath
*****************************************************************************/
{
	float a,oa,v0,v,cr=0.0,sr=0.0,r,or,c0,mx,temp;
	
	/* if (x,z) same as (x0,z0) */
	if (x==x0 && z==z0) {
		*c = 1.0;
		*s = 0.0;
		*t = 0.0;
		*q = 0.0;;
		return;
	}

	/* if constant velocity */
	if (dvdx==0.0 && dvdz==0.0) {
		x -= x0;
		z -= z0;
		r = sqrt(x*x+z*z);
		or = 1.0/r;
		*c = z*or;
		*s = x*or;
		*t = r/v00;
		*q = r*v00;
		return;
	}

	/* if necessary, rotate coordinates so that v(x,z) = v0+a*(z-z0) */
	if (dvdx!=0.0) {
		a = sqrt(dvdx*dvdx+dvdz*dvdz);
		oa = 1.0/a;
		cr = dvdz*oa;
		sr = dvdx*oa;
		temp = cr*x0-sr*z0;
		z0 = cr*z0+sr*x0;
		x0 = temp;
		temp = cr*x-sr*z;
		z = cr*z+sr*x;
		x = temp;
	} else {
		a = dvdz;
	}

	/* velocities at beginning and end of ray */
	v0 = v00+a*z0;
	v = v00+a*z;
	
	/* if either velocity not positive */
	if (v0<0.0 || v<0.0) {
		*c = 1.0;
		*s = 0.0;
		*t = FLT_MAX;
		*q = FLT_MAX;
		return;
	}

	/* translate (x,z) */
	x -= x0;
	z -= z0;
	
	/* if raypath is parallel to velocity gradient */
	if (x*x<0.000001*z*z) {

		/* if ray is going down */
		if (z>0.0) {
			*s = 0.0;
			*c = 1.0;
			*t = log(v/v0)/a;
			*q = 0.5*z*(v+v0);
		
		/* else, if ray is going up */
		} else {
			*s = 0.0;
			*c = -1.0;
			*t = -log(v/v0)/a;
			*q = -0.5*z*(v+v0);
		}
	
	/* else raypath is circular with finite radius */
	} else {

		/* save translated -x to avoid roundoff error below */
		mx = -x;
		
		/* translate to make center of circular raypath at (0,0) */
		z0 = v0/a;
		z += z0;
		x0 = -(x*x+z*z-z0*z0)/(2.0*x);
		x += x0;
		
		/* signed radius of raypath; if ray going clockwise, r > 0 */
		if ( (a>0.0 && mx>0.0) || (a<0.0 && mx<0.0) )
			r = sqrt(x*x+z*z);
		else
			r = -sqrt(x*x+z*z);
		
		/* cosine at (x0,z0), cosine and sine at (x,z) */
		or = 1.0/r;
		c0 = x0*or;
		*c = x*or;
		*s = -z*or;

		/* time along raypath */
		*t = log((v*(1.0+c0))/(v0*(1.0+(*c))))/a;
		if ((*t)<0.0) *t = -(*t);

		/* integral of velocity along raypath */
		*q = a*r*mx;
	}

	/* if coordinates were rotated, unrotate cosine and sine */
	if (dvdx!=0.0) {
		temp = cr*(*s)+sr*(*c);
		*c = cr*(*c)-sr*(*s);
		*s = temp;
	}
}

void addsinc (float time, float amp,
	int nt, float dt, float ft, float *trace)
/*****************************************************************************
Add sinc wavelet to trace at specified time and with specified amplitude
******************************************************************************
Input:
time		time at which to center sinc wavelet
amp		peak amplitude of sinc wavelet
nt		number of time samples
dt		time sampling interval
ft		first time sample
trace		array[nt] containing sample values

Output:
trace		array[nt] with sinc added to sample values
*****************************************************************************/
{
	static float sinc[101][8];
	static int nsinc=101,madesinc=0;
	int jsinc;
	float frac;
	int itlo,ithi,it,jt;
	float tn,*psinc;

	/* if not made sinc coefficients, make them */
	if (!madesinc) {
		for (jsinc=1; jsinc<nsinc-1; ++jsinc) {
			frac = (float)jsinc/(float)(nsinc-1);
			mksinc(frac,8,sinc[jsinc]);
		}
		for (jsinc=0; jsinc<8; ++jsinc)
			sinc[0][jsinc] = sinc[nsinc-1][jsinc] = 0.0;
		sinc[0][3] = 1.0;
		sinc[nsinc-1][4] = 1.0;
		madesinc = 1;
	}
	tn = (time-ft)/dt;
	jt = tn;
	jsinc = (tn-jt)*(nsinc-1);
	itlo = jt-3;
	ithi = jt+4;
	if (itlo>=0 && ithi<nt) {
		psinc = sinc[jsinc];
		trace[itlo] += amp*psinc[0];
		trace[itlo+1] += amp*psinc[1];
		trace[itlo+2] += amp*psinc[2];
		trace[itlo+3] += amp*psinc[3];
		trace[itlo+4] += amp*psinc[4];
		trace[itlo+5] += amp*psinc[5];
		trace[itlo+6] += amp*psinc[6];
		trace[itlo+7] += amp*psinc[7];
	} else if (ithi>=0 && itlo<nt) {
		if (itlo<0) itlo = 0;
		if (ithi>=nt) ithi = nt-1;
		psinc = sinc[jsinc]+itlo-jt+3;
		for (it=itlo; it<=ithi; ++it)
			trace[it] += amp*(*psinc++);
	}
}

void makericker (float fpeak, float dt, Wavelet **w)
/*****************************************************************************
Make Ricker wavelet
******************************************************************************
Input:
fpeak		peak frequency of wavelet
dt		time sampling interval

Output:
w		Ricker wavelet
*****************************************************************************/
{
	int iw,lw,it,jt;
	float t,x,*wv;
	
	iw = -(1+1.0/(fpeak*dt));
	lw = 1-2*iw;
	wv = ealloc1float(lw);
	for (it=iw,jt=0,t=it*dt; jt<lw; ++it,++jt,t+=dt) {
		x = PI*fpeak*t;
		x = x*x;
		wv[jt] = exp(-x)*(1.0-2.0*x);
	}
	*w = ealloc1(1,sizeof(Wavelet));
	(*w)->lw = lw;
	(*w)->iw = iw;
	(*w)->wv = wv;
}