auxil.c

su 的源代码库

			sigmau1[ip]=cneg(cmul(cadd(cadd(cmul(d1,s1),cmul(d2,s2)),
				cadd(cmul(d3,s3),cmul(d4,s4))),pwin[ip]));
			sigmau2[ip]=cneg(cmul(cadd(cadd(cmul(d5,s1),cmul(d6,s2)),
				cadd(cmul(d7,s3),cmul(d8,s4))),pwin[ip]));
			sigmad1[ip]=cmul(cadd(csub(cmul(d1,s1),cmul(d2,s2)),
				csub(cmul(d4,s4),cmul(d3,s3))),pwin[ip]);
			sigmad2[ip]=cmul(cadd(csub(cmul(d6,s2),cmul(d5,s1)),
				csub(cmul(d7,s3),cmul(d8,s4))),pwin[ip]);

		} else if (wtype==2) {
		
			/* avoid p and gt being zero */
			if ((p[ip].r==0.)&&(p[ip].i==0.)) p[ip]=cmplx(0.0001,0.0);
			if ((gt[ijk1].r==0.)&&(gt[ijk1].i==0.)) gt[ijk1]=cmplx(0.0001,0.0);

			/* compute sigmas */
			sigmau2[ip]=cmplx(0.0,0.0);
			sigmad2[ip]=cmplx(0.0,0.0);
			sigmau1[ip]=cneg(cmul(cadd(crmul(cinv(cmul(meu,p[ip])),m2/2.),
				crmul(cinv(cmul(meu,gt[ijk1])),m1/2.)),pwin[ip]));
			sigmad1[ip]=sigmau1[ip];
		}
	}
}

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

			Subroutine to compute moment tensor components 

******************************************************************************/
void compute_moment_tensor (int wtype, float phi, float lambda, float delta, 
	float phis, float m0, float *m1, float *m2, float *m3)
/******************************************************************************
Input:
phi			azimuth of the receiver location (degrees)
lambda		rake (degrees)
delta		dip (degrees)
phis		azimuth of the fault plane (degrees)

Output
m1			[1][1] component of moment tensor
m2			[1][2] and [2][1] components of moment tensor
m3			[2][2] component of moment tensor
******************************************************************************/ 
{
	float a,b,c,d,e,f;			/* auxiliary variables */

	/* update fault plane azimuth */
	phis -=phi;

	/* convert angles to radians */
	delta *=PI/180.;
	lambda *=PI/180.;
	phis *=PI/180.;

	if (wtype==1) {
		/* compute auxiliary variables */
		a=sin(phis);
		b=a*a;
		c=sin(delta)*cos(lambda)*sin(2.*phis);
		d=sin(2.*delta)*sin(lambda);
		e=cos(delta)*cos(lambda)*cos(phis);
		f=cos(2.*delta)*sin(lambda)*a;

		/* compute the moment tensor components */
		*m1=-m0*(c+d*b);
		*m2=-m0*(e+f);
		*m3=m0*d;
	} else if (wtype==2) {
		*m1=m0*(sin(delta)*cos(lambda)*cos(2.*phis)+0.5*sin(2.*delta)*
			sin(lambda)*sin(2.*phis));
		*m2=-m0*(cos(delta)*cos(lambda)*sin(phis)-cos(2.*delta)*
			sin(lambda)*cos(phis));
		*m3=0.0;
	}
}

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

		Subroutine to expand interpolated layers and update
				input parameters

******************************************************************************/
void parameter_interpolation (int nlayers, int *intlayers, int *nintlayers,
	float *intlayth, float *cl, float *ql, float *ct, float *qt, float *rho, 
	float *t) 
