suea2df.c

su 的源代码库

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

/* SUEA2DF: $Revision: 1.19 $ ; $Date: 2006/11/07 22:58:42 $		 */

#include "su.h"
#include "segy.h"
#include "header.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUEA2DF - SU version of (an)elastic anisotropic 2D finite difference 	",
"		forward modeling					",
"									",
"									",
" suea2df > outfile c11file= c55file  [optional parameters]		",
"									",
" Required Parameters:							",
" c11file=c11_file	c11 voigt elasticity parameter filename		",
" c55file=c55_file	c55 voigt elasticity parameter filename		",
"									",
" Optional Parameters:							",
" rhofile=rho_file	density filename				",
"									",
" Anisotropy parameters:						",
" aniso=0	 	 =1 - include anisotropy parameters		",
" mode=0		=0 output particle velocity, =1 output stresses ",
"			(snapshots only)				",
"									",
" ... the next 3 parameters become active only when aniso=1....		",
" c13file=c13_file	c13 voigt elasticity parameter filename		",
" c33file=c33_file	c33 voigt elasticity parameter filename		",
" c15file=c15_file	c15 voigt elasticity parameter filename		",
" c35file=c35_file	c35 voigt elasticity parameter filename		",
"									",
" dt=0.001		time sampling interval (s)			",
" ft=0.0 		first time (s)				 	",
" lt=1.0 		last time (s)					",
"									",
" nx=200		number of values in slow (x-direction)		",
" dx=10.0	 	spatial sampling interval (m) x-coor		",
" fx=-1000		first x coor (m)				",
"									",
" nz=100		number of values in fast (z)-dimension		",
" dz=dx			spatial sampling interval (m) z-coor		",
" fz=0			firstz coor (m)  				",
"									",
" Source parameters:							",
" sx=0			source x position (m)				",
" sz=500		source location (m)  				",
" stype='p'		source type					",
"			  p: P-wave					",
"			  v: velocity					",
"			 pw: P plane-wave				",
" sang=0		for stype='pw': plane wave angle		",
" wtype='dg'		wave type					",
" 			 dg: Gaussian derivative 			",
" 			 ga: Gaussian		 			",
" 			 ri: Ricker					",
" 			 sp: spike, sp2: double spike   		",
" ts=0.05		source duration (s)				",
" favg=50		source average frequency			",
"									",
" Attenuation parameters:						",
" ql=		 	model q at xl,zl				",
" qsw=0		 	switch to include attenuation =1 - include	",
"									",
" Boundary condition parameters:					",
" bc=10,10,10,10 	Top,left,bottom,right boundary condition	",
" 			=0 none						",
" 			=1 symmetry 					",
" 			=2 free surface (top only)			",
" 			>2 absorbing (value indicates width of absorbing",
"			layer	 					",
" bc_a=0.95;		bc initial taper value for absorbing boundary  ",
" bc_r=0.;		bc exponential factor for absorbing boundary  	",
" 			variables are scaled by bc_a*pow(i,-bc_r)	",
"									",
" Output parameters:							",
" sofile=		name of source file				",
" snfile=		name of file containing for snapshots		",
" snaptime=		times of snapshots i.e. snaptime=0.1,0.2,0.3	",
"									",
" vsx=			x coordinate of vertical line of seismograms	",
" hsz=			z coordinate of horizontal line of seismograms	",
" vrslfile=\"vsp.su\"	output file for vertical line of seismograms[nz][nt]",
" hsfile=\"hs.su\" 	output file for horizontal line of seismograms[nx][nt]",
" tsw=0		 	switch to use shear stress only in non-fluid	",
"			media - may help reduce dispersion tsw=1. If	",
"			tsw=0 then standard calculation	  		",
" verbose=0		=1 to print progress on screen			",
"									",
" Notes:								",
" 1) The outfile contains information generated by the input parameters,",
"    such as memory allocation, stability, etc. If your input file does	",
"    not work, check this file first.					",
"									",
" The model is specified as binary files of stiffness parameters and    ",
" densities. These may be created any way the user desires. The program ",
" unif2 or makevel may be used to generate densities, and the program	",
" unif2aniso may be used to generate the stiffnesses. You will need to	",
" transpose these files, as the input format for suea2df assumes that the",
" fast dimesion is the horizontal or the x-dimension. You may do this via",
" 									",
"  transp n1=nz < c11_file > transp_c11_file				",
"									",
" If qsw=1 a ql array specifies the values of the quality factor Q, to	",
" specify the attenuation of the medium.				",
"									",
" If aniso=1 then the program will expect the additional stiffnesss files.",
"									",
" The epl and dsl are epsilon and delstar in Thomson (1986, Geophys.)	",
" anl is the angle the anisotropy makes with the horizontal.		",
" 									",
"									",
"Output files (always generated)					",
"	hsfile								",
"	vrslfile							",
"	hsfile.chd - header for hsfile					",
"	vrslfile.chd - header for vrslfile				",
"	hsfile.mod - model file						",
"									",
" Output files (if requested)						",
"	sofile - ascii source file					",
"	snfile  - su format snapshots file				",
"									",
"									",
NULL};

