elaray.c

su 的源代码库

    } else if((x1-x2)*(x1-x2) < eps*eps && (x1-xs)*(x1-xs) < eps*eps){
             *x = (xs < xmax-eps ? xs + eps : xs - eps);
	     return 1;
    } else if((x1-x3)*(x1-x3) < eps*eps && (x1-xs)*(x1-xs) < eps*eps){
             *x = (xs < xmax-eps ? xs + eps : xs - eps);
	     return 1;

    } else if((z1-z2)*(z1-z2) < eps*eps && (z1-zs)*(z1-zs) < eps*eps){
             *z = (zs < zmax-eps ? zs + eps : zs - eps);
	     return 1;

    } else if((z1-z3)*(z1-z3) < eps*eps && (z1-zs)*(z1-zs) < eps*eps){
             *z = (zs < zmax-eps ? zs + eps : zs - eps);
	     return 1;

    } else if((z2-z3)*(z2-z3) < eps*eps && (z2-zs)*(z2-zs) < eps*eps){
             *z = (zs < zmax-eps ? zs + eps : zs - eps);
	     return 1;


    /* now the hard test */
    } else if((x1-x2)*(x1-x2) < eps*eps){ 
	     m12=FLT_MAX;
    } else if((x1-x3)*(x1-x3) < eps*eps){ 
	     m13=FLT_MAX;
    } else if((x2-x3)*(x2-x3) < eps*eps){ 
	     m23=FLT_MAX;
    } 

    if(m13!=FLT_MAX) m13=(z1-z3)/(x1-x3);
    if(m23!=FLT_MAX) m23=(z2-z3)/(x2-x3);
    if(m12!=FLT_MAX) m12=(z1-z2)/(x1-x2);
    b12=z1-m12*x1;
    b13=z1-m13*x1;
    b23=z2-m23*x2;


    /* source on edge? */
    if((zs-m12*xs-b12)*(zs-m12*xs-b12)<eps*eps ||
       (zs-m13*xs-b13)*(zs-m13*xs-b13)<eps*eps || 
       (zs-m23*xs-b23)*(zs-m23*xs-b23)<eps*eps){
             *z = (zs < zmax-eps ? zs + eps : zs - eps);
             return 1;

    /* source is not on edge */
    } else {
	return 0;
    }

}


void writeWaves (FILE *wavefp, int nt, int nray, RayStep *rs[])
/* for each ray, write nt z(t),x(t) pairs uniformly sampled in time t */
{
	int iray,it;
	float tmax,dt,t1,t2,x1,z1,x,z,xstart,zstart;
	float vgx,vgz,ti,ddt;
	RayStep *rsi,*rs1,*rs2;
	
	/* determine maximum time for all rays */
	for (iray=0,tmax=0.0; iray<nray; ++iray)
		for (rsi=rs[iray]; rsi!=NULL; rsi=rsi->rsNext)
			tmax = MAX(tmax,rsi->t);
	
	/* determine time sampling interval */
	dt = tmax/(nt>1?nt-1:1);

	fprintf(stderr,"wavefront time increment = %g\n",dt);

	/* loop over rays */
	for (iray=0; iray<nray; ++iray) {
	
		/* initialize */
		rs1 = rs[iray];
		rs2 = rs1->rsNext;
		t1 = rs1->t;
		t2 = rs2->t;
		x1 = rs1->x;
		z1 = rs1->z;
		vgx = rs1->vgx;
		vgz = rs1->vgz;
	
		/* remember start x,z */
		xstart = x1;
		zstart = z1;
		
		/* loop over times */
		for (it=0,ti=0.0; it<nt; ++it,ti+=dt) {
			
			/* if necessary, go to next ray step */
			if (t2<ti) {

				do {
					rs1 = rs2;
					rs2 = rs1->rsNext;
					if (rs2==NULL) break;
				} while (rs2->t<ti);

				if (rs2==NULL) break;
				t1 = rs1->t;
				t2 = rs2->t;
				x1 = rs1->x;
				z1 = rs1->z;
				vgx = rs1->vgx;
				vgz = rs1->vgz;
			}
			
			ddt = ti-t1;
	
			/* compute x and z */
			x = x1+ddt*vgx;
			z = z1+ddt*vgz;
		
			/* write x and z */
			fwrite(&z,sizeof(float),1,wavefp);
			fwrite(&x,sizeof(float),1,wavefp);
		}
		
		/* finish writing x and z */
		while (it<nt) {
			fwrite(&z,sizeof(float),1,wavefp);
			fwrite(&x,sizeof(float),1,wavefp);
			++it;
		}
	}
}

