elasyn.c

su 的源代码库

{
	int iti,itf,it,ix,itmin,itmax,temp,sub;
	float frac,x1,x2,odg,tempf,aa1,aa2;


	itmin = 0;
	itmax = 0;
	odg   = 1/dg;

	for(it=0; it<ng; ++it)
		im[it]=0;

	x1  = x[0]-fg;
	iti = x1*odg;
	aa1 = a[0];

	for(ix=1; ix<na; ++ix){
		x2  = x[ix]-fg;
		itf = x2*odg;
		aa2 = a[ix];
		sub = 1;

		if(itf<iti){
			temp = itf;
			itf  = iti;
			iti  = temp;
			tempf= x1;
			x1   = x2;
			x2   = tempf;
			tempf= aa1;
			aa1  = aa2;
			aa2  = tempf;
			sub  = 0;
		}

		if((x2-x1)>tol){
			for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
				ag[it][im[it]] = 0.0;
				++im[it];
			};

		} else{
			for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
				itmin = MIN(itmin,it);
				itmax = MAX(itmax,it);
				frac  = (it*dg-x1)/(x2-x1);
				ag[it][im[it]]= (1.0-frac)*aa1+frac*aa2;
				++im[it];
			};
		}
		
		if(sub) {
			iti = itf;
			x1  = x2;
			aa1 = aa2;
		}
	}

	for(it=0; it<itmin; ++it){
		ag[it][0] = 0.0;
		++im[it];
	}
	for(it=itmax+1; it<ng; ++it){
		ag[it][0] = 0.0;
		++im[it];
	}
}


static void makesyn (int reftrans, int ng, float *xg, int compon,
        float *zg, float fpeak,	int nt, float dt, float ft, 
	FILE *ifp, float **axg, float **azg, float **phg, 
	float **tg, float dangle, int *im, FILE *xfp, FILE *zfp)
{

	int ntfft,nw,iw,ig,it;
	float dw,w,wpeak,delay,*syn,cosf,sinf;
	double campr,campi,cosw,sinw,cosd,sind,tempd;
	complex *cwave,**csynx,**csynz;

	/* constants */
	ntfft = npfaro(nt,2*nt);
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	wpeak = 2*PI*fpeak;
	delay = 0.0;

	/* allocate workspace */
	syn = ealloc1float(ntfft);
	cwave = ealloc1complex(nw);
	csynx = ealloc2complex(nw,ng);
	csynz = ealloc2complex(nw,ng);

	/* initialize synthetics */
	for (ig=0; ig<ng; ++ig){
		for (iw=0; iw<nw; ++iw){
			csynx[ig][iw] = cmplx(0.0,0.0);
			csynz[ig][iw] = cmplx(0.0,0.0);
		}
	}

	/********** make complex ricker wavelet ************/
	for (iw=0,w=0.0; iw<nw; ++iw,w+=dw) {
		cwave[iw] = cricker(w,wpeak,delay);
		cwave[iw] = crmul(cwave[iw],sqrt(w));
	}

	
	/* loop over receiver locations */
	for (ig=0; ig<ng; ++ig) {
         /* fprintf(stderr,"ig=%i im=%i tg=%f a=%f \n",
		  ig,im[ig],tg[ig][0],axg[ig][0]); */
	  
          /* account for multiple arrivals */
	  for(it=0; it<im[ig]; ++it){
	    cosf  = cos(phg[ig][it]);
	    sinf  = sin(phg[ig][it]);

	    /* compute cos, sin, and exp increments */
	    cosd = cos(dw*(tg[ig][it]-ft));
	    sind = sin(dw*(tg[ig][it]-ft));

	    /* initialize cos, sin, and exp */
	    cosw = 1.0;
	    sinw = 0.0;

 	    /* for all frequencies */
	    for(iw=0, w=0.0; iw<nw; ++iw, w+=dw){
              
		campr = cosf*cosw-sinf*sinw;
		campi = cosf*sinw+sinf*cosw;

		csynx[ig][iw].r += campr*axg[ig][it];
		csynx[ig][iw].i += campi*axg[ig][it];

		csynz[ig][iw].r += campr*azg[ig][it];
		csynz[ig][iw].i += campi*azg[ig][it];

		/* update cos, sin */
		tempd = cosw*cosd-sinw*sind;
		sinw = cosw*sind+sinw*cosd;
		cosw = tempd;
	    }

	  }

	  
	  /* x component seismograms */
	  if(compon == 1 || compon == 0){

          	/********** apply wavelet to synthetics *************/
	   	for (iw=0; iw<nw; ++iw)
		  csynx[ig][iw] = cmul(csynx[ig][iw],cwave[iw]);

 	   	/********** inverse Fourier transform ***************/
	   	pfacr(-1,ntfft,csynx[ig],syn);

	  	fwrite(syn,sizeof(float),nt,xfp);
           }

	  /* z component seismograms */
	  if(compon ==3 || compon == 0){

          	/********** apply wavelet to synthetics *************/
	   	for (iw=0; iw<nw; ++iw)
		  csynz[ig][iw] = cmul(csynz[ig][iw],cwave[iw]);

 	   	/********** inverse Fourier transform ***************/
	   	pfacr(-1,ntfft,csynz[ig],syn);

	  	fwrite(syn,sizeof(float),nt,zfp);
          }

     }
        
}	



