auxil.c

su 的源代码库

sigp
sigs
cl			array of compressional wave velocities
ct			array of shear wave velocities
ql			array of compressional Q-values
qt			array of shear Q-values
rho			array of densities

Output:
al			array of ??
at			array of ??
******************************************************************************/
{
	int i;  					/* loop counter */
	float qwif,qwifs,qwifp;		/* auxilairy variables */
	float wref;    				/* reference frequency in rad/sec */	
	float wrefpp,wrefss;		/* reference p ands frequencies in rad/sec */	
	complex eiw,eiwp,eiws;		/* auxiliary variables */

	/* compute reference frequencies for p and s waves in rad/s */
	wrefpp=*wrefp;
	wrefss=*wrefs;
	wref=2.0*PI*fw;
	wrefpp *=2.0*PI;
	wrefss *=2.0*PI;

	/* compute the new p and s q values */
	qwif=(pow((eps*eps+wref*wref),sigma/2.)/(sin(sigma*PI/2.)*2.));
	qwifp=(pow((epsp*epsp+wrefpp*wrefpp),sigp/2.)/(sin(sigp*PI/2.)*2.));
	qwifs=(pow((epss*epss+wrefss*wrefss),sigs/2.)/(sin(sigs*PI/2.)*2.));

	/* compute wie */
	eiw=crpow(csub(cmplx(eps,0.0),cmul(cmplx(0.0,1.0),wpie)),sigma);

	/* check to see if layern was provided as an input */
	if (layern != 0) {
		eiwp=crpow(csub(cmplx(epsp,0.0),cmul(cmplx(0.0,1.0),wpie)),sigp);
		eiws=crpow(csub(cmplx(epss,0.0),cmul(cmplx(0.0,1.0),wpie)),sigs);

		/* compute al and at until just before layern */
		for (i=0; i<layern-1; i++) {
			al[i]=cmul(cmplx(1./cl[i],0.),cadd(cdiv(cmplx(qwif/ql[i],0.),
				eiw),cmplx(1.0,0.0)));
			if (ct[i]==0.0) {
				at[i]=cmplx(0.0,0.0);
			} else {
				at[i]=cmul(cmplx(1./ct[i],0.),cadd(cdiv(cmplx(qwif/qt[i],0.),
				eiw),cmplx(1.0,0.0)));
			}
		}

		/* compute al and at for the remaining layers */
		for (i=layern-1; i<nlayers; i++) {
			al[i]=cmul(cmplx(1./cl[i],0.),cadd(cdiv(cmplx(qwifp/ql[i],0.),
				eiwp),cmplx(1.0,0.0)));
			if (ct[i]==0.0) {
				at[i]=cmplx(0.0,0.0);
			} else {
				at[i]=cmul(cmplx(1./ct[i],0.),cadd(cdiv(cmplx(qwifs/qt[i],0.),
				eiws),cmplx(1.0,0.0)));
			}
		}

	} else {

		/* compute al and at for all layers */
		for (i=0; i<nlayers; i++) {
			al[i]=cmul(cmplx(1./cl[i],0.),cadd(cdiv(cmplx(qwif/ql[i],0.),
				eiw),cmplx(1.0,0.0)));
			if (ct[i]==0.0) {
				at[i]=cmplx(0.0,0.0);
			} else {
				at[i]=cmul(cmplx(1./ct[i],0.),cadd(cdiv(cmplx(qwif/qt[i],0.),
				eiw),cmplx(1.0,0.0)));
			}
		}
	}
	
	/* update pointers */
	*wrefp=wrefpp;
	*wrefs=wrefss;
}

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

		Subroutine to compute the prs array

******************************************************************************/
void compute_prs (int nor, int *lobs, complex *al, complex *at, float *rho,
	complex *prs)
/******************************************************************************
Input:
nor			number of receivers
lobs		array[nor] layers on top of which the receivers are locate
al			arrar....
at			arrar....
rho			array of densities

Output
prs			array  ....
******************************************************************************/
{
	int ijk;
	int ijk1;
	complex ctemp;
	float rp;

	/* loop over the receivers */
	for (ijk=0; ijk<nor; ijk++) {
		ijk1=lobs[ijk]-1;
		rp=at[ijk1].r;
		if (rp==0.0) {
			prs[ijk]=crmul(cdiv(cmplx(1.0,0.0),cmul(al[ijk1],al[ijk1])),rho[ijk1]);
		} else {
			ctemp=cmul(cdiv(al[ijk1],at[ijk1]),cdiv(al[ijk1],
				at[ijk1]));
			prs[ijk]=crmul(csub(cmplx(1.0,0.0),crmul(ctemp,4./3.)),rho[ijk1]);
		}
	}
}

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

	Subroutine to compute the number of blocks and block size for
			matrix computations 

