psvreflect.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***************************/
/******************************************************************************
REFLECTIVITIES -	Subroutines to do matrix propagation as the hard core of 
						a reflectivity modeling code for PSV and SH wavefields
psv_reflectivities		applies matrix propagation to compute the PSV 
							reflectivity response of a horizontally stratified 
							earth model
*******************************************************************************
Function Prototypes:
void psv_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 *acoustic, int *flag,
	int *lobs, float *rho, float *t, float *cl, float *ct, float *ql, float *qt,
	complex ***response1, complex ***response2, complex ***response3,
	FILE *outfp);
*******************************************************************************
psv_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
acoustic		array[nor] of flags for receiver type
flag			array[nor] of flags for source type
rho 			array[nlayers] of densities
t   			array[nlayers] of layer thicknesses

Output
response1   	array[nw][nx][nor] of computed pressure wavefield
response2   	array[nw][nx][nor] of computed radial wavefield
response3   	array[nw][nx][nor] of computed vertical wavefield
*******************************************************************************
Credits:


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





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

		Subroutine to compute the propagation matrices for a single
							frequency step

******************************************************************************/
void psv_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 *acoustic, int *flag, int *lobs, float *rho,
		float *t, float *cl, float *ct, float *ql, float *qt,
		complex ***response1, complex ***response2, 
		complex ***response3, FILE *outfp)