/******************************************************************************
Input:
nlayers		number of total reflecting layers
nlint		number of times layer interpolation is required
intlayers	array[nlint] of layers on top of which interpolation is required
nintlayers	array[nlint] of number of layers to interpolate each time
intlayth	array[nlint] of layer thicknesses to interpolate each time
cl			array[nlayers] of compressional velocities (Km/s)
ql			array[nlayers] of compressional q values
ct			array[nlayers] of shear velocities (Km/s)
qt			array[nlayers] of shear q values
rho			array[nlayers] of densities
t			array[nlayers] of absolute depths (Km)

Output:
			same arrays with interpolated layers
*******************************************************************************
Note:
******************************************************************************/
{
	int i=0,ii=0,j;			/* loop counters */
	int nintlay;			/* number of layer to interpolate */
	float tintlay;			/* thickness at which to interpolate */
	float incrcl;			/* interpolation step for cl */
	float incrql;			/* interpolation step for ql */
	float incrct;			/* interpolation step for ct */
	float incrqt;			/* interpolation step for qt */
	float incrrho;			/* interpolation step for rho */
	float crust=0.0;		/* auxiliary variable */

	while (i<nlayers) {
	
		if (intlayers[ii] == i) { /* interpolation requested for this layer */

			/* get number of layers to interpolate */
			nintlay=nintlayers[ii];
			tintlay=intlayth[ii];

			/* check if thickness is negative */
			if (tintlay<0.0) {
				tintlay=-tintlay-crust;
			}
			crust +=tintlay;

			/* if required, adjust input parameters */
			if (ct[i+nintlay]==-1.0) 
				ct[i+nintlay]=cl[i+nintlay]/sqrt(3.0);
			if (qt[i+nintlay]==-1.0) 
				qt[i+nintlay]=ql[i+nintlay]/2.;
			if (rho[i+nintlay-1]==-1.0) 
				rho[i+nintlay]=(cl[i+nintlay]+1.5)/3.;

			/* compute the interpolation steps */
			incrcl=(cl[i+nintlay]-cl[i-1])/(nintlay+1);
			incrql=(ql[i+nintlay]-ql[i-1])/(nintlay+1);
			incrct=(ct[i+nintlay]-ct[i-1])/(nintlay+1);
			incrqt=(qt[i+nintlay]-qt[i-1])/(nintlay+1);
			incrrho=(rho[i+nintlay]-rho[i-1])/(nintlay+1);

			/* do the actual interpolation */
			for (j=0; j<nintlay; j++) {
				cl[j+i-1]=cl[i-1]+(j+1)*incrcl;
				ql[j+i-1]=ql[i-1]+(j+1)*incrql;
				ct[j+i-1]=ct[i-1]+(j+1)*incrct;
				qt[j+i-1]=qt[i-1]+(j+1)*incrqt;
				rho[j+i-1]=rho[i-1]+(j+1)*incrrho;
				t[j+i-1]=tintlay/nintlay;
			}

			/* update counter by number of interopolated layers */
			i +=nintlay;
			ii++;

		} else {		/* no layer interpolation requested */

			/* if required, adjust input parameters */	
			if (ct[i]==-1.0) 
				ct[i]=cl[i]/sqrt(3.0);
			if (qt[i]==-1.0)
				qt[i]=ql[i]/2.;
			if (rho[i]==-1.0)
				rho[i]=(cl[i]+1.5)/3.;
			if (t[i]<0.0)
				t[i]=-t[i]-crust;
			if (i<nlayers-1)
				crust +=t[i];
			else 
				t[i]=10000.0;	/* make last layer a subspace */
			i++;
		}
	}
}

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

					Compute random velocity layers