******************************************************************************/
void compute_block_size (int np, int *block_size, int *nblock, int *left)
/******************************************************************************
Input:
np				Number of ray parameters
block_size		pointer to block size
nblock			pointer to number of blocks
left			pointer to remainder blocks

Output
				computed values for block_size, nblock and left
******************************************************************************/
{
	if (np<=BLK) {		 

		/* one block is enough and nothing is left */
		*block_size=np;
		*nblock=1;
		*left=0;

	} else {
		
		/* define block size and compute number of required blocks */
		*block_size=BLK;
		*left=np%BLK;
		*nblock=np/BLK+1;
	}
}

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

		Subroutine to compute the pwin array 

*****************************************************************************/
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)
/******************************************************************************
Input:
block_size	block size for matrix computations
bp			begin ray parameter
dp			ray parameter increment
pw1			left lower window ray parameter for tapering
pw2			right lower window ray parameter for tapering
pw3			left upper window ray parameter for tapering
pw4			right upper window ray parameter for tapering

Output:
p			array[block_size] of complex ray parameters
pp			array[block_size] of complex squared ray parameters
pwin		array[block_size] of windowed complex ray parameters
******************************************************************************/
{
	int ip;				/* loop counter			*/
	float rp;			/* current ray parameter	*/
	float win=0.0,win1,win2;	/* auxiliary variables		*/

	/* loop over blocks */
	for (ip=0; ip<block_size; ip++) {

		/* compute array of complex slownesses */
		rp=bp+ip*dp;
		if (wtype==1) {
			p[ip]=cmplx(rp,0.0);
		} else if (wtype==2) {
			p[ip]=cmplx(rp, -decay*rp/w);
		}
		pp[ip]=cmul(p[ip],p[ip]);

		/* compute Hanning tapered window */
		if ((rp>=pw1)&&(rp<pw2)) {
			win1=0.5*(1.+cos(PI*((rp-pw1)/(pw2-pw1)-1.)));
			win=win1;
		}
		if ((rp>=pw2)&&(rp<=pw3)) win=1.0;
		if ((rp>pw3)&&(rp<pw4)) {
			win2=0.5*(1.+cos(PI*((rp-pw3)/(pw4-pw3))));
			win=win2;
		}
		if ((rp<pw1)||(rp>=pw4)) win=0.0;

		/* multiply complex ray parameter by window to compute pwin */
		pwin[ip]=crmul(csqrt(p[ip]),win);
	}
}

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

	Subroutine to compute arrays gl, gt and gam

******************************************************************************/
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)
/******************************************************************************
Input:
nlayers 		number of reflecting layers
block_size	block size for matrix computations
al			array[nlayers] of ...
at			array[nlayers] of ...
rho			array[nlayers] of densities
pp			array[nlayers] of complex squared ray parameters

Output:
gl			array[block_size] of ...
gt			array[block_size] of ...
gam			array[block_size] of ...
*****************************************************************************/
{
	int il,ip,ijk=0;	/* loop counters */
	int ik1,ijk1;		/* auxiliary indices */
	float rho1;		/* current density */
	complex p1,p2;

	/* loop over layers */
	for (il=0; il<nlayers; il++) {
		p1=cmul(al[il],al[il]);
		p2=cmul(at[il],at[il]);
		rho1=rho[il];
		ik1=il*block_size;

		/* if p2 is not complex number zero */
		if ((p2.r !=0.0)||(p2.i !=0.0)) {
			for (ip=0; ip<block_size; ip++) {
				ijk1=ik1+ip;

				/* compute gl and gt */
				gl[ijk1]=csqrt(csub(p1,pp[ip]));
				gt[ijk1]=csqrt(csub(p2,pp[ip]));

				/* make imaginary parts of gl and gt positive */
				if (gl[ijk1].i<0.0) gl[ijk1].i *=-1.0;
				if (gt[ijk1].i<0.0) gt[ijk1].i *=-1.0;
				if (wtype==2) {
					/* make imaginary parts of gl and gt positive */
					if (gl[ijk1].i<0.0) gl[ijk1].r *=-1.0;
					if (gt[ijk1].i<0.0) gt[ijk1].r *=-1.0;
				}
	
				/* compute gam */
				gam[ijk1]=crmul(csub(cmplx(1.0,0.0),crmul(cdiv(pp[ip],p2),2.)),rho1);
			}

		} else {		/* p2 is complex number zero */

			for (ip=0; ip<block_size; ip++) {
				ijk1=ik1+ip;

				/* compute gl and set gt to zero */
				gl[ijk1]=csqrt(csub(p1,pp[ip]));
				gt[ijk1]=cmplx(0.0,0.0);

				/* if necessary, make imaginary part of gl positive */
				if (gl[ijk].i<0.0) gl[ijk1]=cneg(gl[ijk1]);
				gam[ijk1]=cmplx(rho1,0.0);
			}
		}
	}
}

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

		Subroutine to compute alpha and betha

