sureflpsvsh.c

su 的源代码库

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

#include "su.h"
#include "segy.h"
#include "header.h"
#include "Reflect/reflpsvsh.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUREFLPSVSH - REFlectivity modeling of either PSV or SH waves for layered",
"		earth model						",
"									",
" surelfpsvsh  required parameters [optional parameters]	 		",
"									",
" Required Parameters:							",
" m0=			Seismic moment					",
" p2w=			maximum ray parameter value to which the	",
"			synthetics are computed				",
"									",
" Optional Parameters:							",
" Main Flags:								",
" wtype=1		=1 for PSV, =2 for SH				",
" stype=1		=1 if the moemnt tensor components are given	",
"			=2 if they are to be computed from fault	",
"			 plane mechanism parameters			",
" wfield=2		=1 for displacement, =2 for particle velocity,	",
"			=3 for particle acceleration			",
" flt=0			=1 to apply earth flattening correction		",
"			(important for earthquake data)			",
" vsp=0			=0 for surface data, =1 for VSP data		",
" int_type=1		=1 to compute the slowness integration using	",
"			the trapezoidal rule. =2 to use a first order	",
"			Filon scheme (faster but maybe noisier)		",
" verbose=0		=0 no processing information is output (makes	",
"			the program run a litle faster), =1 to output	",
"			processing information to the screen, =2 to	",
"			output processing information to a user		",
"			supplied file, =3 output processing		",
"			information to both, the screen and a file	",
"									",
" Flags for special non-standard options (used only in rare ocasions)	",
" rand=0		=1 to include random velocity and q layers	",
" qopt=0		=1 if a q-correction is desired			",
"									",
" Flags for output data							",
" win=0			=1 to apply a frequency domain Hanning		", 
"			window prior to inverse FFT to time domain	",
" wavelet_type		=1 for a spike, =2 for a Ricker wavelet		",
"			=3 for an akb wavelet				",
"									",
" Main input parameters							",
" tsec=2.048		length of computed traces (in seconds)		",
" dt=0.004		time sampling interval (in seconds)		",
" nt=tsec/dt		number of samples per trace			",
" nx=60			number of traces (ranges) per shot		",
" nw=100		number of frequencies to process		",
" nor=1			number of depths at which receivers are		",
"			located (different from one for VSP's)		",
" nlayers=10		number of horizontal layers in the model	",
" fref=1.0		reference frequency (Hz)			",
" bx=0.0 		first range (first trace offset) (km)		",
" dx=0.05		range increment (trace to trace distance) km	",
" pw1=0.0 pw2=0.1	apply Hanning window tapering to lower end	",
"			of ray parameter computations (if set to zero	",
"			good default values are used)			",
" pw3=6.7 pw4=7.0	apply Hanning window tapering to higher end	",
"			of ray parameter computations (if set to zero	",
"			good default values are used)			",
" fs=0.07		Filon sampling parameter, usually between	",
"			0.07 and 0.12. Sampling is finer as fs		",
"			increases. For short range synthetics		",
"			(<100km) use 0.07, for medium range (<1000km)	",
"			use 0.09 and for large	ranges (>1000km) use	",
"			0.12 (this parameter is ignored if int_type=1)	",
" np=1300		number of ray parameters used to compute the	",
"			seismograms. Automatically set if int_type=2	",
" bp=0.0		slowest ray parameter used to compute the	",
"			seismograms. Set to zero if int_type=2		",
" decay=50.0		decay factor use to avoid time series		",
"			wraparound. A value 'n' for decay means that	",
"			the wrapped around signals are diminished by	",
"			a factor of 'n'. The default of 50 is the	",
"			recomended value				",
" lobs=			array[nor] of layers on top of which the	",
"			receivers are located. If any receiver is	",
"			within a layer, introduce a ficticious layer	",
" cl=			array[nlayers] of compressional velocities	",
"			for each layer in km/s				",
" ct=			array[nlayers] of shear velocities for each	",
"			layer in km/s. A value of -1 assumes		",
"			ct=cl/sqrt(3.					", 
" ql=			array[nlayers] of compressional q values	",
"			for each layer					",
" qt=			array[nlayers] of shear q values for each	",
"			layer. A value of -1 assumes qt=ql/2.		",
" rho=			array[nlayers] of densities for each layer	",
"			in gr/cc					",
" t=			absolute thickness of the layer in kms		",
"									",
" Note: any of these arrays can also be input in a file as lobsfile=,clfile=,",
" ctfile=, etc.								",
"									",
" Source parameters:							", 
" lsource=1		layer on top of which the source is located	",
" h1=1.0 h2=0.0		vertical and horizontal linear part of the	",
"			source. For unit point source set h1=0, h2=0	",
"  			and all other source parameters to zero. For	",
"			earthquake source both h1 and h2 are usually	",
"			zero					 	",
" m1,m2,m3		[1][1],[1][2] and [2][2] components of the	",
"			moment tensor ([1][2]=[2][1] component)		",	
" delta=0.0		dip in degrees (necessary only if stype=1)	",
" lambda=0.0		rake in degrees (necessary only if stype=1)	",
" phis=0.0		fault plane azimuth in degrees (necessary	",
"			only if stype=1)				",
" phi=0.0		azimuth of the receiver location in degrees	",
"			(necessary only if stype=1)		  	",
" Parameters for output data:						",
" tlag=0.0		time lag to appy to the  seismograms (sec)	",
" nf=nw			number of frequencies in output data		",
" fpeak=25.0		peak frequency for Ricker or akb wavelet (hz)	",
" red_vel=0.0		reducing velocity (km/s). If set to zero, the	",
"			maximum compressional velocity is used		",
"			(this parameter is garanteed not to be <1.5)	",
" w1=0.0		A Hanning window will be applied to freqs	",
"			less than w1.If set to zero, the default of	",
"			0.15 of maximum frequency is used (provided	",
"			that win=1)					",
" w2=0.0		A Hanning window will be applied to freqs	",
"			greater than w2. If set to zero, the default	",
" 			0.85 of maximum frequency is used (provided	",
"			that win=1)					",
" nfilters=0		number of filters to apply to the synthetics	",
" filters_phase=	array[nfilters] of: 0 for zero phase filters	", 
" 			or 1 for minimum phase filters(can also be	",
"			input in a file via fphfile=			",
" filters_type=		array[nfilters] of: 1 for high cut filters,	",
"			input in a file via filphasefile=		",
" 			2 for low cut filters, 3 for notch filters	",
"			(can also be input in a file via filtypefile=)	",
" dbpo=			array[nfilters] of: filter slopes in db/oct	",
"			(can also be input in a file via dbpofile=)	",
" f1=			array[nfilters] of: frequency to start filter	",
"			action (Hz)(can also be input in a file via	",
"			f1file=)					",
" f2=			array[nfilters] of: frequency to end filter	",
"			action (Hz). Only for notch filters		",
"			(can also be input in a file via f2file=)	",
" wfp=			name of output pressure seismogram file		",
"			(only if wtype=1)				",
" wfr=			name of output radial comp seismogram file	",
"			(only if wtype=1)				",
" wfz=			name of output vertical comp seismogram file	",
"			(only if wtype=1)				",
" wft=			name of output tangential comp seismogram file	",
"			(only if wtype=1)			  	",
" outf=info		name of output processing information file	",
"			(only if verbos=2 or 3)				",
"									",
" Interpolation parameters (required only if layers with gradients desired)",
" nlint=0		number of times layer interp is required	",
" nintlayers=		array[nlint] of number of layers to interpol	",
"			each time (can be input as file nintlayfile=)	", 
" intlayers=		array[nlint] of layers on top of which to	",
"			start each interpolation			",
"			(can be input as file intlayfile=)		",
" intlayth=		array[nlint] of layer thicknesses to interp	",
"			(can be input as file intlaythfile=)		",
"									",
" Other parameters (required only under very special circumstances)	",
" nrand_layers=0	maximum number of random layers allowed		",
"			(only if rand=1)				",
" layer=0		layer on top of which the random velocity	",
"			layers are inserted (only if rand=1)		",
" zlayer=0.0		thickness of random layers (only if rand=1)	",
"   			(if zlayer<t[il], then zlayer=t[il])		",
" sdcl=0.0		standard deviation for compressional vels	",
"			(only if rand=1)				",
" sdct=0.0		standard deviation for shear velocities		",
"			(only if rand=1)				",
" layern=0		layer on top of which the q-option is invoked	",
"			(only if qopt=1)				",
" wrefp=1.0		reference frequency for compressional vels	",
"			(only if qopt=1)				",
" wrefs=1.0		reference frequency for shear velocities	",
"			(only if qopt=1)				",
" epsp=0.001		reference amplitude for comporessional vels	",
"			(only if qopt=1)				",
" epss=0.001		reference amplitude for shear velocities	",
"			(only if qopt=1)				",
" sigp=0.1		xxxxxx for comporessional vels			",
"			(only if qopt=1)				",
" sigs=0.1		xxxxxx for shear velocities		  	",
"			(only if qopt=1)				",
" Notes:								",
" Gradient zones between two layers can be handled with the use of the	",
" layer interpolation option. The program will automatically compute and",
" insert the required number of layers with the appropriate thicknesses	",
" between the layers above and below. There is no restriction as to how	",
" many layers the program can handle, provided enough computer power is	",
" available.								",
"									",
" The number of frequencies to be processed is the most critical parameter",
" to determine how long the program will take to run. The maximum frequency",
" that can be present in the seismogram is nw/tsec Hz, however, a lower	",
" frequency can be selected for the output via the nf param.		",
"									",
" The number of computed time samples is tsec/dt, however, the user may	",
" choose a smaller number of samples for the output traces by setting	",
" nt this can be useful when a short seismogram is required with a	",
" broadband frequency range (directly computing a small seismogram can	",
" be hazardous)								",
"									",
" The decay parameter should be chosen with care. A value of 50 seems	",
" to give good results, but depending on the data this parameter can	",
" boost up late wrapped around noncausal energy.			", 
"									",
" The integration flag is important for short range data, in particular ",
" for oil exploration, for which the first order Filon scheme, though	",
" faster, can produce noisier seismograms, specialy in the late parts of",
" the record. A standard trapezoidal rule seems to work better but is	",
" slower by about 25% in a normal situation.				",
"									",
" When arrays are required as input, they can also be input as files,	",
" however, the program will only check that the number of parameters is ",
" the same if the arrays are used or if the number of elements in the	",
" files are set. A combination of arrays and files is also permited.	",
"									",
" Examples of Source Parameters:					",
" For a vertical point force, set stype=1, h1=1, h2=0 and all moment	",
" tensor components to zero. Ignore all other source parameters.	",
" For an explosion, set stype=1, h1=h2=0, m1=m3=A, m2=0 where A is some	",
" constant (normally one). If A is negative, an implosion is generated	",
" instead.  Ignore all other source parameters.				",
" For a fault slip, set stype=2, set phis,phi,lambda,alpha and m0 to their",
" corresponding values and ignore the moment tensor components		",
"									",
NULL};

