suea2df.c

su 的源代码库

		 		+ F2*(w[j][i+2]-w[j][i-1]);

		duz = F1*(u[j][i]-u[j-1][i])
	 				+ F2*(u[j+1][i]-u[j-2][i]);

	 			dwz = F1*(w[j+1][i]-w[j][i])
			+ F2*(w[j+2][i]-w[j-1][i]);

	 			dux = F1*(u[j][i]-u[j][i-1])
			+ F2*(u[j][i+1]-u[j][i-2]);

		txx[j][i] = txx[j][i] 
			+ c11[j][i]*dtx*dux 
			+ c13[j][i]*dtz*dwz 
			+ c15[j][i]*(dtz*duz + dtx*dwx);

	 			tzz[j][i] = tzz[j][i] 
			+ c33[j][i]*dtz*dwz 
			+ c13[j][i]*dtx*dux 
			+ c35[j][i]*(dtz*duz + dtx*dwx); 
	}
}


void add_p_source(float **txx, float **tzz, float amp, int i, int j)
{
	txx[j][i] = txx[j][i] + amp;
	tzz[j][i] = tzz[j][i] + amp;
}

void add_v_source(float **u, float **w, float amp, int i, int j)
{
	w[j][i] = w[j][i] + amp*0.000001;
}


void add_pw_source_V(float **u, float **w, float **txx, float **tzz, 
	float **txz, float **c11, float **c55, float *a, float sang,
	int il, int ir, int jt, int jb, int ns, int k, float dx,
	float dz, float dt, float *vl, float *vb, float *vr,
	float *ttl, float *ttb, float *ttr)
{
	int i,j;
	float sa, ca, ak, p;
	float tw=0, tu=0, tiw, tiu, tt;
	int ksw, ksu;
	double pi=D_PI;

	/* sang = direction of wave propagation relative to vertical */
	tt=k*dt;

	/* introduce plane wave on left boundary */
	p=sin(sang*pi/180.)/vl[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));

		for (i=il-2; i<il; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttl[j]-(dx*(i-il+0.75)*sa 
					+ dz*(-0.25)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.75)*ca)/vl[j];
	   		}
	   		if(sang>-90&&sang<0) { 
					tw=tt-ttl[j]-(dx*(i-ir-1.75)*sa 
					+ dz*(-0.25)*ca)/vl[j];
					tu=tt-ttl[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.75)*ca)/vl[j];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
	   		ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;
				if (ksw>=0 && ksw <ns-1) {
					ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i]=-vl[j]*ak*ca;
			} 

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vl[j]*ak*sa;
			}
		}
	}

	/* introduce plane wave on right boundary */
	p=sin(sang*pi/180.)/vr[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));

		for (i=ir+1; i<ir+3; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttr[j]-(dx*(i-il+1.75)*sa 
					+ dz*(-0.25)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.75)*ca)/vr[j];
			}
			if(sang>-90&&sang<0) { 
				tw=tt-ttr[j]-(dx*(i-ir-0.75)*sa 
					+ dz*(-0.25)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.75)*ca)/vr[j];
		   	}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i+1]=-vr[j]*ak*ca;
			}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vr[j]*ak*sa;
			}
		}
	}
	 
	/* introduce plane wave on bottom boundary */
	 if(sang>=0&&sang<90) p=sin(sang*pi/180.)/vb[il];
	 if(sang>-90&&sang<0) p=sin(sang*pi/180.)/vb[ir];
	 for (i=il; i<ir+2; ++i){
		sa=sin(asin(p*vb[i]));
		ca=cos(asin(p*vb[i]));

	  	for (j=jb-1; j<jb+3; ++j) {
			if(sang>=0&&sang<90) { 
					tw=tt-ttb[i]-(dx*(+0.75)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
					tu=tt-ttb[i]-(dx*(+1.25)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
	   		}
	   		if(sang>-90&&sang<0) { 
					tw=tt-ttb[i]-(dx*(-1.75)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
					tu=tt-ttb[i]-(dx*(-1.25)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
					ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i]=-vb[i]*ak*ca;
	   		}
			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vb[i]*ak*sa;
	   		}
		}	
	}

}

void add_pw_source_S(float **u, float **w, float **txx, float **tzz,
	float **txz, float **c11, float **c55, float *a, float sang,
	int il, int ir, int jt, int jb, int ns, int k, float dx,
	float dz, float dt, float *vl, float *vb, float *vr,
	float *ttl, float *ttb, float *ttr)
{
	int i,j;
	float sa, ca, ak, p;
	float tw=0, tu=0, tiw, tiu, tt;
	int ksw, ksu;
	double pi=D_PI;

	/* sang = direction of wave propagation relative to vertical */
	tt=k*dt;

	/* introduce plane wave on left boundary */
	p=sin(sang*pi/180.)/vl[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));

		for (i=il-2; i<il; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttl[j]-(dx*(i-il+0.75)*sa 
					+ dz*(-0.75)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.25)*ca)/vl[j];
	   		}
			if(sang>-90&&sang<0) { 
				tw=tt-ttl[j]-(dx*(i-ir-1.75)*sa 
					+ dz*(-0.75)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.25)*ca)/vl[j];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				txx[j][i]=(c11[j][i]*sa*sa 
					   +(c11[j][i]-2*c55[j][i])*ca*ca)*ak;
				tzz[j][i]=(c11[j][i]*ca*ca
					   +(c11[j][i]-2*c55[j][i])*sa*sa)*ak;
	   		}
			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
	   		}
		}	
	}

	/* introduce plane wave on right boundary */
	 p=sin(sang*pi/180.)/vr[jb];

	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));
		for (i=ir+1; i<ir+3; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttr[j]-(dx*(i-il+1.75)*sa 
					+ dz*(-0.75)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.25)*ca)/vr[j];
	   		}
			if(sang>-90&&sang<0) { 
				tw=tt-ttr[j]-(dx*(i-ir-0.75)*sa 
					+ dz*(-0.75)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.25)*ca)/vr[j];
	   		}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;

			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
					txx[j][i+1]=(c11[j][i]*sa*sa
							+(c11[j][i]
							-2*c55[j][i])*ca*ca)*ak;
					tzz[j][i+1]=(c11[j][i]*ca*ca
							+(c11[j][i]
							-2*c55[j][i])*sa*sa)*ak;
			} else  {
				txx[j][i+1]=0;
				tzz[j][i+1]=0;
			}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
			} else {
				txz[j][i]=0;
			}
		}
	}
	 
	/* introduce plane wave on bottom boundary */
	p=sin(sang*pi/180.)/vb[il];
	for (i=il; i<ir+2; ++i){
		sa=sin(asin(p*vb[i])); ca=cos(asin(p*vb[i]));
		for (j=jb-1; j<jb+3; ++j) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttb[i]-(dx*(+0.75)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
				tu=tt-ttb[i]-(dx*(+1.25)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
			}
			if(sang>-90&&sang<0) { 
				tw=tt-ttb[i]-(dx*(-1.75)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
				tu=tt-ttb[i]-(dx*(-1.25)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
	   		}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;
			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				txx[j][i]=(c11[j][i]*sa*sa
						+(c11[j][i]
							-2*c55[j][i])*ca*ca)*ak;
				tzz[j][i]=(c11[j][i]*ca*ca
						+(c11[j][i]
							-2*c55[j][i])*sa*sa)*ak;
	   		}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
	   		}
		}	
	}

}