void writeWaves2 (FILE *wavefp, float tw, int nray, RayStep *rs[])
/* for each ray, write z(tw),x(tw) pairs uniformly sampled in time t */
{
	int iray;
	float tmax,t1,t2,x1,z1,x,z,xstart,zstart;
	float vgx,vgz,ddt;
	RayStep *rsi,*rs1,*rs2;

	fprintf(stderr,"\n wavefront plotted for tw=%g\n",tw);

	/* determine maximum time for all rays */
	for (iray=0,tmax=0.0; iray<nray; ++iray)
		for (rsi=rs[iray]; rsi!=NULL; rsi=rsi->rsNext)
			tmax = MAX(tmax,rsi->t);
	
	/* loop over rays */
	for (iray=0; iray<nray; ++iray) {
	
		/* initialize */
		rs1 = rs[iray];
		rs2 = rs1->rsNext;
		t1 = rs1->t;
		t2 = rs2->t;
		x1 = rs1->x;
		z1 = rs1->z;
		vgx = rs1->vgx;
		vgz = rs1->vgz;
	
		/* remember start x,z */
		xstart = x1;
		zstart = z1;
		
		/* if necessary, go to next ray step */
		if (t2<tw) {

				do {
					rs1 = rs2;
					rs2 = rs1->rsNext;
					if (rs2==NULL) break;
				} while (rs2->t<tw);
				if (rs2==NULL){
					/* write xstart and zstart */
					fwrite(&zstart,sizeof(float),1,wavefp);
					fwrite(&xstart,sizeof(float),1,wavefp);
					 continue;
				}
				t1 = rs1->t;
				t2 = rs2->t;
				x1 = rs1->x;
				z1 = rs1->z;
				vgx = rs1->vgx;
				vgz = rs1->vgz;
		}
			
			ddt = tw-t1;
	
			/* compute x and z */
			x = x1+ddt*vgx;
			z = z1+ddt*vgz;
			/* write x and z */
			fwrite(&z,sizeof(float),1,wavefp);
			fwrite(&x,sizeof(float),1,wavefp);
	}
		
}
	
void writeRays (Model *m,FILE *rayfp, int nxz, int nray, RayStep *rs[],RayEnd re[],int krecord,int prim,FILE *ifp,FILE *outparfp)

/* for each ray, write x,z pairs */
{
	int iray,nxzw,icount;
	float x,z;
	RayStep *rsi;

	/* initialize counting */
        icount=0;

	/* loop over rays */
	for (iray=0; iray<nray; ++iray) {

	   if(( krecord == INT_MAX) || (krecord ==re[iray].kend &&
		 ( prim == INT_MAX || prim == re[iray].nref))){ 
                
 		/* count # of rays */
                icount +=1;
			
		/* loop through ray steps */
		for (rsi=rs[iray],nxzw=0; rsi!=NULL; rsi=rsi->rsNext,++nxzw){
			if (nxzw>=nxz) break;
			x = rsi->x;
			z = rsi->z;
			fwrite(&z,sizeof(float),1,rayfp);
			fwrite(&x,sizeof(float),1,rayfp);
		}

 		/* if necessary, repeat last (x,z) */
		while (nxzw<nxz) {
			fwrite(&z,sizeof(float),1,rayfp);
			fwrite(&x,sizeof(float),1,rayfp);
			++nxzw;
		}

	    } else {
                /* manipulate rayend information */
                re[iray].kend=-1;

            } /* closes if statement */

       } /* end loop over rays */

      if (ifp!=NULL)
   	fprintf(ifp,"\n Number of rays written into output : %i \n",icount);

      fprintf(outparfp,"\n %i\n",icount);

}