******************************************************************************/
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)
/******************************************************************************
Input:
nlayers		pointer to number of reflecting layers
lsource		pointer to layer on top of which the source is located
sdcl		compressional wave velocity standar deviation
sdct		shear wave velocity standar deviation
layer		layer number to star computing random velocity layers
zlayer		layer thickness above which to insert the random layers
cl			array[nlayers] of compressional wave velocities
ql			array[nlayers] of compressional wave Q-values
ct			array[nlayers] of shear wave velocities
qt			array[nlayers] of shear wave Q-values
rho			array[nlayers] of densities
t			array[nlayers] of layer thicknesses

Output:	
			same arrays with the random velocity layers information
			included and nlayers and lsource updated to reflect the
			insertion of the random velocity layers
*******************************************************************************
Note:
This subroutine computes a "layer" number of layers with random velocities
with means cl(compressional) and ct(shear) and standard deviations sdcl and 
sdct. The layers are inserted after the layer thickness zlayer and numbered
"layer" and onwards
******************************************************************************/
{
	int i=0,il,j=0;		/* loop counters */
	int icnt;		/* count number of layers to insert */
	int ijk;
	int ksource=0;		/* index to update source layer number*/
	int nl=*nlayers;	/* current number of layers */
	int nlmax;
	float d1,d2,tlast,x,zl;			/* auxiliary varibles */
	float cls,cts,qls,qts,rhos;		/* auxiliary variables */
	float *cll,*ctt,*qll;			/* scratch arrays */
	float *qtt,*rhoo,*tt;			/* more scratch arrays */

	fprintf (stderr, "before rand: nlayers=%d ",nl);

	/* allocate working space */
	cll=alloc1float(nl);
	ctt=alloc1float(nl);
	qll=alloc1float(nl);
	qtt=alloc1float(nl);
	rhoo=alloc1float(nl);
	tt=alloc1float(nl);
	
	/* initialize variables */
	d1=2.0*sqrt(3.0)*(sdcl/100.0);
	d2=2.0*sqrt(3.0)*(sdct/100.0);

	/* save information from last layer */
	cls=cl[nl-1];
	cts=ct[nl-1];
	qls=ql[nl-1];
	qts=qt[nl-1];
	rhos=rho[nl-1];
	nlmax=nrand_layers+nl-layer+1;

	/* main loop to compute the random layers */
	for (il=layer-1; il<nl-1; il++) {

		if (t[il]<=zlayer) {
			icnt=1;		/* compute just one layer */
			zl=t[il];	/* with thickness t[il] */
		} else {
			tlast=fmod(t[il],zlayer);
			icnt=t[il]/zlayer;				/* make icnt layers with */
			zl=zlayer+tlast/(float)icnt;	/* thickness t[il]/zlayer */
		}

		if (il==*lsource-1) ksource=j;
		for (ijk=0; ijk<icnt; ijk++) {
			tt[i]=zl;					/* all layers, same thickness */
			x=franuni()-0.5;
			cll[i]=cl[il]*(1.0+d1*x);	/* random velocities */
			ctt[i]=ct[il]*(1.0+d2*x);
			qll[i]=ql[il];				/* same q's and rho */
			qtt[i]=qt[il];
			rhoo[i]=rho[il];
			
			/* if maximum number of layers reached, exit inner loop */
			if (i==nlmax) break; 

			/* update counters */
			i++;
			if (il<*lsource-1) j++;
		}

		/* if maximum number of layers reached, exit outer loop */
		if (i==nlmax) break; 	
	}

	/* update number of layers and source layer index */ 
	nl=i+layer-1;
	*lsource +=ksource-1;

	/* restore information for last layer */
	cl[nl-1]=cls;
	ql[nl-1]=qls;
	ct[nl-1]=cts;
	qt[nl-1]=qts;
	rho[nl-1]=rhos;

	/* main loop to insert the computed information for the random layers */
	for (il=layer-1, i=0; il<nl-1; il++, i++) {
		cl[il]=cll[i];
		ct[il]=ctt[i];
		rho[il]=rhoo[i];
		ql[il]=qll[i];
		qt[il]=qtt[i];
		t[il]=tt[i];
		if (il==*lsource-1) {
			cl[il]=cl[il-1];
			ct[il]=ct[il-1];
			ql[il]=ql[il-1];
			qt[il]=qt[il-1];
			rho[il]=rho[il-1];
		}
	}
	t[nl-1]=10000.0;	/* make last layer a subspace */

	/* update pointer for number of layers */
	*nlayers=nl;
	fprintf (stderr, "after rand: nlayers=%d\n",nl);

	/* free allocated space */
	free1float(cll);
	free1float(ctt);
	free1float(qll);
	free1float(qtt);
	free1float(rhoo);
	free1float(tt);
}


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

		Subroutine to apply flattening earth approximation

******************************************************************************/
void apply_earth_flattening (int nlayers, float z0, float *cl, float *ct, 
	float *rho, float *t)
/******************************************************************************
Input:
nlayer		number of reflecting layers
z0			first layer depth
cl			array[nlayers] of compressional wave velocities
ct			array[nlayers] of shear wave velocities
rho			array[nlayers] of densities
t			array[nlayers] of layer thicknesses

Output:	
			same arrays with the correction applied
			(the thicksnesses array, t, is not modified)
*******************************************************************************
Note:
This subroutine applies the earth flattening correction after Chapman, 1973
******************************************************************************/
{
	int il;				/* loop counter */
	float y0,w0;		/* auxiliary variables */
	float z=z0;			/* layer depth */
	float *zz;			/* scratch array */

	/* allocate working space */
	zz=alloc1float(nlayers);

	/* create and array of absolute depths from the thicknesses */
	zz[0]=z;
	for (il=0; il<nlayers-1; il++) {
		zz[il+1]=z+t[il];
		z=zz[il+1];
	}

	/* apply the earth flattening correction */
	for (il=0; il<nlayers; il++) {

		/* compute the correction */
		y0=RSO/(RSO-zz[il]);
		zz[il]=RSO*log(y0);
		w0=RSO/(RSO-zz[il]);

		/* apply the correction */
		cl[il] *=w0;
		ct[il] *=w0;
		rho[il] /=w0;
	}

	/* update the thicknesses */
	for (il=1; il<nlayers; il++) t[il-1]=zz[il]-zz[il-1];

 	/* make last layer a half space*/
	t[nlayers-1]=10000.0;

	/* free allocated space */
	free1float(zz);
}