void make_snap(float **u, float **w, float sx, float sz, 
		float dx, float dz, int nxpadded, int nzpadded,
		int nx, int nz, int k, float t, float fx, float fz,
		segy tr, FILE *sneisfp, int *wbc)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr;		/* trace number within a reel */
	int i,j;
	
	tr.sx = sx;
	tr.selev = -sz;
	tr.trid = 1;
	tr.fldr = k;
	tr.ns = nz;
	tr.dt = 1000 * dz;
	tr.d2 = dx ;
	tr.f1 = fz ;
	tr.f2 = fx ;

	if(&tracl==NULL) tracl = 0;
	tracr = 0;
	tr.trid = TINLIN; /* inline particle velocity component*/
	for (i=0 ; i < nxpadded-wbc[3]-wbc[1] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nzpadded-wbc[2]-wbc[0]; ++j){
			tr.data[j] = u[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.trid = TVERT ; /* vertical particle velocity component */
	for (i=0 ; i < nxpadded-wbc[3]-wbc[1] ; ++i){
		++tracl;
		++tracr;

		tr.gx = fx+dx*(i-wbc[1]);
		tr.offset = tr.gx-tr.sx ;
		tr.tracl = tracl;
		tr.tracr = tracr;

		for(j=0; j<nzpadded-wbc[2]-wbc[0]; ++j){
			tr.data[j] = w[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
}

void make_stress_snap(float **txx, float **tzz, float **txz, 
	float sx, float sz, 
	float dx, float dz, int nxpadded, int nzpadded, int k, float t,
	float fx, float fz,
	segy tr, FILE *sneisfp, int *wbc)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr=0;		/* trace number within a reel */
	int i,j;
	
	tr.sx = sx;
	tr.selev = -sz;
	tr.trid = 1;
	tr.fldr = k;
	tr.ns = nzpadded ;
	tr.dt = 1000 * dz;
	tr.d2 = dx ;
	tr.f1 = fz ;
	tr.f2 = fx ;
	tr.delrt = t ;	/*store time of snapshot*/

	if(&tracl==NULL) tracl = 0;
	tracr = 0;
	tr.trid = TINLIN ; /*txx-component in line component*/
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = txx[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.trid = TVERT ; /*tzz-component*/
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.offset = tr.gx-tr.sx ;
		tr.tracl = tracl; tr.tracr = tracr;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = tzz[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.duse = TXLIN ; /* txz-component using crossline label */
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = txz[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
}


void make_seis(float **ut, float **wt, float sx, float sz, 
	float dn, float dt, int nn, int nt, float dd,
	 float fx, float fz,
	segy tr, FILE *fp, int b1, int b2, int ivh, int verbose)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr=0;		/* trace number within a reel */
	int i,j,k;
	
	if (verbose) warn("%d %d %d %d",tr.ns,tr.tracl,tr.tracr,tr.ntr);

	tr.sdepth = 0;
	tr.sx = sx;
	tr.selev = -sz ;
	tr.ns = nt ;
	tr.d2 = dn ;
	tr.dt = dt * 1000000.0 ;
	tr.ntr = 2*(nn-b2-b1);

	if(ivh==1) {  /* surface seismic */
		tr.gelev = -dd ;
		tr.f2 = fx ;
	}
	if(ivh==2) {  /* VSP seismic */
	  tr.gx = dd ;
	  tr.f2 = fz ;
	}
	
	tracl = tracr = 0; k=0;