static void gvelreal (float a1111, float a3333, float a1133, float a1313,
	float a1113, float a3313, float px, float pz, float *vgx,  float *vgz,
	float *g11n, float *g13n, float *g33n)
/*****************************************************************************
compute group velocity and polarization
******************************************************************************
Input:
aijkl		stiffness elements
px,pz		slowness elements

Output:
*vgx, *vgz	group velocities
*g11,*g13,*g33  polarizations
******************************************************************************/
{
	float gamm11,gamm13,gamm33;
	float px2,pz2,pxz,den,g11,g13,g33;

	px2   = px*px;
	pz2   = pz*pz;
	pxz   = px*pz;


	/*anisotropy parameters*/
	gamm11 = a1111*px2+ a1313*pz2 +2*a1113*pxz;
	gamm33 = a3333*pz2 + a1313*px2+2*a3313*pxz;
	gamm13 = (a1133+a1313)*pxz+ a1113*px2+ a3313*pz2;
	den     = 1/(gamm11+gamm33-2);
	g11     = (gamm33-1)*den;
	g33     = (gamm11-1)*den;
	g13     = -gamm13*den;

	/* computing ray (group) velocities */
	*vgx =  (a1111*px*g11+(a1133+a1313)*pz*g13+a3313*pz*g33+
		a1113*(pz*g11+2*px*g13)+a1313*g33*px);
	*vgz =  (a3333*pz*g33+(a1133+a1313)*px*g13+a1113*px*g11+
		+a3313*(px*g33+2*pz*g13)+a1313*g11*pz);

	/* kill round-off */
	if( g11 < -10*FLT_EPSILON || g33 < -10*FLT_EPSILON) 
		err("\n ERROR: g11 or g33 <0 in initialization \n");

	else if( g11 < 0.0 )
		g11 = 0;

	else if( g33 < 0.0 )
		g33 = 0;


	*g11n=g11;
	*g13n=g13;
	*g33n=g33;


}


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;
}

int findnewmode (int mode , int* newmode, int conv, int mindex)
/*****************************************************************************
find the next ray mode
******************************************************************************
Input:
mode		incidence mode
conv		mode conversion=1
mindex		medium index

Output:
*newmode	pointer to new mode
******************************************************************************
Notes:
routine returns 1 if new mode found
******************************************************************************
Credits:  Andreas Rueger, Colorado School of Mines, 02/06/94
******************************************************************************/
{

	/* P isotropic */
	if( (mode == 0 && mindex == 0 && conv == 0) ||
	    (mode == 3 && mindex == 0 && conv == 0) ||
	    (mode == 4 && mindex == 0 && conv == 1) ||
	    (mode == 1 && mindex == 0 && conv == 1))
		*newmode = 0;

	/* SV isotropic */
	else if( (mode == 1 && mindex == 0 && conv == 0) ||
	    (mode == 4 && mindex == 0 && conv == 0) ||
	    (mode == 0 && mindex == 0 && conv == 1) ||
	    (mode == 3 && mindex == 0 && conv == 1))
		*newmode = 1;

	/* SH isotropic */
	else if( (mode == 2 && mindex == 0 && conv == 0) ||
	    (mode == 2 && mindex == 0 && conv == 1) ||
	    (mode == 5 && mindex == 0 && conv == 0) ||
	    (mode == 5 && mindex == 0 && conv == 1))
		*newmode = 2;

	/* qP anisotropic */
	else if( (mode == 0 && mindex > 0 && conv == 0) ||
	    (mode == 1 && mindex > 0 && conv == 1) ||
	    (mode == 3 && mindex > 0 && conv == 0) ||
	    (mode == 4 && mindex > 0 && conv == 1))
		*newmode = 3;

	/* qSV anisotropic */
	else if( (mode == 0 && mindex > 0 && conv == 1) ||
	    (mode == 1 && mindex > 0 && conv == 0) ||
	    (mode == 3 && mindex > 0 && conv == 1) ||
	    (mode == 4 && mindex > 0 && conv == 0))
		*newmode = 4;

	/* qSH anisotropic */
	else if( (mode == 2 && mindex > 0 && conv == 1) ||
	    (mode == 2 && mindex > 0 && conv == 0) ||
	    (mode == 5 && mindex > 0 && conv == 1) ||
	    (mode == 5 && mindex > 0 && conv == 0))
		*newmode = 5;
	else 
		return -1;

	return 1;
}
int findqPqSV(float s, float c, float pl, float a1111, float a3333,
	float a1133,float a1313, float a1113, float a3313, int mode, float
	*pxnew, float *pznew, float *vgx, float *vgz, float *g11, float *g13,  
	float *g33, float rt, FILE *ifp)
