sureflpsvsh.c

su 的源代码库

	int nf;			/* number of frequencies in output traces */
	int *filters_phase=NULL;	/* =0 for zero phase, =1 for minimum phase fil*/
	int nfilters;		/* number of required filters */
	int wavelet_type;	/* =1 spike =2 ricker1 =3 ricker2 =4 akb */

	float dt;		/* time sampling interval */
	float tsec;		/* trace length in seconds */
	float fpeak;		/* peak frequency for output wavelet */
	float fref;		/* first frequency */
	float p2w;		/* maximum ray parameter value */
	float bp;		/* smallest ray parameter (s/km) */
	float bx;		/* beginning of range in Kms. */
	float fx;		/* final range in Kms. */
	float dx;		/* range increment in Kms. */
	float pw1,pw2,pw3,pw4;	/* window ray parameters (to apply taper) */
	float h1;		/* horizontal linear part of the source */ 
	float h2;		/* vertical linear part of the source */ 
	float m0;		/* seismic moment */
	float m1,m2,m3;		/* components of the moment tensor */

	float delta;		/* dip */
	float lambda;		/* rake */
	float phis;		/* azimuth of the fault plane */
	float phi;		/* azimuth of the receiver location */

	float sdcl,sdct;	/* standar deviation for p and s-wave vels */
	float z0=0.0;		/* reference depth */
	float zlayer;		/* thickness of random layers */
	int layer;		/* layer over on top of which to compute rand*/
	float tlag;		/* time lag in output traces */
	float red_vel;		/* erducing velocity */

	float w1=0.0;		/* low end frequency cutoff for taper */
	float w2=0.0;		/* high end frequency cutoff for taper */
	float wrefp;		/* reference frequency for p-wave velocities */
	float wrefs;		/* reference frequency for s-wave velocities */

	float epsp;		/* .... for p-wave velocities */
	float epss;		/* .... for p-wave velocities */
	float sigp;		/* .... for p-wave velocities */
	float sigs;		/* .... for s-wave velocities */
	float fs;		/* sampling parameter, usually 0.07<fs<0.12 */
	float decay;		/* decay factor to avoid wraparound */

	int *lobs;		/* layers on top of which lay the receivers */
	int *nintlayers=NULL;	/* array of number of layers to interpolate */
	int *filters_type;	/* array of 1 lo cut, 2 hi cut, 3 notch */

	float *dbpo=NULL;	/* array of filter slopes in db/octave */
	float *f1=NULL;		/* array of lo frequencies for filters */
	float *f2=NULL;		/* array of high frequencies for filters */
	float *cl;		/* array of compressional wave velocities */
	float *ql;		/* array of compressional Q values */
	float *ct;		/* array of shear wave velocities */
	float *qt;		/* array of shear Q values */
	float *rho;		/* array of densities */
	float *t;		/* array of absolute layer thickness */

	int *intlayers=NULL;	/* array of layers to interpolate */

	float *intlayth=NULL;	/* array of thicknesses over which to interp */
	float **wavefield1;	/* array for pressure wavefield component */
	float **wavefield2=NULL;/* array for radial wavefield component */
	float **wavefield3=NULL;/* array for vertical wavefield component */

	char *lobsfile="";	/* input file receiver layers */
	char *clfile="";	/* input file of p-wave velocities */
	char *qlfile="";	/* input file of compressional Q-values */
	char *ctfile="";	/* input file of s-wave velocities */
	char *qtfile="";	/* input file of shear Q-values */
	char *rhofile="";	/* input file of density values */
	char *tfile="";		/* input file of absolute layer thicknesses */
	char *intlayfile="";	/* input file of layers to interpolate */
	char *nintlayfile="";	/* input file of number of layers to interp */
	char *intlaythfile="";	/*input file of layer thickness where to inter*/
	char *filtypefile="";	/* input file of filter types to apply */
	char *fphfile="";	/* input file of filters phases */
	char *dbpofile="";	/* input file of filter slopes in db/octave */
	char *f1file="";	/* input file of lo-end frequency */
	char *f2file="";	/* input file of hi-end frequency */

	char *wfp="";		/* output file of pressure */
	char *wfr="";		/* output file of radial wavefield */
	char *wfz="";		/* output file of vertical wavefield */
	char *wft="";		/* output file of tangential wavefield */
	char *outf="";		/* output file for processing information */

	FILE *wfp_file;		/* file pointer to output pressure */
	FILE *wfr_file;		/* file pointer to output radial wavefield */
	FILE *wfz_file;		/* file pointer to output vertical wavefield */
	FILE *wft_file;		/* file pointer to output tangential wavefield*/
	FILE *outfp=NULL;	/* file pointer to processing information */
	FILE *infp;		/* file pointer to input information */

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

	/* get required parameter, seismic moment */
	if (!getparfloat("m0",&m0))	
		err("error: the seismic moment, m0, is a required parameter\n");

	/*********************************************************************/
	/* get general flags and set their defaults */
	if (!getparint("rand",&rand))			rand	= 0;
	if (!getparint("qopt",&qopt))			qopt	= 0;
	if (!getparint("stype",&stype))			stype	= 1;
	if (!getparint("wtype",&wtype))			wtype	= 1;
	if (!getparint("wfield",&wfield))		wfield	= 1;
	if (!getparint("int_type",&int_type))		int_type= 1;
	if (!getparint("flt",&flt))			flt	= 0;
	if (!getparint("vsp",&vsp))			vsp	= 0;
	if (!getparint("win",&win))			win	= 0;
	if (!getparint("wavelet_type",&wavelet_type))	wavelet_type = 1;
	if (!getparint("verbose",&verbose))		verbose	= 0;

	/* get model parameters and set their defaults */
	if (!getparint("lsource",&lsource))		lsource = 0;
	if (!getparfloat("fs",&fs)) 			fs	= 0.07;
	if (!getparfloat("decay",&decay))		decay	= 50.0;
	if (!getparfloat("tsec",&tsec))			tsec	= 2.048;

	/* get response parameters and set their defaults */
	if (!getparfloat("fref",&fref))			fref	= 1.0;
	if (!getparint("nw",&nw))			nw	= 100;
	if (!getparint("nor",&nor))			nor	= 100;
	if (!getparint("np",&np))			np	= 1300;
	if (!getparfloat("p2w",&p2w))			p2w	= 5.0;
	if (!getparfloat("bx",&bx))			bx	= 0.005;
	if (!getparfloat("bp",&bp))			bp	= 0.0;
	if (!getparfloat("fx",&fx))			fx	= 0.1;
	if (!getparfloat("dx",&dx))			dx	= 0.001;
	if (!getparfloat("pw1",&pw1))			pw1	= 0.0;
	if (!getparfloat("pw2",&pw2))			pw2	= 0.1;
	if (!getparfloat("pw3",&pw3))			pw3	= 6.7;
	if (!getparfloat("pw4",&pw4))			pw4	= 7.0;
	if (!getparfloat("h1",&h1))			h1	= 1.0;
	if (!getparfloat("h2",&h2))			h2	= 0.0;

	/* get output parameters and set their defaults */
	if (!getparint("nx",&nx))			nx	= 100;
	if (!getparfloat("dt",&dt))			dt	= 0.004;
	if (!getparint("nt",&nt))			nt	= tsec/dt;
	if (!getparint("nf",&nf))			nf	= 50;
	if (!getparfloat("red_vel",&red_vel))		red_vel	= 5;
	if (!getparfloat("fpeak",&fpeak))		fpeak	= 25.;
	if (!getparfloat("tlag",&tlag))			tlag	= 0.;

	/* get names of output files */
	if (wtype==1) {
		getparstring("wfp",&wfp);
		getparstring("wfr",&wfr);
		getparstring("wfz",&wfz);
	} else if (wtype==2) {
		getparstring("wft",&wft);
	} else err ("wtype has to be zero or one");

	/*********************************************************************/
	/* get or compute moment tensor components */
	if (stype==1) {

		/* get source parameters */
		if (!getparfloat("delta",&delta))	
			err("if stype==1, delta is a required parameter\n");
		if (!getparfloat("lambda",&lambda))	
			err("if stype==1, lambda is a required parameter\n");
		if (!getparfloat("phis",&phis))	
			err("if stype==1, phis is a required parameter\n");
		if (!getparfloat("phi",&phi))	
			err("if stype==1, phi is a required parameter\n");

		/* compute moment tensor components */
		compute_moment_tensor (wtype, phi, lambda, delta, phis, m0, 
			&m1, &m2, &m3);

	} else if (stype==2) {

		/* get moment tensor components from input */	
		if (!getparfloat("m1",&m1))	
			err("if stype==2, m1 is a required parameter\n");
		if (!getparfloat("m2",&m2))	
			err("if stype==2, m2 is a required parameter\n");
		if (!getparfloat("m3",&m3))	
			err("if stype==2, m3 is a required parameter\n");

	} else err("error, stype flag has to be one or two\n");

	/*********************************************************************/
	/* if q-option is not requesed, set corresponding parameters to zero */
	if (!getparint("layern",&layern))		layern	=0;	
	if (!getparfloat("wrefp",&wrefp))		wrefp	=0.0;
	if (!getparfloat("wrefs",&wrefs))		wrefs	=0.0;
	if (!getparfloat("epsp",&epsp))			epsp	=0.0;
	if (!getparfloat("epss",&epss))			epss	=0.0;
	if (!getparfloat("sigp",&sigp))			sigp	=0.0;
	if (!getparfloat("sigs",&sigs))			sigs	=0.0;

	/*********************************************************************/
	/* get number of layers and check input parameters */
	if (*clfile=='\0') {	/* p-wave vels input from the comand line */
		nlayers=countparval("cl");
	} else  {		/* p-wave vels input from a file */
		getparint("nlayers",&nlayers);
	}
	if (*ctfile=='\0') {	/* s-wave vels input from the comand line */
		if (nlayers !=countparval("cl")) 
			err("number of p-wave and s-wave velocities"
				"has to be the same");
	}
	if (*qlfile=='\0') { 	/* compressional q-values from comand line */
		if (nlayers !=countparval("ql")) 
			err("number of p-wave velocities and q-values"
				"has to be the same");
	}
	if (*qtfile=='\0') { 	/* shear q-values input from comand line */
		if (nlayers !=countparval("qt")) 
			err("number of p-wave velocities and shear q-values"
				"has to be the same");
	}
	if (*rhofile=='\0') { 	/* densities input from comand line */
		if (nlayers !=countparval("rho")) 
			err("number of p-wave velocities and densities"
				"has to be the same");
	}
	if (*tfile=='\0') { 	/* layer thicknesses input from comand line */
		if (nlayers !=countparval("t")) 
			err("number of p-wave velocities and thicknesses"
				"has to be the same");
	}
	if (int_type!=1 && int_type!=2) err("int_type flag has to be one or two");

	/*********************************************************************/
	/* if layer interpolation is requested, get parameters */
	if (*intlayfile !='\0') {
		getparint("nlint",&nlint);
		if ((infp=efopen(intlayfile,"r"))==NULL)
			err("cannot open file of layer interp=%s\n",intlayfile);
		intlayers=alloc1int(nlint);
		fread (intlayers,sizeof(int),nlint,infp);
		efclose(infp);
	} else if (countparval("intlayers") !=0) {
		nlint=countparval("intlayers");
		intlayers=alloc1int(nlint);
		getparint("intlayers",intlayers);
	}
	if (*nintlayfile !='\0') {
		if ((infp=efopen(nintlayfile,"r"))==NULL)
			err("cannot open file of layer inter=%s\n",nintlayfile);
		nintlayers=alloc1int(nlint);
		fread (nintlayers,sizeof(int),nlint,infp);
		efclose(infp);
	} else if (countparval("nintlayers") !=0) {
		if (nlint !=countparval("nintlayers")) 
			err("number of values in intlay and nintlay not equal");
		nintlayers=alloc1int(nlint);
		getparint("nintlayers",nintlayers);
	}
	if (*intlaythfile !='\0') {
		if ((infp=efopen(intlaythfile,"r"))==NULL)
			err("cannot open file=%s\n",intlaythfile);
		intlayth=alloc1float(nlint);
		fread (intlayth,sizeof(int),nlint,infp);
		efclose(infp);
	} else if (countparval("intlayth") !=0) {
		if (nlint !=countparval("intlayth")) 
			err("# of values in intlay and intlayth not equal");
		intlayth=alloc1float(nlint);
		getparfloat("intlayth",intlayth);
	}
	/* update total number of layers */
	if (nlint!=0) {
		for (i=0; i<nlint; i++) nlayers +=intlayers[i]-1;
	}
		
	/*********************************************************************/
	/* if random velocity layers requested, get parameters */
	if (rand==1) {
		getparint("layer",&layer);
		getparint("nrand_layers",&nrand_layers);
		getparfloat("zlayer",&zlayer);
		getparfloat("sdcl",&sdcl);
		getparfloat("sdct",&sdct);
	} else nrand_layers=0;	

	/*********************************************************************/
	/* allocate space */