susynlv.c

seismic software,very useful

Output:
r		array[nr] of reflectors
******************************************************************************
Notes:
Space for the ar, nu, xu, and zu arrays is freed by this function, since
they are only used to construct the reflectors.

This function is meant to be called only once, so it need not be very
efficient.  Once made, the reflectors are likely to be used many times, 
so the cost of making them is insignificant.
*****************************************************************************/
{
	int ir,iu,nuu,iuu,ns,is;
	float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
		*u,*s,(*xud)[4],(*zud)[4],*us;
	ReflectorSegment *rs;
	Reflector *rr;
	
	/* allocate space for reflectors */
	*r = rr = ealloc1(nr,sizeof(Reflector));

	/* loop over reflectors */
	for (ir=0; ir<nr; ++ir) {

		/* compute cubic spline coefficients for uniformly sampled u */
		u = ealloc1float(nu[ir]);
		for (iu=0; iu<nu[ir]; ++iu)
			u[iu] = iu;
		xud = (float(*)[4])ealloc1float(4*nu[ir]);
		csplin(nu[ir],u,xu[ir],xud);
		zud = (float(*)[4])ealloc1float(4*nu[ir]);
		csplin(nu[ir],u,zu[ir],zud);

		/* finely sample x(u) and z(u) and compute length s(u) */
		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=0,uu=0.0,s[0]=0.0; 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);
}

static void makeone (float v00, float dvdx, float dvdz, 
	int ls, int er, int ob, Wavelet *w,
	float xs, float zs, float xg, float zg,
	int nr, Reflector *r, int nt, float dt, float ft, float *trace)
/*****************************************************************************
Make one synthetic seismogram for linear velocity v(x,z) = v00+dvdx*x+dvdz*z
******************************************************************************
Input:
v00		velocity v at (x=0,z=0)
dvdx		derivative dv/dx of velocity v with respect to x
dvdz		derivative dv/dz of velocity v with respect to z
ls		=1 for line source amplitudes; =0 for point source
er		=1 for exploding, =0 for normal reflector amplitudes
ob		=1 to include cos obliquity factors; =0 to omit
w		wavelet to convolve with trace
xs		x coordinate of source
zs		z coordinate of source
xg		x coordinate of receiver group
zg		z coordinate of receiver group
nr		number of reflectors
r		array[nr] of reflectors
nt		number of time samples
dt		time sampling interval
ft		first time sample

Output:
trace		array[nt] containing synthetic seismogram
*****************************************************************************/
{
	int it,ir,is,ns;
	float ar,ds,xd,zd,cd,sd,vs,vg,vd,cs,ss,ts,qs,cg,sg,tg,qg,
		ci,cr,time,amp,*temp;
	ReflectorSegment *rs;
	int lhd=LHD,nhd=NHD;
	static float hd[NHD];
	static int madehd=0;
	
	/* if half-derivative filter not yet made, make it */
	if (!madehd) {
		mkhdiff(dt,lhd,hd);
		madehd = 1;
	}

	/* zero trace */
	for (it=0; it<nt; ++it)
		trace[it] = 0.0;
	
	/* velocities at source and receiver */
	vs = v00+dvdx*xs+dvdz*zs;
	vg = v00+dvdx*xg+dvdz*zg;

	/* loop over reflectors */
	for (ir=0; ir<nr; ++ir) {

		/* amplitude, number of segments, segment length */
		ar = r[ir].a;
		ns = r[ir].ns;
		ds = r[ir].ds;
		rs = r[ir].rs;
	
		/* loop over diffracting segments */
		for (is=0; is<ns; ++is) {
		
			/* diffractor midpoint, unit-normal, and length */
			xd = rs[is].x;
			zd = rs[is].z;
			cd = rs[is].c;
			sd = rs[is].s;
			
			/* velocity at diffractor */
			vd = v00+dvdx*xd+dvdz*zd;

			/* ray from shot to diffractor */
			raylv2(v00,dvdx,dvdz,xs,zs,xd,zd,&cs,&ss,&ts,&qs);

			/* ray from receiver to diffractor */
			raylv2(v00,dvdx,dvdz,xg,zg,xd,zd,&cg,&sg,&tg,&qg);

			/* cosines of incidence and reflection angles */
			if (ob) {
				ci = cd*cs+sd*ss;
				cr = cd*cg+sd*sg;
			} else {
				ci = 1.0;
				cr = 1.0;
			}

			/* if either cosine is negative, skip diffractor */
			if (ci<0.0 || cr<0.0) continue;

			/* two-way time and amplitude */
			time = ts+tg;
			if (er) {
				amp = sqrt(vg*vd/qg);
			} else {
				if (ls)
					amp = sqrt((vs*vd*vd*vg)/(qs*qg));
				else
					amp = sqrt((vs*vd*vd*vg)/
						(qs*qg*(qs+qg)));
			}
			amp *= (ci+cr)*ar*ds;
				
			/* add sinc wavelet to trace */
			addsinc(time,amp,nt,dt,ft,trace);
		}
	}
	
	/* allocate workspace */
	temp = ealloc1float(nt);
	
	/* apply half-derivative filter to trace */
	conv(nhd,-lhd,hd,nt,0,trace,nt,0,temp);

	/* convolve wavelet with trace */
	conv(w->lw,w->iw,w->wv,nt,0,temp,nt,0,trace);
	
	/* free workspace */
	free1float(temp);
}

static 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,sr,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;
	}
}

static 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++);
	}
}

static 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;
}