suea2df.c

su 的源代码库

	tr.duse = 2 ; /*u-component*/
	for (i=b1 ; i < nn-b2 ; ++i){
		++tracl; ++tracr;			
		tr.tracl = tracl; tr.tracr = tracr; tr.tracf=i;
		if(ivh==1) {  /* surface seismic */
		  tr.gx = fx+dn*k ;
		  tr.offset = tr.gx-tr.sx ;
		}
		if(ivh==2) {  /* VSP seismic */
		  tr.gelev = fz-dn*k ;
		  tr.offset = tr.sx-tr.gelev ;
		}
		 for(j=0; j<nt; ++j){tr.data[j] = ut[j][i];}
		fputtr(fp , &tr); k++;
	}
	k=0;
	tr.duse = 1 ; /*w-component*/
	for (i=b1 ; i < nn-b2 ; ++i){
		++tracl; ++tracr;			
		tr.tracl = tracl; tr.tracr = tracr; tr.tracf=i;

		if(ivh==1) {  /* surface seismic */
			tr.gx = fx+dn*k ;
			tr.offset = tr.gx-tr.sx ;
		}
		if(ivh==2) {  /* VSP seismic */
			tr.gelev = fz-dn*k ;
			tr.offset = tr.sx-tr.gelev ;
		}
		for(j=0; j<nt; ++j){
			tr.data[j] = wt[j][i];
		}

		fputtr(fp , &tr);
		k++;
	}
}


/* absorbing bc for top*/	
void abs_bc_top(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int ifx, int ilx)
{
	int i,j; float fac;
		
	for (j=0; j<wbc[0]; ++j){
		fac=r[(wbc[0]-j-1)];
		for (i=ifx; i<ilx; ++i){
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}
}

/* absorbing bc for left side*/	
void abs_bc_left(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int jfz, int jlz)
{
	int i,j;
	float fac;
		
	for (j=jfz; j<jlz; ++j){
		for (i=0; i<wbc[1]; ++i){
			fac=r[(wbc[1]-i-1)];
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}
}

/* absorbing bc for bottom*/	
void abs_bc_bot(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int ifx, int ilx)
{
	int i,j;
	float fac;
		
	for (j=nzpadded-wbc[2]; j<nzpadded; ++j){
		fac=r[(j-nzpadded+wbc[2])];
	  	for (i=ifx; i<ilx; ++i){
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}		
}

/* absorbing bc for right side*/	
void abs_bc_right(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r, 
			int *bc, int *wbc, int jfz, int jlz)
{
	int i,j;
	float fac;
		
	for (j=jfz; j<jlz; ++j){
		for (i=nxpadded-wbc[3]; i<nxpadded; ++i){
				fac=r[(i-nxpadded+wbc[3])];
				u[j][i]*=fac;
			w[j][i+1]*=fac;
				txx[j][i+1]*=fac;
			tzz[j][i+1]*=fac;
			txz[j][i]*=fac;
		}
	}
}
	
/* Free surface boundary condition where stresses are set to zero
   for the order 4 FD operater on a staggered grid */
void fs4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		txz[1][i]=0.; txz[0][i]=0.; txz[2][i]=0.;
	}
	for (i=2; i<nxpadded-1; ++i){
	 		txx[1][i]=0.; txx[0][i]=0.; 
	 		tzz[1][i]=-tzz[2][i]; tzz[0][i]=-tzz[3][i];
	}
}

/* Free surface boundary condition where velocities satisfy stress free
  surface for the order 4 FD operater on a staggered grid

** requires anisotropy symmetry parallel to free surface in 
** current implementation */

void fs4v_bc_top(float **u, float **w, float **c11,
			float **c55, int nxpadded, int nzpadded)
{
	int i;
	
	/*  	010514-split loop into two four parts as it was 
		in original FORTRAN CODE, previous versions 
		were calculating all vaulues at once	 */

	for (i=2; i<nxpadded-2; ++i){
		u[0][i]=u[2][i]
			+ 0.5*( F1*(w[2][i+1]-w[2][i]) 
			+ F2*(w[2][i+2]-w[2][i-1]) );
	}

	for (i=2; i<nxpadded-1; ++i){
		w[0][i]=w[2][i]
			+ 0.5*(1-2*(c55[2][i]/c11[2][i]))
			*( F1*(u[0][i]-u[0][i-1]) 
			+ F2*(u[0][i+1]-u[0][i-2]) );
	}
	for (i=2; i<nxpadded-2; ++i){
		u[1][i]=u[0][i]
			+ 0.5*( F1*(w[2][i+1]-w[2][i])
				+ F2*(w[2][i+2]-w[2][i-1]) );
	}
	for (i=2; i<nxpadded-1; ++i){
		w[1][i]=w[0][i]
			+ 0.5*(1-2*(c55[2][i]/c11[2][i]))
			*( F1*(u[1][i]-u[1][i-1]) 
				+ F2*(u[1][i+1]-u[1][i-2]) );
	}
}

void sym4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the TOP boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		txx[1][i]=txx[2][i];  txx[0][i]=txx[3][i]; 
	 		tzz[1][i]=tzz[2][i];  tzz[0][i]=tzz[3][i];
	 		txz[1][i]=-txz[3][i]; txz[0][i]=-txz[4][i];
	}
}

