auxil.c

su 的源代码库

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


#include "par.h"
#include "Reflect/reflpsvsh.h"
#define BLK 100

/*********************** self documentation **********************/
/******************************************************************************
AUX_REFLEXTIVITY - Set of subroutines to compute auxiliary stuff to be
	used by the reflectivity modeling code in producing synthetic seismograms
*******************************************************************************
Function prototypes:
void compute_w_aux_arrays (int wtype, int layern, int nlayers, int nor, 
	int lsource, int *np, int *block_size, int *nblock, int *left, float fw,
	float wrefp, float wrefs, int *lobs, float tsec, float p2w, float fs, 
	float xmax, float w, float decay, float epsp, float epss, float sigp, 
	float sigs, float *pw1, float *pw2, float *pw3, float *pw4, float *dp, 
	float *fp, complex wpie, complex *cdp, float *cl, float *ct, float*ql, 
	float *qt, float *rho, float *t, complex *al, complex *at, complex *prs);

void compute_p_aux_arrays (int wtype, int nlayers, int lsource, int block_size,
	float bp, float dp, float w, float decay, float pw1, float pw2,
	float pw3, float pw4, complex *pwin, float m1, float m2, float m3, float h1,
	float h2, complex wpie, complex *divfac, complex *p, complex *pp,
	complex *al, complex *at, float *rho, complex *gl, complex *gt,
	complex *gam, complex *alpha, complex *betha, complex *sigmad1,
	complex *sigmad2, complex *sigmau1, complex *sigmau2) ;

void source_receiver_type (int nor, int nlayers, int lsource, int  *lobs, 
	float *cl, float *ct, int *acoust, int *flag);

void compute_slowness (int wtype, int lsource, float w, float p2w, int *np, 
	float *fp, float fs, float e4, float dk, float decay, float xmax, float *dp,
	float *pw1, float *pw2, float *pw3, float *pw4, complex *cdp, float *cl, 
	float *ct, float *t);
		
void compute_al_at (int wtype, int nlayers, int layern, float eps, float epsp, 
	float epss, float sigma, float sigp, float sigs, float *wrefp, float *wrefs,
	float fw, float *cl, float *ct, float *ql, float *qt, float *rho, 
	complex wpie, complex *al, complex *at);

void compute_prs (int nor, int *lobs, complex *al, complex *at, float *rho,
	complex *prs);

void compute_block_size (int np, int *block_size, int *nblock, int *left);

void compute_pwin (int wtype, int block_size, float bp, float dp, float w, 
	float decay, float pw1, float pw2, float pw3, float pw4, complex *pwin, 
	complex *p, complex *pp);

void compute_gl_gt_gam (int wtype, int nlayers, int block_size, complex *al, 
	complex *at, float *rho, complex *pp, complex *gl, complex *gt, 
	complex *gam);

void compute_alpha_betha (int lsource, int block_size, complex ats, float rhos, 
	complex *gl, complex *gt, complex *alpha, complex *betha);

void compute_sigmas (int wtype, int block_size, int lsource, float m1, 
	float m2, complex meu, complex *alpha, complex *betha, complex *p, 
	complex *gl, complex *gt, complex *gam, complex a2, complex a3, 
	complex *pwin, complex s1, complex s2, complex s4, complex *sigmau1, 
	complex *sigmau2, complex *sigmad1, complex *sigmad2);

void compute_moment_tensor (int wtype, float phi, float lambda, float delta, 
	float phis, float m0, float *m1, float *m2, float *m3);

void parameter_interpolation (int nlayers, int *intlayers, int *nintlayers,
	float *intlayth, float *cl, float *ql, float *ct, float *qt, float *rho, 
	float *t) ;

void random_velocity_layers (int *nlayers, int *lsource, int nrand_layers,
	float sdcl, float sdct, float layer, float zlayer, float *cl, float *ql, 
	float *ct, float *qt, float *rho, float *t);

void apply_earth_flattening (int nlayers, float z0, float *cl, float *ct, 
	float *rho, float *t);
*******************************************************************************
Notes:

******************************************************************************/
/********************************End Self-Documentation***********************/


/******************************************************************************

				Subroutine to compute auxiliary arrays

******************************************************************************/
void compute_w_aux_arrays (int wtype, int layern, int nlayers, int nor, 
	int lsource, int *np, int *block_size, int *nblock, int *left, float fref,
	float wrefp, float wrefs, int *lobs, float tsec, float p2w, float fs, 
	float xmax, float w, float decay, float epsp, float epss, float sigp, 
	float sigs, float *pw1, float *pw2, float *pw3, float *pw4, float *dp, 
	float *fp, complex wpie, complex *cdp, float *cl, float *ct, float*ql, 
	float *qt, float *rho, float *t, complex *al, complex *at, complex *prs)