/* 
 * Credits: UU GEOPHY Chris Juhlin 15 May 1999
 * Copyright (c) Uppsala University, 1998.
 * All rights reserved.			
 * Parts of program use Seismic Unix Package - CSM
 * Changes - C. Juhlin
 *
 * 1. Fixed upgrading of stresses. There was an error in the coding for
 * the Tzz term, c15 was being used instead of c35. This only caused
 * problems for dipping anisotropic layers
 *
 * 2. Added some header information for hutput of snapshots.
 *
 * 3. 2001-01-30: Added option to set absorbing bc constants bc_a and bc_r 
 *
 * 4. 2001-02-23: Corrected bug in outputting model boundaries to standard
 * output in 
 * routine get_econst
 *
 * 5. 2001-04-26: Added option for updating velocities to only use 
 * shear stress if material is non-fluid, this appears to reduce dispersion at 
 * near grazing angles for fluid-solid boundary. Set tsw=1 to invoke
 *
 * 6. 2001-05-14: Changed loop in free-surface boundary condition for velocties
 * Thanks goes to Mike Holzrichter for pointing out this problem and the wrong
 * scaling factor in the updating.
 *
 * 7. 2001-05-16: Changed set_layers function to avoid negative indexing.
 * Thanks goes to Mike Holzrichter for pointing out this problem
 *
 * 8. 2001-05-17: Modified make_seis to take into account VSP geometry
 * correctly and not store unnecessary data.
 *
 * 9. 2001-08-21: Fixed set_layers so mode fills properly in depth. Earlier
 * versions were accessing incorrect array locations at last defined depth.
 *
 * 10. 2003-04-21: Fixed boundary conditions.
 *
 * 11. 2003-05-02: Extended the model area by half the grid spacing on the RHS.
 * This makes the model area symmetric allowing a plane wave source to be
 * introduced into the model (stype=pw). The w, txx and tzz grids contain now
 * one more column than the u and txz grids.
 *
 * 12. 2003-05-02: Added routines to allow plane waves to be introduced at a 
 * specified angle (sang) into the model with functions add_pw_source_V and
 * add_pw_source_S.
 *
 * 13. 15 Oct 2005 -- tossed all the model building stuff. Read models
 *      from binary files made by  unif2aniso
 *
 * Alogrithm based on Juhlin (1995, Geophys. Prosp.)
 *	and Levander (1988, Geophysics)
 * Attenuation included as in Graves (1996, BSSA)
 *
 */

/**************** end self doc ********************************/

#define	F1	1.125
#define	F2	-0.04166667

segy tr, trv, trh;

/* prototypes for functions defined and used below */
int get_source(float dt, float ts, float favg, char *wtype, float *source);
void update_vel(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz);
void update_vel_tsw(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz);
void update_stress_iso(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float dtx, float dtz);
void update_stress_ani(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float **c33, float **c13, 
	float **c15, float **c35, float dtx, float dtz);
void add_p_source(float **txx, float** tzz, float amp, int i, int j);

void add_v_source(float **u, float **w, float amp, int i, int j);

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);
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);
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);
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);
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);
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);
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);
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);
void fs4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void fs4v_bc_top(float **u, float **w, float **c11, float **c55, 
			int nxpadded, int nzpadded);