void rayJacob(float dg, float *x, int *num, float *pz, float *ph,
	float *ax, float *az, int nrays)
{
	int iray,diffnum,caustic;
	float jac1,jac2,jac;

	caustic=0;

	/****** first ray *******/
	jac = (x[0]-x[1])*pz[0];
	diffnum = num[1]-num[0];
	if(diffnum !=1) jac=jac/diffnum;

	ax[0] /= sqrt(ABS(jac));
	az[0] /= sqrt(ABS(jac));
	ph[0]=0;

	/****** all but the last ray **************/
	for(iray=1; iray<nrays-1;iray++){
		jac1 = (x[iray]-x[iray+1])*pz[iray];
		diffnum = num[iray+1]-num[iray];
		if(diffnum != 1)
			jac1 /= diffnum;

		jac2 = (x[iray-1]-x[iray])*pz[iray];
		diffnum = num[iray]-num[iray-1];
		if(diffnum != 1)
			jac2 /= diffnum;

		/* detect caustics 
		if(jac2*jac1 < 0 && caustic==0)
			caustic=1;
		else if(jac2*jac1 <0 && caustic==1)
			caustic=0;*/
		
		/* cusps */
		if(jac1*jac1 < 0.001*dg*dg && jac1*jac1 < 0.001*dg*dg)
			jac=jac;
		else if(jac1*jac1 < 0.001*dg*dg)
			jac=jac2;
		else if(jac2*jac2 < 0.001*dg*dg)
			jac=jac1;
		else
			jac=(jac1+jac2)/2.0;

		ax[iray] /= sqrt(ABS(jac));
		az[iray] /= sqrt(ABS(jac));

		if(caustic==1) ph[iray] -= PI/2.0;

	}

	/****** last ray *******/
	jac = (x[nrays-1]-x[nrays-2])*pz[nrays-1];
	diffnum = num[1]-num[0];
	if(diffnum !=1) jac=jac/diffnum;


	ax[nrays-1] /= sqrt(ABS(jac));
	az[nrays-1] /= sqrt(ABS(jac));
}