/******************************************************************************
Iunput:
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
acoustic		array[nor] of flags for receiver type
flag			array[nor] of flags for source type
rho 			array[nlayers] of densities
t   			array[nlayers] of layer thicknesses

Output
response1   	array[nw][nx][nor] of computed pressure wavefield
response2   	array[nw][nx][nor] of computed radial wavefield
response3   	array[nw][nx][nor] of computed vertical wavefield
******************************************************************************/
{
	int iw,iz,ip,ix;		/* loop counters */
	int ijk1=0,ik1,ik2,il,jl;	/* auxiliary indices */
	int ijk=0,ijk2,jj;		/* more auxiliary indices */
	int flg1,lrec;				
	float bp=0.0;			/* smallest ray parameter */
	float fp;			/* final ray parameter */
	float dp;			/* ray parameter increment */
	int nblock;			/* number of blocks to process */
	int block_size;			/* size of current block */
	int left;			/* remainder */
	int *prv,*ibl;			/* scratch arrays */
	float x1,w;
	float wbeg=0.0,wfin=0.0;	/* aux variables for output */
	complex cdp,wpie;		/* complex cdp and frequency */

	/* complex auxiliary variables */
	complex divfac,g,h,e,f,mu,x,y,z,a,b,c,d,t1,tem,temr,sig,sigh,rho1,rho2;
	complex det,psq,gl1,gt1,gl2,gt2,gam1,gam2,at1,at2,al1,al2;
	complex rvrb111,rvrb112,rvrb121,rvrb122,rvrb211,rvrb212,rvrb221;
	complex ewigh11,ewigh22,rvrb222;
	complex wrn011,wrn012,wrn021,wrn022,ewrd11,ewrd12,ewrd21,ewrd22;
	complex ewd0211,ewd0212,ewd0221,ewd0222,ewtu211,ewtu212;
	complex ewtu221,ewtu222,rtb111,rtb112,rtb121,rtb122,rtb211;
	complex rtb212,rtb221,rtb222,rtb311,rtb312,rtb321,rtb322;
	complex func1,func2,func3,func4,func5,func6;
	complex tun0p11,tun0p12,tun0p21,tun0p22,rd0np11,rd0np12,rd0np21;
	complex rd0np22,td0np11,td0np12,td0np21,td0np22;
	complex rd11,rd12,rd21,rd22,td11,td12,td21,td22,tu11,tu12,tu21,tu22;
	complex ru11,ru12,ru21,ru22,ruf11,ruf12,ruf21,ruf22;
	complex fact11,fact12,vuzs1,vuzs2,vdzs1,vdzs2;

	/* complex scratch arrays */
	complex *rd0n11=NULL,*rd0n12=NULL,*rd0n21=NULL,*rd0n22=NULL,
		*td0n11=NULL,*td0n12=NULL;
	complex *td0n21=NULL,*td0n22=NULL,*tun011=NULL,*tun012=NULL,
		*tun021=NULL,*tun022=NULL;
	complex *run011=NULL,*run012=NULL,*run021=NULL,*run022=NULL,
		*tusf11=NULL,*tusf12=NULL;
	complex *tusf21=NULL,*tusf22=NULL,*rdfs11=NULL,*rdfs12=NULL,
		*rdfs21=NULL,*rdfs22=NULL;
	complex *tdfs11=NULL,*tdfs12=NULL,*tdfs21=NULL,*tdfs22=NULL,
		*rurf11=NULL,*rurf12=NULL;
	complex *rurf21=NULL,*rurf22=NULL,*turf11=NULL,*turf12=NULL,
		*turf21=NULL,*turf22=NULL;
	complex *rdfr11=NULL,*rdfr12=NULL,*rdfr21=NULL,*rdfr22=NULL,
		*tdfr11=NULL,*tdfr12=NULL;
	complex *tdfr21=NULL,*tdfr22=NULL,*rnr11=NULL,*rnr12=NULL,
		*rnr21=NULL,*rnr22=NULL;
	complex *rdsr11=NULL,*rdsr12=NULL,*rdsr21=NULL,*rdsr22=NULL;
	/* complex *r11=NULL,*r12=NULL,*r21=NULL,*r22=NULL; */
	complex *texp=NULL,*vos1=NULL,*vos2=NULL,*vos3=NULL,*ux=NULL,
		*uz=NULL,*up=NULL;
	complex *alpha=NULL,*betha=NULL,*sigmad1=NULL,*sigmad2=NULL,
		*sigmau1=NULL,*sigmau2=NULL;
	complex *al=NULL,*at=NULL,*gl=NULL,*gt=NULL,*gam=NULL,*p=NULL,
		*pp=NULL,*prs=NULL,*pwin=NULL;

	/* complex pointers */	
	complex *rusf11=NULL,*rusf12=NULL,*rusf21=NULL,*rusf22=NULL;
	complex *tusr11=NULL,*tusr12=NULL,*tusr21=NULL,*tusr22=NULL;
	complex *rusr11=NULL,*rusr12=NULL,*rusr21=NULL,*rusr22=NULL;
	complex *tdsr11=NULL,*tdsr12=NULL,*tdsr21=NULL,*tdsr22=NULL;
	complex *rdnr11=NULL,*rdnr12=NULL,*rdnr21=NULL,*rdnr22=NULL;

	/* allocate working space */
	ux=alloc1complex(BLK);uz=alloc1complex(BLK);
	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);tdfs11=alloc1complex(BLK);
	tdfs12=alloc1complex(BLK);tdfs21=alloc1complex(BLK);
	tdfs22=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);rnr11=alloc1complex(BLK);
	rnr12=alloc1complex(BLK);rnr21=alloc1complex(BLK);
	rnr22=alloc1complex(BLK); rdsr11=alloc1complex(BLK);
	rdsr12=alloc1complex(BLK);rdsr21=alloc1complex(BLK);
	rdsr22=alloc1complex(BLK);vos1=alloc1complex(BLK);
	rusf21=alloc1complex(BLK);rusf22=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 */
	rdnr11=alloc1complex(BLK);rdnr12=alloc1complex(BLK);
	rdnr21=alloc1complex(BLK);rdnr22=alloc1complex(BLK);
	prv=alloc1int(nor);
	prs=alloc1complex(nor);
	ibl=alloc1int(BLK);
	p=alloc1complex(BLK);
	pp=alloc1complex(BLK);
	pwin=alloc1complex(BLK);
	al=alloc1complex(nlayers);
	at=alloc1complex(nlayers);
	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);

	/* set 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);
				response2[iw][ix][iz]=cmplx(0.0,0.0);
				response3[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 */
		if (iw==0) wbeg=w;
		if (iw==nw-1) wfin=w;
		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);

		/* 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);
		}


		/**********************************************************************
		* 		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);


			/******************************************************************
			*			start computation of reflectivity series			  *
			******************************************************************/
			/* loop over reflecting layers */
			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) {				/* initialize, first layer only */
					if ((at1.r==0.0)&&(at1.i==0.0)) {
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							gl1=gl[ijk1];
							ewigh11=cexp(cmul(cmplx(0.0,1.0),cmul(wpie,cmul(gl1,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);
							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]=cneg(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 {
						for (ip=0; ip<block_size; ip++) {
							ijk1=ik1+ip;
							psq=pp[ip];
							gl1=gl[ijk1];
							gt1=gt[ijk1];
							gam1=gam[ijk1];
							mu=cmul(cdiv(cmplx(1.0,0.0),cmul(at1,at1)),rho1);
							rtb211=crmul(cmul(p[ip],cmul(mu,gl1)),2.);
							rtb212=gam1;
							rtb221=cneg(gam1);
							rtb222=crmul(cmul(p[ip],cmul(mu,gt1)),2.);
							det=csub(cmul(rtb211,rtb222),cmul(rtb212,rtb221));
							rtb111=cdiv(rtb222,det);
							rtb112=cdiv(rtb212,det);
							rtb121=cneg(cdiv(rtb221,det));
							rtb122=cdiv(rtb211,det);
							ruf11=cadd(cmul(rtb211,rtb222),cmul(rtb212,rtb221));
							ruf12=cmul(rtb222,crmul(rtb212,2.));
							ruf21=cmul(rtb211,crmul(rtb221,2.));
							ruf22=ruf11;
							ruf11=cdiv(ruf11,det);
							ruf12=cdiv(ruf12,det);
							ruf21=cdiv(ruf21,det);