/*****************************************************************************
Continue slowness across interface 
******************************************************************************
Input:
s,c		slope of interface measured with respect to horizontal
pl		tangential component of slowness
aijkl		density-normalized stiffness elements
mode		Ray mode
rt		=-1 transmission
		=1 reflection
******************************************************************************
Output:
*pxnew		address pointing to new px
*pznew 		address pointing to new pz
-1		no root found
1		root found
*vgx,*vgz  	group velocity pointers
*g11,*g33,*g13  polarizations

******************************************************************************
Author:  Andreas Rueger, Colorado School of Mines, 02/01/94
******************************************************************************/
{
	int check1,check2;
	float vgm,vgxd,vgzd,pxnewd,pznewd,plold;
 
	check1 = check2 = 0;
	plold = pl;

	if(mode == 3){

		check1=solveForSlowqP(s,c,pl,a1111,a3333,a1133,a1313,a1113,
		a3313,pxnew,pznew,rt);
		
		
	}else if (mode ==4)

		check1=solveForSlowqSV(s,c,pl,a1111,a3333,a1133,a1313,a1113,
		a3313,pxnew,pznew,rt);

	else 
		err(" wrong mode in findqP \n ");


	if (check1 == 1){
		gvelreal (a1111,a3333,a1133,a1313,a1113,a3313,*pxnew,
				*pznew,vgx,vgz,g11,g13,g33);
	
		vgxd=*vgx;
		vgzd=*vgz;

		vgm=-s*vgxd+c*vgzd;

		if(rt == 1 && vgm > 0 ){
			check2 = 1;
		} else if( rt == -1 && vgm <0 ){
			check2 = 1;
		} else
			check2 = 0;

		if (ifp!=NULL && check2 == 1)
			fprintf(ifp,"\n first try in rootfinder "
				"was successfull \n");
	}
	if (check1 != 1 || check2 != 1){

		if (ifp!=NULL)
		  fprintf(ifp,"\n need second try in rootfinder \n");

		if(mode == 3)
		  check1=solveForSlowqP(s,c,pl,a1111,a3333,a1133,a1313,a1113,
		  a3313,pxnew,pznew,-1*rt);
		
		
		else if (mode ==4)
		  check1=solveForSlowqSV(s,c,pl,a1111,a3333,a1133,a1313,a1113,
		  a3313,pxnew,pznew,-1*rt);

	

		
		
		if(check1 != 1){
		  
	                fprintf(ifp," rootfinder not sucessful !!! \n");
			return 0;
		}
		
	        else {
		
			gvelreal (a1111,a3333,a1133,a1313,a1113,a3313,*pxnew,
				*pznew,vgx,vgz,g11,g13,g33);
	
			vgxd=*vgx;
			vgzd=*vgz;

			vgm=-s*vgxd+c*vgzd;

			if(rt == 1 && vgm < 0){
				return 0;
			} else if(rt == -1 && vgm >0){
				return 0;
			} else {
				check2 = 1;
				if (ifp!=NULL)
				fprintf(ifp,"\n second try "
				"in rootfinder was successfull \n");
			}
		}
	}

	pxnewd=*pxnew;
	pznewd=*pznew;

	pl = c*pxnewd+pznewd*s;
	if((pl-plold)*(pl-plold)<0.000000001 && check2 == 1 && check1 == 1)
		return 1;
	else
		return 0;		

}