void interpara(int na, float *a, float *x, int ng, float dg, float fg,
		float tol, int *im, float **ag)
{
	int iti,itf,it,ix,itmin,itmax,temp,sub,ixl;
	float frac,x1,x2,odg;
	float a1=0.0,a2=0.0,b1=0.0,b2=0.0,c1=0.0,c2=0.0,t1,t2,xl,xl2;
	float x0=0.0,x3,u0=0.0,u1,u2,u3,dxx,dxx3,eps,epsx;
	double den3,num1,den2,den1,num2,num3;


	itmin = 0;
	itmax = 0;
	odg   = 1/dg;
	eps   = 0.001;
	epsx  = 0.01;

	for(it=0; it<ng; ++it)
		im[it]=0;

	x1  = x[0]-fg;
	iti = x1*odg;
	u1 = a[0];
	ixl=-1;

	for(ix=1; ix<na; ++ix){
		x2  = x[ix]-fg;
		itf = x2*odg;
		u2 = a[ix];
		sub = 1;
		x3=x[ix+1]-fg;
		u3=a[ix+1];

		if(itf<iti){
			temp = itf;
			itf  = iti;
			iti  = temp;
			sub  = 0;
		}

		if(ABS(x2-x1)>tol){
			for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
				ag[it][im[it]] = 0.0;
				im[it] += 1;
			};

		} else{
			if(ix==1 || ix==na-1){
				dxx=x2-x1;
				if(ABS(dxx)<eps){
					for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
						itmin = MIN(itmin,it);
						itmax = MAX(itmax,it);
						ag[it][im[it]]= 0.5*(u1+u2);
						im[it] += 1;
					}
				} else{
					
					for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
						itmin = MIN(itmin,it);
						itmax = MAX(itmax,it);
						frac  = (it*dg-x1)/dxx;
						ag[it][im[it]]= (1.0-frac)*u1+frac*u2;
						im[it] += 1;
					}
				}
			} else{
				dxx=x2-x1;
				dxx3=dxx*(x3-x2)*(x1-x0);
				if(ABS(dxx)<eps){
					for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
						itmin = MIN(itmin,it);
						itmax = MAX(itmax,it);
						ag[it][im[it]]= 0.5*(u1+u2);
						im[it] += 1;
					}
				} else if(ABS(dxx3)<epsx){
					
					for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
						itmin = MIN(itmin,it);
						itmax = MAX(itmax,it);
						frac  = (it*dg-x1)/dxx;
						ag[it][im[it]]= (1.0-frac)*u1+frac*u2;
						im[it] += 1;
					}
				} else{
					if(ixl!=ix-1){
						num1=(x1-x0)*(x2-x0);
						den1=1/num1;
						num2=(x1-x0)*(x2-x1);
						den2=1/num2;
						num3=(x2-x1)*(x2-x0);
						den3=1/num3;
						a1=u0*x1*x2*den1-u1*x2*x0*den2+
							u2*x0*x1*den3;
						b1=-u0*(x1+x2)*den1+u1*(x2+x0)*den2-
							u2*(x1+x0)*den3;
						c1=u0*den1-u1*den2+u2*den3;
					}
					num1=(x2-x1)*(x3-x1);
					den1=1/num1;
					num2=(x2-x1)*(x3-x2);
					den2=1/num2;
					num3=(x3-x2)*(x3-x1);
					den3=1/num3;
					a2=u1*x2*x3*den1-u2*x3*x1*den2+
						u3*x1*x2*den3;
					b2=-u1*(x2+x3)*den1+u2*(x3+x1)*den2-
						u3*(x2+x1)*den3;
					c2=u1*den1-u2*den2+u3*den3;
					for(it=MAX(iti+1,0); it<=MIN(itf,ng-1); ++it){
						itmin = MIN(itmin,it);
						itmax = MAX(itmax,it);
						xl=it*dg;
						xl2=xl*xl;
						t1= a1+b1*xl+c1*xl2;
						t2= a2+b2*xl+c2*xl2;
						ag[it][im[it]]=.5*(t1+t2);
						im[it] += 1;
					}
					a1=a2;
					b1=b2;
					c1=c2;
					ixl=ix;
				}
			}
		}
		
		if(sub)
			iti = itf;
		u0=u1;
		x0=x1;
		x1=x2;
		u1=u2;
	}

	for(it=0; it<itmin; ++it){
		ag[it][0] = 0.0;
		im[it] += 1;
	}
	for(it=itmax+1; it<ng; ++it){
		ag[it][0] = 0.0;
		im[it] += 1;
       	}
}
void polar(float px, float pz, float g11, float g13,
	float g33, float *polx, float *polz, int mode )
/*****************************************************************************
compute polarization oriented along slowness. This convention is identical to that used in Aki&Richards, pages 148ff.                     
*******************************************************************************
Input:
px,pz		slowness components
g11,g13,g33 	polarizations squared
mode=0,1,3,4	ray-mode(qP.qSV)

Output:
polx,polz	polarization components

Author:  Andreas Rueger, Colorado School of Mines, 03/14/94
******************************************************************************/
{
	float polarx,polarz;

	if(g11 < FLT_EPSILON){
		polarx=0;
		polarz=sqrt(g33);
	} else if(g33 < FLT_EPSILON) {
		polarz=0;
		polarx=sqrt(g11);
	} else {
		polarx=sqrt(g11)*SGN(g13);
		polarz=sqrt(g33);
        }
	
	if((mode==0 || mode==3) && px*polarx+pz*polarz <0 ){
		polarx=-polarx;
		polarz=-polarz;
	} else if((mode==4 || mode==1) && px*polarx <0){
		polarx=-polarx;
		polarz=-polarz;
	}


	*polx=polarx;
	*polz=polarz;
}