/*
 * Credits:
 *
 * Original Fortran 77 PSV version written by Subshashis Mallick (1988)
 * Original Fortran 77 SH version written by Mrinal Sen (1988) based on the 
 * PSV version by Subshashis Mallick 
 * Translated to C, expanded and reformatted for SU by Gabriel Alvarez (1995)
 *
 * References:
 * The reflectivity method: a tutorial. G Muller. J. Geophysics(1985)
 *	v. 58.  153-174
 * Practical aspects of reflectivity modeling. Mallick and Frazer. Geophysics
 *	v. 52 No. 10. October 1987. 1355-1364
 *
 */
/************************** end self doc *************************************/

#define RSO 6371.0

segy tr1,tr2,tr3;


int main(int argc, char **argv)
{
	int i,ix,it;		/* loop counters */
	int wtype;		/* =1 psv. =2 sh wavefields */
	int wfield;		/* =1 displcement =2 velocity =3 acceleration */
	int stype;		/* source type */
	int int_type;		/* =1 for trapezoidal rule. =2 for Filon */
	int flt;		/* =1 apply earth flattening correction */
	int rand;		/* =1 for random velocity layers */
	int qopt;		/* some flag ???? */
	int vsp;		/* =1 for vsp, =0 otherwise */
	int win;		/* =1 if frequency windowing required */
	int verbose;		/* flag to output processing information */
	int nt;			/* samples per trace in output traces */
	int ntc;		/* samples per trace in computed traces */
	int nx;			/* number of output traces */
	int np;			/* number of ray parameters */
	int nlint=0;		/* number of times layer interp is required */
	int lsource;		/* layer on top of which the source is located*/
	int nw;			/* number of frequencies */
	int nor;		/* number of receivers */
	int nlayers;		/* number of reflecting layers */
	int layern;
	int nrand_layers;	/* maximum number of random layers permitted */