******************************************************************************/
void compute_alpha_betha (int lsource, int block_size, complex ats, float rhos, 
	complex *gl, complex *gt, complex *alpha, complex *betha)
/******************************************************************************
Input:
lsource		layer on top of which the source is located
block_size	block size for matrix computaions
ats			al[lsource-1]
rhos		rho[lsource-1] (density at the source layer)
gl			complex array[nlayers] of ...
gt			complex array[nlayers] of ...

Output
alpha		complex array[block_size] of ...
betha		complex array[block_size] of ...
******************************************************************************/
{
	int ip;				/* loop counter */
	int ik1,ijk1;		/* auxiliary indices */

	ik1=(lsource-1)*block_size;

	if ((ats.r==0.0)&&(ats.i==0.0)) {
		for (ip=0; ip<block_size; ip++) {
			ijk1=ik1+ip;
			alpha[ip]=cdiv(cmplx(1.0,0.0),crmul(crmul(gl[ijk1],rhos),2.0));
			betha[ip]=cmplx(0.0,0.0);
		}
	} else {
		for (ip=0; ip<block_size; ip++) {
			ijk1=ik1+ip;
			alpha[ip]=cdiv(cmplx(1.0,0.0),crmul(crmul(gl[ijk1],rhos),2.0));
			betha[ip]=cdiv(cmplx(1.0,0.0),crmul(crmul(gt[ijk1],rhos),2.0));
		}
	}
}	

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

	Subroutine to compute arrays for sigmau1, sigmau2, sigmad1
						and sigmad2

******************************************************************************/
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) 
/******************************************************************************
Input:
block_size	block size for matrix computations
meu
alpha
betha
p			array of complex ray parameters
gl
gt
gam
a2			
a3
pwin
s1
s2
s4
sigmau1		pointer to sigmau1
sigmau2		pointer to sigmau2
sigmad1		pointer to sigmad1
sigmad2		pointer to sigmad2

Output:
			Computed values for sigmau1, sigmau2, sigmad1, sigmad2
******************************************************************************/
{
	int ip;								/* loop counter */
	int ijk1,ik1;						/* auxiliary indices */
	complex d1,d2,d3,d4,d5,d6,d7,d8,s3;	/* auxiliary variables */	

	/* compute auxiliary index */
	ik1=(lsource-1)*block_size;

	/* loop over layers in a block */	
	for (ip=0; ip<block_size; ip++) {
		ijk1=ik1+ip;

		if (wtype==1) {
			/* compute auxiliary variables */
			d1=crmul(cmul(cmul(meu,gl[ijk1]),cmul(alpha[ip],p[ip])),2.0);
			d2=cneg(cmul(gam[ijk1],alpha[ip]));
			d3=cmul(p[ip],alpha[ip]);
			d4=cneg(cmul(gl[ijk1],alpha[ip]));
			d5=cmul(gam[ijk1],betha[ip]);
			d6=crmul(cmul(cmul(meu,gt[ijk1]),cmul(betha[ip],p[ip])),2.0);
			d7=cmul(gt[ijk1],betha[ip]);
			d8=cmul(p[ip],betha[ip]);
			s3=cadd(a2,cmul(p[ip],a3));

			/* compute sigmas */