shreflect.c

su 的源代码库

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

#include "par.h"
#include "Reflect/reflpsvsh.h"

#define BLK 100
#define IMAX 300
#define SZERO 0.000000000001

/*********************** self documentation **********************/
/******************************************************************************
SH_REFLECTIVITIES -	Subroutines to do matrix propagation as the hard core of 
						a reflectivity modeling code for SH wavefields

sh_reflectivities		applies matrix propagation to compute the SH 
						reflectivity response of a horizontally stratified 
						earth model

*******************************************************************************
Function Prototypes:
void sh_reflectivities (int int_type, int verbose, int wtype, int nw,
	int nlayers, int nx, int layern, int nor, int np, float bp1, float m1,
	float m2, float m3, float h1, float h2, int lsource, float bx, float dx,
	float xmax, float decay, float fref, float wrefp, float wrefs, float tsec,
	float p2w, float fs, float epsp, float epss, float sigp, float sigs,
	float pw1, float pw2, float pw3, float pw4, int *flag, int *lobs,
	float *rho, float *t, float *cl, float *ct, float *ql, float *qt,
	complex ***response1, FILE *outfp);
*******************************************************************************
sh_reflectivities
Input:
nw  			number of frequencies
nlayers 		number of reflecting layers
nx  			number of ranges
nor 			number of receivers
lsource 		layer on top of which the source is located
bx  			beginning range
fx  			final range
dx  			range increment
flag			array[nor] of flags for source type
rho 			array[nlayers] of densities
t   			array[nlayers] of layer thicknesses

Output
response1   	array[nw][nor][nx] of updated tangential field response
*******************************************************************************
Credits:



******************************************************************************/
/**************** end self doc ********************************/



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

		Subroutine to compute the SH propagation matrices for a single
							frequency step

******************************************************************************/
void sh_reflectivities (int int_type, int verbose, int wtype, int nw,
	int nlayers, int nx, int layern, int nor, int np, float bp1, float m1,
	float m2, float m3, float h1, float h2, int lsource, float bx, float dx,
	float xmax, float decay, float fref, float wrefp, float wrefs, float tsec,
	float p2w, float fs, float epsp, float epss, float sigp, float sigs,
	float pw1, float pw2, float pw3, float pw4, int *flag, int *lobs,
	float *rho, float *t, float *cl, float *ct, float *ql, float *qt,
	complex ***response1, FILE *outfp)