void sym4v_bc_top(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the TOP boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		u[1][i]=u[2][i];  u[0][i]=u[3][i]; 
	 		w[1][i]=-w[3][i]; w[0][i]=-w[4][i];
	}
}

void sym4s_bc_right(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the RHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		txx[j][nxpadded]=txx[j][nxpadded-4]; txx[j][nxpadded-1]=txx[j][nxpadded-3]; 
	 		tzz[j][nxpadded]=tzz[j][nxpadded-4]; tzz[j][nxpadded-1]=tzz[j][nxpadded-3];
	 		txz[j][nxpadded-1]=-txz[j][nxpadded-4]; txz[j][nxpadded-2]=-txz[j][nxpadded-3];
	}
}

void sym4v_bc_right(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the RHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		u[j][nxpadded-1]=-u[j][nxpadded-4]; u[j][nxpadded-2]=-u[j][nxpadded-3]; 
	 		w[j][nxpadded]=w[j][nxpadded-4]; w[j][nxpadded-1]=w[j][nxpadded-3];
	}
}

void sym4s_bc_left(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the LHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		txx[j][1]=txx[j][3]; txx[j][0]=txx[j][4]; 
	 		tzz[j][1]=tzz[j][3];  tzz[j][0]=tzz[j][4];
	 		txz[j][1]=-txz[j][2]; txz[j][0]=-txz[j][3];
	}
}

void sym4v_bc_left(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the LHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		u[j][1]=-u[j][2]; u[j][0]=-u[j][3]; 
	 		w[j][1]=w[j][3]; w[j][0]=w[j][4];
	}
}

void sym4s_bc_bot(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the BOTTOM boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
/*	 		txx[nzpadded-2][i]=txx[nzpadded-4][i];  txx[nzpadded-1][i]=txx[nzpadded-5][i];  */
/*	 		tzz[nzpadded-2][i]=-tzz[nzpadded-4][i]; tzz[nzpadded-1][i]=-tzz[nzpadded-5][i]; */
/*	 		txz[nzpadded-2][i]=-txz[nzpadded-3][i]; txz[nzpadded-1][i]=-txz[nzpadded-4][i];} */
	 		txx[nzpadded-2][i]=txx[nzpadded-4][i];  txx[nzpadded-1][i]=txx[nzpadded-5][i]; 
	 		tzz[nzpadded-2][i]=tzz[nzpadded-4][i]; tzz[nzpadded-1][i]=tzz[nzpadded-5][i];
	 		txz[nzpadded-2][i]=-txz[nzpadded-3][i]; txz[nzpadded-1][i]=-txz[nzpadded-4][i];
	}
}

void sym4v_bc_bot(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the BOTTOM boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
/*	 		u[nzpadded-2][i]=u[nzpadded-4][i]; u[nzpadded-1][i]=u[nzpadded-5][i];  */
/*	 		w[nzpadded-2][i]=-w[nzpadded-3][i];w[nzpadded-1][i]=-w[nzpadded-4][i];} */
	 		u[nzpadded-2][i]=u[nzpadded-4][i]; u[nzpadded-1][i]=u[nzpadded-5][i]; 
	 		w[nzpadded-2][i]=-w[nzpadded-3][i]; w[nzpadded-1][i]=-w[nzpadded-4][i];
	}
}

void calc_area(float vmax, float dtx, float dtz, 
		int *limits, int i, int j, int k, int nxpadded, int nzpadded)
{
	 limits[0] = i - NINT(vmax*k*dtx) - 10;
	 limits[1] = i + NINT(vmax*k*dtx) + 10;
	 limits[2] = j - NINT(vmax*k*dtz) - 10;
	 limits[3] = j + NINT(vmax*k*dtz) + 10;

	 if (limits[0] < 2) limits[0] = 2;

	 if (limits[1] > nxpadded-3) limits[1] = nxpadded-3;

	 if (limits[2] < 2) limits[2] = 2;

	 if (limits[3] > nzpadded-3) limits[3] = nzpadded-3;	 
}