void sym4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_top(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_left(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_left(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_bot(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_bot(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_right(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_right(float **u, float **w, int nxpadded, int nzpadded);

void calc_area(float vmax, float dtx, float dtz, int *limits, int i,
			int j, int k, int nxpadded, int nzpadded);
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 fwrite_chd(FILE *fp, char *s, int k);
void write_grid(float **u, float **w, float **txx, float **tzz, 
			float **txz, int is, int js);
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);
void pad_float_2D(int n1, int n2, int *padval, float **in , float **out); 
		
/* the main program */
int
main (int argc, char **argv)
{
	float dt;		/* time sampling interval 	*/
	float ft;		/* first time sample		*/
	float lt;		/* last time sample		*/

	float dx;		/* spatial sampling interval x 	*/
	float dz;		/* spatial sampling interval z	*/
	float fx;		/* min x coor		 	*/

	float xmax;		/* max x coor		 	*/
	float fz;		/* min z coor		 	*/
	float zmax;		/* max z coor		 	*/

	char *stype=NULL;	/* source type			*/
	char *wtype=NULL;	/* wave type			*/
	float sang;		/* plane wave angle (deg)	*/

	float ts;		/* source length (s)		*/
	float favg;		/* source average frequency (Hz)*/
	float sx,sz;		/* current source location (m)	*/

	int nsnap;		/* number of snapshots		*/
	float *snaptime=NULL;	/* snapshot times		*/
	int *isnap=NULL;	/* snapshot samples		*/

	int nt=0;		/* number of time samples 	*/
	int nz=0;		/* number of z samples of model */
	int nx=0;		/* number of x samples of model	*/
	int nzpadded=0;		/* nz padded by bc		*/
	int nxpadded=0;		/* nx padded by bc		*/

	int nbytes;		/* number of bytes allocated	*/
	int ngrids;		/* number of bytes allocated	*/
	int i,j,k,ix,iz;	/* loop counters 		*/

	float **c11=NULL;	/* elastic constant c11(x,z)	*/
	float **c11temp=NULL;	/* ...related temporary array	*/
	float **c55=NULL;	/* elastic constant c55(x,z)	*/
	float **c55temp=NULL;	/* ...related temporary array	*/
	float **rho=NULL;	/* density rho(x,z), initially	*/
				/* ... later rho(x,z)=1/density	*/
	float **rhotemp=NULL;	/* ...related temporary array	*/
	float **q=NULL;		/* Q-factor q(x,z)		*/

	float rhomin=0.0;	/* minimum density		*/
	float rhomax=0.0;	/* maximum density		*/
	float c11min=0.0;	/* minimum c11			*/
	float c11max=0.0;	/* maximum c11			*/
	float c55min=0.0;	/* minimum c55			*/
	float c55max=0.0;	/* maximum c55			*/

	float c13min=0.0;	/* minimum c13			*/
	float c13max=0.0;	/* maximum c13			*/
	float c33min=0.0;	/* minimum c33			*/
	float c33max=0.0;	/* maximum c33			*/
	float c15min=0.0;	/* minimum c15			*/
	float c15max=0.0;	/* maximum c15			*/
	float c35min=0.0;	/* minimum c35			*/
	float c35max=0.0;	/* maximum c35			*/

	int qsw;		/* switch to include attenuation*/

	int aniso;		/* switch to include anisotropy */
	float **c13=NULL;	/* elastic constant c13(x,z)	*/
	float **c13temp=NULL;	/* ...associated temp file	*/
	float **c33=NULL;	/* elastic constant c33(x,z)	*/
	float **c33temp=NULL;	/* ...associated temp file	*/
	float **c15=NULL;	/* elastic constant c13(x,z)	*/
	float **c15temp=NULL;	/* ...associated temp file	*/
	float **c35=NULL;	/* elastic constant c35(x,z)	*/
	float **c35temp=NULL;	/* ...associated temp file	*/

	float **u=NULL,**w=NULL;	/* particle velocities		*/
	float **txx=NULL,**tzz=NULL;	/* normal stresses		*/
	float **txz=NULL;		/* shear stress			*/
	float **ut=NULL,**wt=NULL;	/* time sections 		*/

	int tsw;		/* switch for shear stress updating */

	float *source=NULL;	/* source waveform		*/
	float pfac;		/* P-wave source scaling factor	*/
	int ns;			/* number of samples in source	*/

	float energy,efac;	/* energy in model		*/
	float dtx,dtz;		/* fd factors			*/
	int ifx,ilx,jfz,jlz;	/* grid area for operators	*/
	int is,js;		/* source location in grid	*/
	float vlim[2];		/* min and max velocity in grid	*/
	int *limits=NULL;	/* calculation area in grid	*/

	int verbose;		/* is verbose? 			*/

	int vs2;		/* depth in samples of horiz rec line */
	int hs1;		/* horiz sample of vert rec line */
	float hsz;		/* z position of horiz receiver line */
	float vsx;		/* x position of vertical receiver line */

	int bc[4];		/* boundary condition flags 	*/
	int wbc[4];		/* boundary width	 	*/
	int nbc;		/* number of values given 	*/
	float *bc0=NULL;	/* taper at top boundary 	*/
	float *bc1=NULL;	/* taper at left boundary 	*/
	float *bc2=NULL;	/* taper at bottom boundary 	*/
	float *bc3=NULL;	/* taper at right boundary 	*/
	float bc_a;		/* bc initial taper value 	*/
	float bc_r;		/* bc exponential factor 	*/
	float *vl=NULL,*vr=NULL,*vb=NULL;/* velocity at boundaries 	*/
	float *ttl=NULL,*ttr=NULL,*ttb=NULL;/* traveltimes at boundaries */

	float t;		/* current time of propagation	*/


	char *sofile="";	/* source filename		*/
	FILE *soeisfp=NULL;	/* ... its file pointer  */

	char *snfile="";	/* snapshot filename		*/
	FILE *sneisfp=NULL;	/* ... its file pointer  */

	char *mfile=NULL;	/* model filename*/
	FILE *mfp=NULL;		/* ... its file pointer   */
	char *vrslfile="";	/* vert receiver seismogram line file  name */
	FILE *vseisfp=NULL;	/* ... its file pointer  */
	char *vrshfile=NULL;	/* vert receiver seismogram header filename*/
	FILE *vchdfp=NULL;	/* ... its file pointer  */

	char *hsfile="";	/* horiz receiver seismogram line filename*/
	FILE *hseisfp=NULL;	/* ... its file pointer  */
	char *hcfile=NULL;	/* horiz receiver seismogram header filename*/
	FILE *hchdfp=NULL;	/* ... its file pointer   */
	char tfile[256];	/* temporary filename*/

	char *rhofile="";	/* density filename*/
	FILE *rhofilefp=NULL;	/* ... its file pointer */
	char *c11file="";	/* c11 stiffness filename */
	FILE *c11filefp=NULL;	/* ... its file pointer */
	char *c55file="";	/* c55 stiffness filename */
	FILE *c55filefp=NULL;	/* ... its file pointer */
	char *c13file="";	/* c13 stiffness filename */
	FILE *c13filefp=NULL;	/* ... its file pointer */
	char *c15file="";	/* c15 stiffness filename */
	FILE *c15filefp=NULL;	/* ... its file pointer */
	char *c33file="";	/* c33 stiffness filename */
	FILE *c33filefp=NULL;	/* ... its file pointer */
	char *c35file="";	/* c35 stiffness filename */
	FILE *c35filefp=NULL;	/* ... its file pointer */

	float sa,ca,p;

	int mode;		/* output mode =0 for particle velocity */
				/* =1 for stresses */

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

	/* get optional parameters */
	if (!getparfloat("dt",&dt)) dt = 0.001;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("lt",&lt)) lt = 1.0;

	/* get model size in slow direction */
	if (!getparint("nx",&nx)) nx = 200;
	if (!getparfloat("dx",&dx)) dx = 10;
	if (!getparfloat("fx",&fx)) fx = -1000;

	/* get model size in fast direction */
	if (!getparint("nz",&nz)) nz = 100;
	if (!getparfloat("dz",&dz)) dz = dx;
	if (!getparfloat("fz",&fz)) fz = 0;

	/* calculate maximum model coordinates */
	xmax = fx + dx*nx;
	zmax = fz + dz*nz;

	/* get source location */
	if (!getparfloat("sx",&sx)) sx = (fx+xmax)/2;
	if (!getparfloat("sz",&sz)) sz = (fz+zmax)/2;