/******************************************************************************
Input:
nw  			number of frequencies
nlayers 		number of reflecting layers
nx  			number of ranges
nor 			number of receivers
lsource 		layer on top of which the source is located
bx  			beginning range
fx  			final range
dx  			range increment
flag			array[nor] of flags for source type
rho 			array[nlayers] of densities
t   			array[nlayers] of layer thicknesses

Output
response1   	array[nw][nor][nx] of updated tangencialfield response
******************************************************************************/
{
	int iw,iz,ip,ix;		/* loop counters */
	int ijk1=0,ik1,ik2,il,jl;	/* auxiliary indices */
	int ijk=0,ijk2,jj=0;		/* more auxiliary indices */
	int lrec;			
	float bp=0.0;		/* beginning ray parameter */
	float fp;		/* final ray parameter */
	float dp;		/* ray parameter increment */
	int nblock;		/* number of blocks to process */
	int block_size;		/* size of processing blocks */
	int left;		/* remainder */
	int *prv,*ibl;		/* scratch arrays */
	float x1;
	float w;
	complex wpie,cdp;	/* complex frecuency and ray parameter inc */
		
	/* complex auxiliary variables */
	complex x,y,t1,divfac;
	complex det,psq,gl1,gt1,gl2,gt2,gam1,gam2,at1,at2,al1,al2;
	complex tem,temr,sig,sigh,rho1,rho2;
	complex rvrb111,rvrb112,rvrb121,rvrb122,rvrb211;
	complex ewigh11,ewigh22,wrn011,ewrd11;
	complex ewd0211,ewd0212,ewtu211;
	complex rtb111,rtb112,rtb121,rtb122,rtb211;
	complex rtb212,rtb221,rtb222,rtb311,rtb312,rtb321,rtb322;
	complex func1,func2;
	complex tun0p11,tun0p21,rd0np11,rd0np12,rd0np21;
	complex rd0np22,td0np11,td0np12;
	complex rd11,td11,tu11,ru11;
	complex fact11,fact12,vuzs1,vuzs2,vdzs1,vdzs2;

	/* complex scratch arrays */
	complex *rd0n11,*rd0n12,*rd0n21,*rd0n22,*td0n11,*td0n12;
	complex *td0n21,*td0n22,*tun011,*tun012,*tun021,*tun022;
	complex *run011,*run012,*run021,*run022;
	complex *rusf21,*rusf22,*tusf11,*tusf12,*tusf21,*tusf22;
	complex *rdfs11,*rdfs12,*rdfs21,*rdfs22,*rurf11,*rurf12;
	complex *rurf21,*rurf22,*turf11,*turf12,*turf21,*turf22;
	complex *rdfr11,*rdfr12,*rdfr21,*rdfr22,*tdfr11,*tdfr12;
	complex *tdfr21,*tdfr22,*r11,*r12,*r21,*r22;
	complex *texp,*vos1,*vos2,*vos3,*up;
	complex *sigmau1,*sigmau2,*sigmad1,*sigmad2,*alpha,*betha;
	complex *al=NULL,*at=NULL,*gl,*gt,*gam,*p,*pp,*prs=NULL,*pwin;
	
	/* complex pointers */
	complex *rusf11,*rusf12;
	complex *tusr11,*tusr12,*tusr21,*tusr22;
	complex *rusr11,*rusr12,*rusr21,*rusr22;
	complex *tdsr11,*tdsr12,*tdsr21,*tdsr22;
	complex *rdnr11,*rdnr12,*rdnr21,*rdnr22;

	/* allocate working space */

	/* allocate working space */
	up=alloc1complex(BLK);
	rd0n11=alloc1complex(BLK);rd0n12=alloc1complex(BLK);
	rd0n21=alloc1complex(BLK);rd0n22=alloc1complex(BLK);
	td0n11=alloc1complex(BLK);td0n12=alloc1complex(BLK);
	td0n21=alloc1complex(BLK);td0n22=alloc1complex(BLK);
	tun011=alloc1complex(BLK);tun012=alloc1complex(BLK);
	tun021=alloc1complex(BLK);tun022=alloc1complex(BLK);
	run011=alloc1complex(BLK);run012=alloc1complex(BLK);
	run021=alloc1complex(BLK);run022=alloc1complex(BLK);
	tusf11=alloc1complex(BLK);tusf12=alloc1complex(BLK);
	tusf21=alloc1complex(BLK);tusf22=alloc1complex(BLK);
	rdfs11=alloc1complex(BLK);rdfs12=alloc1complex(BLK);
	rdfs21=alloc1complex(BLK);rdfs22=alloc1complex(BLK);
	rurf11=alloc1complex(IMAX);rurf12=alloc1complex(IMAX);
	rurf21=alloc1complex(IMAX);rurf22=alloc1complex(IMAX);
	turf11=alloc1complex(IMAX);turf12=alloc1complex(IMAX);
	turf21=alloc1complex(IMAX);turf22=alloc1complex(IMAX);
	rdfr11=alloc1complex(IMAX);rdfr12=alloc1complex(IMAX);
	rdfr21=alloc1complex(IMAX);rdfr22=alloc1complex(IMAX);
	tdfr11=alloc1complex(IMAX);tdfr12=alloc1complex(IMAX);
	tdfr21=alloc1complex(IMAX);tdfr22=alloc1complex(IMAX);
	rusf21=alloc1complex(BLK);rusf22=alloc1complex(BLK);
	vos1=alloc1complex(BLK);vos2=alloc1complex(BLK);vos3=alloc1complex(BLK);
	r11=alloc1complex(BLK);r12=alloc1complex(BLK);
	r21=alloc1complex(BLK);r22=alloc1complex(BLK);
	texp=alloc1complex(BLK);

	/* allocate other working space */
	prv=alloc1int(nor);
	ibl=alloc1int(nor);
	/* ibl=alloc1int(BLK); */
    p=alloc1complex(BLK);
    pp=alloc1complex(BLK);
    pwin=alloc1complex(BLK);
    /* gl=alloc1complex(BLK); */
    gl=alloc1complex(BLK*nlayers);
    gt=alloc1complex(BLK*nlayers);
    gam=alloc1complex(BLK*nlayers);
    alpha=alloc1complex(BLK);
    betha=alloc1complex(BLK);
    sigmau1=alloc1complex(BLK);
    sigmau2=alloc1complex(BLK);
    sigmad1=alloc1complex(BLK);
    sigmad2=alloc1complex(BLK);
	at=alloc1complex(nlayers);
	al=alloc1complex(nlayers);
	/* prs=ealloc1complex(nlayers); */
	prs=alloc1complex(nor);

	/* initialize pointers */
	rusf11=alpha;
	rusf12=betha;
	tusr11=tun011;
	tusr12=tun012;
	tusr21=tun021;
	tusr22=tun022;
	rusr11=run011;
	rusr12=run012;
	rusr21=run021;
	rusr22=run022;
	tdsr11=td0n11;
	tdsr12=td0n12;
	tdsr21=td0n21;
	tdsr22=td0n22;
	rdnr11=r11;
	rdnr12=r12;
	rdnr21=r21;
	rdnr22=r22;

	/* assign initial value to working variables */
	tem=cdiv(cmplx(1.0,0.0),cexp(crmul(cmplx(0.0,1.0),PI/4)));
	ewigh11=cmplx(0.0,0.0);
	ewigh22=cmplx(0.0,0.0);

	/* initialize complex output arrays */
    for (iw=0; iw<nw; iw++)
	for (ix=0; ix<nx; ix++)
	    for (iz=0; iz<nor; iz++) {
		response1[iw][ix][iz]=cmplx(0.0,0.0);
	    }


	/**************************************************************************
	*					Main loop over frequencies			    *	
	**************************************************************************/
	for (iw=0; iw<nw; iw++) {
		w=(iw+1)*PI*2/tsec;		/* w in radians/sec */
		wpie=cmplx(w,decay);		/* complex w */

		/* compute frequency-dependent auxiliary variables for reflectivity */
		compute_w_aux_arrays (wtype, layern, nlayers, nor, lsource,
			&np, &block_size, &nblock, &left, fref, wrefp, wrefs, 
			lobs, tsec, p2w, fs, xmax, w, decay, epsp, epss, sigp, 
			sigs, &pw1, &pw2, &pw3, &pw4, &dp, &fp, wpie, &cdp, cl,
			ct, ql, qt, rho, t, al, at, prs);

		/**********************************************************************
		* 	loop over slopes to compute the reflectivities (by blocks)		  *
		**********************************************************************/
		while (1) {

			if ((nblock==1) && (left !=0)) block_size=left;

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

			/* if requested, output processing information */
			if (verbose==1||verbose==3) {
				fprintf(stderr,"%3d%6.1f%7.3f%7.3f%6d%7.3f%7.3f%7.3f"
				"%7.3f%7.3f\n",iw,w,bp,fp,np,dp,pw1,pw2,pw3,pw4);
			} if (verbose==2||verbose==3) {
				fprintf(outfp,"%3d%6.1f%7.3f%7.3f%6d%7.3f%7.3f%7.3f"
				"%7.3f%7.3f\n",iw,w,bp,fp,np,dp,pw1,pw2,pw3,pw4);
			}



			/******************************************************************
			* 			start the computation of the reflectivities 		*
			******************************************************************/
			/* loop over reflecting layers */
			/* for (il=0; il<nlayers; il++) { */
			ijk=0;
			for (il=0; il<nlayers-1; il++) {
				jl=il+1;
				ik1=il*block_size;
				ik2=(il+1)*block_size;
				al1=al[il];
				al2=al[jl];
				at1=at[il];
				at2=at[jl];
				rho1=cmplx(rho[il],0.0);
				rho2=cmplx(rho[jl],0.0);
				t1=cmplx(t[il],0.0);
				if (il==0) {
					for (ip=0; ip<block_size; ip++) {
						ijk1=ik1+ip;
						gt1=gt[ijk1];
						ewigh11=cexp(cmul(cmplx(0.0,1.0),cmul(wpie,cmul(gt1,t1))));
						rd0n11[ip]=cmplx(0.0,0.0);
						rd0n12[ip]=cmplx(0.0,0.0);
						rd0n21[ip]=cmplx(0.0,0.0);
						rd0n22[ip]=cmplx(0.0,0.0);

						/* free surface */
						td0n11[ip]=ewigh11;
						td0n12[ip]=cmplx(0.0,0.0);
						td0n21[ip]=cmplx(0.0,0.0);
						td0n22[ip]=cmplx(0.0,0.0);
						tun011[ip]=ewigh11;
						tun012[ip]=cmplx(0.0,0.0);
						tun021[ip]=cmplx(0.0,0.0);
						tun022[ip]=cmplx(0.0,0.0);
						run011[ip]=cmul(ewigh11,ewigh11);
						run012[ip]=cmplx(0.0,0.0);
						run021[ip]=cmplx(0.0,0.0);
						run022[ip]=cmplx(0.0,0.0);
					}	
				} else {
					if ((at1.r==at2.r)&&(at1.i==at2.i)) {
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							gt1=gt[ijk1];
							ewigh11=cexp(cmul(cmplx(0.0,1.0),cmul(wpie,cmul(gt1,t1))));
							tun0p11=cmul(tun011[ip],ewigh11);
							tun0p21=cmul(tun021[ip],ewigh11);
						
							/* -R/T coefficients betha1=betha2 */
							rtb111=cmul(ewigh11,run011[ip]);
							run011[ip]=cmul(rtb111,ewigh11);
							run012[ip]=cmplx(0.0,0.0);
							run021[ip]=cmplx(0.0,0.0);
							run022[ip]=cmplx(0.0,0.0);
							td0np11=cmul(ewigh11,td0n11[ip]);
							td0np12=cmul(ewigh11,td0n12[ip]);
							tun011[ip]=tun0p11;
							tun012[ip]=cmplx(0.0,0.0);
							tun021[ip]=tun0p21;
							tun022[ip]=cmplx(0.0,0.0);
							td0n11[ip]=td0np11;
							td0n12[ip]=td0np12;
							td0n21[ip]=cmplx(0.0,0.0);
							td0n22[ip]=cmplx(0.0,0.0);
						}
					} else {
		
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							ijk2=ik2+ip;
							psq=pp[ip];
							gl1=gl[ijk1];
							gt1=gt[ijk1];
							gam1=gam[ijk1];
							gl2=gl[ijk2];
							gt2=gt[ijk2];
							gam2=gam[ijk2];
							ewigh11=cexp(cmul(wpie,cmul(cmplx(0.0,1.0),cmul(t1,gt1))));
							y=cadd(cmul(rho1,cdiv(gt1,cmul(at1,at1))),
								cmul(rho2,cdiv(gt2,cmul(at2,at2))));
							td11=crmul(cmul(rho1,cdiv(gt1,cmul(y,
								cmul(at1,at1)))),2.);
							rd11=cdiv(csub(cmul(rho1,cdiv(at1,at1)),
								cmul(rho1,cdiv(gt2,cmul(at2,at2)))),y);
							tu11=crmul(cmul(rho2,cdiv(gt2,cmul(y,
								cmul(at2,at2)))),2.);
							ru11=cneg(rd11);
							ewrd11=cmul(ewigh11,rd11);
							wrn011=cmul(ewigh11,run011[ip]);
							rvrb111=cdiv(cmplx(1.0,0.0),csub(cmplx(1.0,0.0),cmul(wrn011,ewrd11)));
							rvrb211=cadd(cmplx(1.0,0.0),cmul(rvrb111,cmul(ewrd11,wrn011)));
							ewd0211=cmul(td0n11[ip],cmul(rvrb111,ewigh11));
							ewd0212=cmul(td0n12[ip],cmul(rvrb111,ewigh11));
							td0n11[ip]=cmul(td11,ewd0211);
							td0n12[ip]=cmul(td11,ewd0212);
							td0n21[ip]=cmplx(0.0,0.0);