void write_chd(int nxpadded, int nzpadded, int nt, float dx, float dz, float dt,
	float fx, float xmax, float fz, float zmax,
	float sx, float sz, float favg, float ts, char *wtype, char *stype,
	int *bc, int qsw, int aniso, float xz,
	FILE *fp, int comp, char *mfile)
{
	int k=1;
	char s[80];
	k=fwrite_chd(fp,"SUEA2DF",k); 
	k=fwrite_chd(fp,"Department of Earth Sciences",k);
	k=fwrite_chd(fp,"Uppsala University",k);
	k=fwrite_chd(fp,"Date:",k);
	k=fwrite_chd(fp," ",k);

	if(comp==1) {
		sprintf(s,"Horizontal line recording at z=%f",xz);
		k=fwrite_chd(fp,s,k);
	}

	if(comp==2) {
		sprintf(s,"Vertical line recording at x=%f",xz);
	  	k=fwrite_chd(fp,s,k);
	}
	sprintf(s,"nxpadded=%d nzpadded=%d nt=%d dx=%f dz=%f dt=%f",nxpadded,nzpadded,nt,dx,dz,dt);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"fx=%f xmax=%f fz=%f zmax=%f",fx,xmax,fz,zmax);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"sx=%f sz=%f favg=%f ts=%f",sx,sz,favg,ts);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"wtype=%s stype=%s",wtype,stype);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"bc=%d,%d,%d,%d",bc[0],bc[1],bc[2],bc[3]);
	  k=fwrite_chd(fp,s,k);
	k=fwrite_chd(fp," ",k);
	k=fwrite_chd(fp,"Model file:",k);

	sprintf(s," %s",mfile);

	k=fwrite_chd(fp,s,k);

	k=fwrite_chd(fp," ",k);

	if(qsw==1) k=fwrite_chd(fp,"Attenuation included",k);

	if(aniso==1) k=fwrite_chd(fp,"Anisotropy included",k);

	while(k<=40) k=fwrite_chd(fp," ",k);
}	
	
int fwrite_chd(FILE *fp, char *s, int k)
{
	int i,n;	

	n=fprintf(fp,"C %s",s);

	for(i=1; i<80-n-1; ++i) fputs(" ",fp); 

	fprintf(fp,"%s\n"," ");

	return(k+1);
}


void write_grid(float **u, float **w, float **txx, float **tzz,
			float **txz, int is, int js)
{
	int i,j;
		
	for (i=is-2; i<is+3; ++i){
		fprintf(stderr,"%d ",i);
	}

	fprintf(stderr,"\n");

	for (j=js-2; j<js+3; ++j){
	   	for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",u[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"w \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",w[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"txx \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",txx[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"tzz \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",tzz[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"txz \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",txz[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

}

void pad_float_2D(int n1, int n2, int *padval, float **in, float **out) 
{
/*************************************************************************
pad_float_2D - padd a 2D array of floats
**************************************************************************

Input:
n1		size of the 1 (fast) dimension
n2		size of the 2 (slow) dimension
padval[4]	paddings  	padval[0] beginning pad in dimension 1 
				padval[1] beginning pad in dimension 2
				padval[2] ending pad in dimension 1
				padval[3] ending pad in dimension 2
in[][]		n1 by n2 array of input floats
Output:
out[][]		n1padded by n2padded array of input floats
		where the padded sizes are given by
		n1padded = n1 + padval[1] + padval[3]
		n2padded = n2 + padval[0] + padval[2]
**************************************************************************
Notes:
**************************************************************************
Author: CWP: John Stockwell, 2005
**************************************************************************/

	int ix1;	/* counter in dimension 1 */
	int ix2;	/* counter in dimension 2 */
	int n1padded;	/* size of dimension 1 for padded array */
	int n2padded;	/* size of dimension 2 for padded array */

	/* calculate padded six2es */
	n1padded = n1 + padval[1] + padval[3];
	n2padded = n2 + padval[0] + padval[2];

	/* zero out the out[][] array */
	memset((void *) out[0], 0, n1padded*n2padded*FSIZE);

	/* copy in[][] into out[][]  */
	for (ix2=padval[0]; ix2 < n2; ++ix2) 
	  	for (ix1=padval[1]; ix1<n1; ++ix1) 
			out[ix2]