/******************************************************************************
Input:
layern			
nlayers		number of reflecting layers
tsec		time series length in seconds (the maximum frequency
			of the seismogram will be (nw/tsec) Hz
w			frequency in rad/sec
decay		decay factor, used to avoid time series wraparound
epsp
epss
sigp
sigs 
pw1			left lower side window ray paramter for tapering
pw2			right lower side window ray paramter for tapering
pw3			left upper side window ray paramter for tapering
pw4			right upper side window ray paramter for tapering
cl			array[nlayers] of compressional wave velocities km/s
ct			array[nlayers] of shear wave velocities km/s
ql			array[nlayers] of compressional Q-values
qt			array[nlayers] of shear Q-values
rho			array[nlayers] of densities g/cc
t			array[nlayers] of layer thicknesses km

Output:
******************************************************************************/
{
	float wref;			/* reference frequency in rads/sec */
	float eps,sigma,e4,dk;		/* auxiliary variables */

	/* initialize varaibles */ 
	eps=0.001;
	sigma=0.1;
	dk=0.8;
	e4=0.8;
	wref=2*PI*fref;

	/* compute al and at */
	compute_al_at (wtype, nlayers, layern, eps, epsp, epss, sigma, sigp, 
		sigs, &wrefp, &wrefs, fref, cl, ct, ql, qt, rho, wpie, al, at);

	/* compute prs for every receiver */
	compute_prs (nor, lobs, al, at, rho, prs);

	/* compute number of slowness, slowness increment and maximum slowness */
	compute_slowness (wtype, lsource, w, p2w, np, fp, fs, e4, dk, decay, 
		xmax, dp, pw1, pw2, pw3, pw4, cdp, cl, ct, t);

	/* compute number of blocks and block size for matrix computations */
	compute_block_size (*np, block_size, nblock, left);
}

/******************************************************************************

		Subroutine to compute auxiliary arrays

******************************************************************************/
void compute_p_aux_arrays (int wtype, int nlayers, int lsource, int block_size,
	float bp, float dp, float w, float decay, float pw1, float pw2,
	float pw3, float pw4, complex *pwin, float m1, float m2, float m3, float h1,
	float h2, complex wpie, complex *divfac, complex *p, complex *pp,
	complex *al, complex *at, float *rho, complex *gl, complex *gt,
	complex *gam, complex *alpha, complex *betha, complex *sigmad1,
	complex *sigmad2, complex *sigmau1, complex *sigmau2) 
/******************************************************************************
Input:







Output:


******************************************************************************/
{
	float rhos;
	complex als,ats;
	complex alphasq,bethasq;
	complex meu,lambda;
	complex s1,s2,s4,a1,a2,a3;
	
	/* save al, at and rho values for the source layer */
	als=al[lsource-1];
	ats=at[lsource-1];
	rhos=rho[lsource-1];

	/* compute complex auxiliary variables */
	alphasq=cdiv(cmplx(1.0,0.0),cmul(als,als));
	if ((ats.r==0.)&&(ats.i==0.)) bethasq=cmplx(0.0,0.0);
	else bethasq=cdiv(cmplx(1.0,0.0),cmul(ats,ats));
	meu=crmul(bethasq,rhos);
	if ((meu.r != 0.)||(meu.i !=0.)) a1=cdiv(cmplx(1.0,0.0),meu);
	else a1=cmplx(0.0,0.0);
	lambda=crmul(cmul(csub(alphasq,cmplx(2,0)),bethasq),rhos);
	*divfac=cmul(cmplx(0.0,1.0),wpie);
	s1=crmul(a1,m2);
	s2=crmul(crmul(cmul(als,als),1./rhos),m3);
	s4=cdiv(cmplx(h2,0.),*divfac);
	a2=cdiv(cmplx(h1,0.),*divfac);
	a3=csub(cmul(lambda,s2),cmplx(m1,0.));

	/* compute pwin */
	compute_pwin (wtype, block_size, bp, dp, w, decay, pw1, pw2, pw3, pw4, 
		pwin, p, pp);

	/* compute gl and gt and gam */
	compute_gl_gt_gam (wtype, nlayers, block_size, al, at, rho, pp, gl, gt, 
		gam);

	/* compute alpha and betha */
	compute_alpha_betha (lsource, block_size, ats, rhos, gl, gt, alpha, betha);

	/* compute sigmau1,sigmau2,sigmad1,sigmad2 */
	compute_sigmas (wtype, block_size, lsource, m1, m2, meu, alpha, betha, p, 
		gl, gt, gam, a2, a3, pwin, s1, s2, s4, sigmau1, sigmau2, sigmad1,
		sigmad2); 
}
	
/******************************************************************************

	Subroutine to determine the type (acoustic or elastic) of sources
				and receivers 

******************************************************************************/
void source_receiver_type (int nor, int nlayers, int lsource, int  *lobs, 
	float *cl, float *ct, int *acoustic, int *flag)
/******************************************************************************
Input:
nor	    	number of recievers
nlayers		number of reflecting layers
cl 	    	array[nlayers] of compressional velocities km/s
ct 	    	array[nlayers] of shear velocities km/s
acoust  	array[nor] of receiver-type flags
		=1 for receiver acoustic
		=2 for first receiver elastic below source
		=3 for receiver elastic below source
		=4 for non-physical conditions ??
flag    	array[nor] of flags ??
		=1 if the compressional and shear velocities in the
		receiver layer are the same as their corresponding
		values for the last layer

Output:
acoust   	array[nor] of updated flags according to reciver type
flag
******************************************************************************/
{
	int k=0,iz;			/* loop counters */
	int jkl,lmn;			/* auxiliary indices */

	/* loop over the receivers */
	for (iz=0; iz<nor; iz++) {
		jkl=lsource-1;
		lmn=lobs[iz]-1;
		acoustic[iz]=0;
		flag[iz]=0;

		/* test condition for flag */
		if ((cl[nlayers-1]==cl[lmn])&&(ct[nlayers-1]==ct[lmn]))
			flag[iz]=1;

		/* test condition for source and receiver type */
		if (ct[jkl]==0.0) {		/* source acoustic */

			if (ct[lmn]==0.0) {	/* receiver acoustic */
				acoustic[iz]=1;	
			} else {
				if (lmn>lsource-1) k++;

				if (k==1) { 	/* 1st elastic rec below sorce*/
					acoustic[iz]=2;	
				} else {    	/* elastic receiver below surface */
					acoustic[iz]=3;	
				}
			}

		} else {    			/* source elastic */

			if (ct[lmn]==0.0) { /* receiver acoustic */

				acoustic[iz]=1;

			}
		}
	}
	
	/* only for Po/So */
	if ((nor<=1)&&(lsource>lobs[0])&&(ct[0]==0.0)) acoustic[0]=4;
}

/******************************************************************************

		Subroutine to compute the number of slowness values, final 
				slownesses, slownes increment and np3

******************************************************************************/
void compute_slowness (int wtype, int lsource, float w, float p2w, int *np, 
	float *fp, float fs, float e4, float dk, float decay, float xmax, float *dp,
	float *pw1, float *pw2, float *pw3, float *pw4, complex *cdp, float *cl, 
	float *ct, float *t)
/*****************************************************************************
Input:
lsource 		layer on top of which the source is located 
w			frequency in rad/sec
p2w			maximum ray parameter value to which the synthetics
			are computed
fs			sampling parameter, usually between 0.07 and 0.12
e4
dk
decay   		decay factor to avoid time series wraparound
xmax    		maximum range
dp			pointer to ray parameter increment
np			pointer to number of ray parameters
pw1			pointer to left lower side window ray parameter
pw2			pointer to right lower side window ray parameter
pw3			pointer to left higher side window ray prameter		
pw4			pointer to right higher side window ray parameter
cdp			pointer to complex ray parameter increment
cl			array of compressional velocities
t			array of layer thicknesses

Output:
			calculated values for np,bp,fp,dp,pw1,pw2,pw3 and pw4
******************************************************************************/
{
	int numw;		/* number of frequencies */
	int index=0;
	float fq;		/* frequency in hertz */
	float fp1,fp2,dkt;   	/* auxiliary variables */
	float bp;		/* beginning and final ray parameters */

	/* initialize variables */
	fq=w/(PI*2.0);
	bp=0.0;
	if (wtype==1) index=lsource-1;
	else if (wtype==2) index=0;

	/* compute fp1 and fp2 */
	if (ct[index]<=0.0005) {
		fp1=sqrt(pow(1.0/cl[lsource-1],2)+pow(50.0/(w*(t[lsource-1]+
				0.00001)),2));
		fp2=sqrt(pow(1.0/cl[lsource-2],2)+pow(50.0/(w*(t[lsource-2]+
				0.00001)),2));
	} else {
		fp1=sqrt(pow(1.0/ct[lsource-1],2)+pow(50.0/(w*(t[lsource-1]+
				0.00001)),2));
		fp2=sqrt(pow(1.0/ct[lsource-2],2)+pow(50.0/(w*(t[lsource-2]+
				0.00001)),2));
	}

	/* compute the slowness parameters */
	*fp=MAX(fp1,fp2);
	dkt=dk*PI/xmax;
	*dp=(dkt/w)*pow((1.0+fq/fs),e4/2.);
	*dp=MIN(0.002,*dp);
	if ((p2w>0.0)&&(*fp>=p2w)) *fp=p2w;
	*np=*fp/ *dp;
	*fp=*np* *dp; 
	*cdp=cmplx(*dp,-(*dp)*(decay/w));

	/* if required, use default values for window slowness parameters */
	if ((*pw2<=0.) && (*pw3<=0.)) {
		numw=*np*0.15-1;
		*pw1=bp;
		*pw2=bp+*dp*numw;
		*pw3=*fp-*dp*numw;
		*pw4=*fp;
	}
}

/******************************************************************************

			Subroutine to compute the al and at arrays 

******************************************************************************/
void compute_al_at (int wtype, int nlayers, int layern, float eps, float epsp, 
	float epss, float sigma, float sigp, float sigs, float *wrefp, float *wrefs,
	float fw, float *cl, float *ct, float *ql, float *qt, float *rho, 
	complex wpie, complex *al, complex *at)
/******************************************************************************
Input:
nlayers 		number of reflecting layers
nlayern	
eps
epsp
spss
wpie    		complex frequency
sigma