unif2aniso.c

su 的源代码库

	vp = ealloc1float(nz);
	vs = ealloc1float(nz);
	rho = ealloc1float(nz);
	epsilon = ealloc1float(nz);
	delta = ealloc1float(nz);
	gamma = ealloc1float(nz);
	phi = ealloc1double(nz);
	x0 = ealloc1float(ninf+1);
	z0 = ealloc1float(ninf+1);
	vp00 = ealloc1float(ninf+1);
	vs00 = ealloc1float(ninf+1);
	rho00 = ealloc1float(ninf+1);
	eps00 = ealloc1float(ninf+1);
	delta00 = ealloc1float(ninf+1);
	gamma00 = ealloc1float(ninf+1);
	phi00 = ealloc1double(ninf+1);

	dvpdx = ealloc1float(ninf+1);
	dvpdz = ealloc1float(ninf+1);
	dvsdx = ealloc1float(ninf+1);
	dvsdz = ealloc1float(ninf+1);
	drdx = ealloc1float(ninf+1);
	drdz = ealloc1float(ninf+1);

	dedx = ealloc1float(ninf+1);
	dedz = ealloc1float(ninf+1);
	dddx = ealloc1float(ninf+1);
	dddz = ealloc1float(ninf+1);
	dgdx = ealloc1float(ninf+1);
	dgdz = ealloc1float(ninf+1);


	xint =  ealloc1float(npmax);
	zint =  ealloc1float(npmax);
	xs =  ealloc1float(nx);
	zs =  ealloc2float(nx,ninf+1);

	c11 = ealloc2float(nz,nx);
	c13 = ealloc2float(nz,nx);
	c33 = ealloc2float(nz,nx);
	c44 = ealloc2float(nz,nx);
	c66 = ealloc2float(nz,nx);
	c15 = ealloc2float(nz,nx);
	c35 = ealloc2float(nz,nx);
	c55 = ealloc2float(nz,nx);
	rhooutput = ealloc2float(nz,nx);

	/* Compute uniformly sampled xs */
	for(ix=0; ix<nx; ++ix)
		xs[ix] = fx+ix*dx; 
	
	/* Input picked interfaces and make interpolation velocity values  */
	for(i=0; i<=ninf; ++i) {
		j= -1;
		do{
			j +=1;
			if(j>=npmax) 
		err("The point number on the %ith interface exceeds %i!",
				i,npmax);
			scanf("%f%f\n", &xint[j], &zint[j]);
		} while( zint[j]>-9999);
			npt = j;

		/* if linear interpolation or only one input sample */
		if (method[0]=='l' || nx==1) {
		    		intlin(npt,xint,zint,zint[0],zint[npt-1],
				nx,xs,zs[i]);
		/* else, if monotonic interpolation */
   		} else if (method[0]=='m') {
				zind = (float (*)[4])ealloc1float(npt*4);
				cmonot(npt,xint,zint,zind);
				intcub(0,npt,xint,zind,nx,xs,zs[i]);

    		/* else, if Akima interpolation */
   		} else if (method[0]=='a') {
				zind = (float (*)[4])ealloc1float(npt*4);
				cakima(npt,xint,zint,zind);
				intcub(0,npt,xint,zind,nx,xs,zs[i]);

    		/* else, if cubic spline interpolation */
    		} else if (method[0]=='s') {
			   	zind = (float (*)[4])ealloc1float(npt*4);
			  	csplin(npt,xint,zint,zind);
			   	intcub(0,npt,xint,zind,nx,xs,zs[i]);

 		/* else, if unknown method specified */
 		} else {
  			err("%s is an unknown interpolation method!\n",method);
 		}

		for(ix=0; ix<nx; ++ix){
 			if(i>0 && zs[i][ix]<zs[i-1][ix])
			 err("the %ith interface is above the %ith one at position \n\t\t(%g,%g)! \n", i,i-1,xs[ix],zs[i][ix]);
		}

	}
	 	
	/* Input layer velocities, densities, anisotropy parameters and */
	/* their derivatives */
	if (!getparfloat("x0",x0))
 		for(i=0;i<=ninf;++i) x0[i] = 0.;
	if (!getparfloat("z0",z0))
 		for(i=0;i<=ninf;++i) z0[i] = 0.;
	if (!getparfloat("vp00",vp00))
 		for(i=0;i<=ninf;++i) vp00[i] = 1500.+ 500*i;
	if (!getparfloat("vs00",vs00))
 		for(i=0;i<=ninf;++i) vs00[i] = sqrt(3)*(1500.+ 500*i);
	if (!getparfloat("rho00",rho00))
 		for(i=0;i<=ninf;++i) rho00[i] = 1000.0 + 10*i;
	if (!getparfloat("dvpdx",dvpdx)) 
		for(i=0;i<=ninf;++i) dvpdx[i] = 0.;
	if (!getparfloat("dvpdz",dvpdz)) 
		for(i=0;i<=ninf;++i) dvpdz[i] = 0.;
	if (!getparfloat("dvsdx",dvsdx)) 
		for(i=0;i<=ninf;++i) dvsdx[i] = 0.;
	if (!getparfloat("dvsdz",dvsdz)) 
		for(i=0;i<=ninf;++i) dvsdz[i] = 0.;
	if (!getparfloat("drdx",drdx)) 
		for(i=0;i<=ninf;++i) drdx[i] = 0.;
	if (!getparfloat("drdz",drdz)) 
		for(i=0;i<=ninf;++i) drdz[i] = 0.;
	if (!getparfloat("dedx",dedx)) 
		for(i=0;i<=ninf;++i) dedx[i] = 0.;
	if (!getparfloat("dedz",dedz)) 
		for(i=0;i<=ninf;++i) dedz[i] = 0.;
	if (!getparfloat("dgdx",dedx)) 
		for(i=0;i<=ninf;++i) dgdx[i] = 0.;
	if (!getparfloat("dgdz",dgdz)) 
		for(i=0;i<=ninf;++i) dgdz[i] = 0.;
	if (!getpardouble("phi00",phi00)) 
		for(i=0;i<=ninf;++i) phi00[i] = 0.;

	/* compute linear parameters */
	for(ix=0; ix<nx; ++ix){ 
		j = 1;
		for(iz=0,z=fz; iz<nz; ++iz,z+=dz){
			if(z>=zs[ninf][ix]) {
				vp[iz] = vp00[ninf]
					  +(xs[ix]-x0[ninf])*dvpdx[ninf]
					  +(z-z0[ninf])*dvpdz[ninf];
				vs[iz] = vs00[ninf]
					  +(xs[ix]-x0[ninf])*dvsdx[ninf]
					  +(z-z0[ninf])*dvsdz[ninf];
				rho[iz] = rho00[ninf]
					  +(xs[ix]-x0[ninf])*drdx[ninf]
					  +(z-z0[ninf])*drdz[ninf];
				epsilon[iz] = eps00[ninf]
					  +(xs[ix]-x0[ninf])*dedx[ninf]
					  +(z-z0[ninf])*dedz[ninf];
				delta[iz] = delta00[ninf]
					  +(xs[ix]-x0[ninf])*dddx[ninf]
					  +(z-z0[ninf])*dddz[ninf];
				gamma[iz] = gamma00[ninf]
					  +(xs[ix]-x0[ninf])*dgdx[ninf]
					  +(z-z0[ninf])*dgdz[ninf];
				phi[iz]	= phi00[ninf];
			} else if(z<=zs[1][ix]) {
				vp[iz] = vp00[0]
					  +(xs[ix]-x0[0])*dvpdx[0]
					  +(z-z0[0])*dvpdz[0];
				vs[iz] = vs00[0]
					  +(xs[ix]-x0[0])*dvsdx[0]
					  +(z-z0[0])*dvsdz[0];
				rho[iz] = rho00[0]
					  +(xs[ix]-x0[0])*drdx[0]
					  +(z-z0[0])*drdz[0];
				epsilon[iz] = eps00[0]
					  +(xs[ix]-x0[0])*dedx[0]
					  +(z-z0[0])*dedz[0];
				delta[iz] = delta00[0]
					  +(xs[ix]-x0[0])*dddx[0]
					  +(z-z0[0])*dddz[0];
				gamma[iz] = gamma00[0]
					  +(xs[ix]-x0[0])*dgdx[0]
					  +(z-z0[0])*dgdz[0];
				phi[iz]	= phi00[0];
			} else {
					for(j=j; z>zs[j][ix]; ++j);
				j -= 1;
				vp[iz] = vp00[j]
					 +(xs[ix]-x0[j])*dvpdx[j]
					 +(z-z0[j])*dvpdz[j];
				vs[iz] = vs00[j]
					 +(xs[ix]-x0[j])*dvsdx[j]
					 +(z-z0[j])*dvsdz[j];
				rho[iz] = rho00[j]
					 +(xs[ix]-x0[j])*drdx[j]
					 +(z-z0[j])*drdz[j];
				epsilon[iz] = eps00[j]
					 +(xs[ix]-x0[j])*dedx[j]
					 +(z-z0[j])*dedz[j];
				delta[iz] = delta00[j]
					 +(xs[ix]-x0[j])*dddx[j]
					 +(z-z0[j])*dddz[j];
				gamma[iz] = gamma00[j]
					 +(xs[ix]-x0[j])*dgdx[j]
					 +(z-z0[j])*dgdz[j];
				phi[iz]	= phi00[j];
			}

		}
		
		/* loop over gridpoints and do calculations */
		if (paramtype==1) { /* Thomsen anisotropy parameters */
			float b=0.0,c=0.0;	/* temporary variables */ 

			for(iz=0; iz<nz; ++iz){
				/* ... calculate stiffnesses */
				c33[ix][iz] = vp[iz]*vp[iz]*rho[iz];
				c44[ix][iz] = vs[iz]*vs[iz]*rho[iz];
				c11[ix][iz] = c33[ix][iz]*(1+2*epsilon[iz]);
				c66[ix][iz] = c44[ix][iz]*(1+2*gamma[iz]);
				c = c33[ix][iz]*(2*c44[ix][iz]*(1+delta[iz])
					- c33[ix][iz]*(1+2*delta[iz]));
				b = 2*c44[ix][iz];
				c13[ix][iz] = (-b+sqrt(b*b-4*c))/2;
				c15[ix][iz] = 0.0;
				c35[ix][iz] = 0.0;
				c55[ix][iz] = c44[ix][iz];

				rhooutput[ix][iz] = rho[iz];

			}
		} else {  /* Sayers anisotropy parameters */
			for(iz=0; iz<nz; ++iz){

				/* ... calculate stiffnesses */
				c33[ix][iz] = vp[iz]*vp[iz]*rho[iz];
				c44[ix][iz] = vs[iz]*vs[iz]*rho[iz];
				c11[ix][iz] = c33[ix][iz]*(1+2*epsilon[iz]);
				c66[ix][iz] = c44[ix][iz]*(1+2*gamma[iz]);
				c13[ix][iz] = (delta[iz]+1)*c33[ix][iz]
							-2*c44[ix][iz];
				c15[ix][iz] = 0.0;
				c35[ix][iz] = 0.0;
				c55[ix][iz] = c44[ix][iz];
				rhooutput[ix][iz] = rho[iz];
			}
		}

		for(iz=0; iz<nz; ++iz) { /* compute rotation */

				spar.a1111 = (double) c11[ix][iz];
				spar.a1133 = (double) c13[ix][iz];
				spar.a3333 = (double) c33[ix][iz];
				spar.a2323 = (double) c44[ix][iz];
				spar.a1212 = (double) c66[ix][iz];
				spar.a1113 = (double) c15[ix][iz];
				spar.a3313 = (double) c35[ix][iz];
				spar.a1313 = (double) c55[ix][iz];

				rottens2D(&spar,phi[iz]);

				c11[ix][iz] = (float) spar.a1111;
				c13[ix][iz] = (float) spar.a1133;
				c33[ix][iz] = (float) spar.a3333;
				c44[ix][iz] = (float) spar.a2323;
				c66[ix][iz] = (float) spar.a1212;
				c15[ix][iz] = (float) spar.a1113;
				c35[ix][iz] = (float) spar.a3313;
				c55[ix][iz] = (float) spar.a1313;
		}
	
		
	}

	
	/* Write the output files to disk */
	efwrite(*c11,sizeof(float),nz*nx,c11fp);
	efwrite(*c13,sizeof(float),nz*nx,c13fp);
	efwrite(*c33,sizeof(float),nz*nx,c33fp);
	efwrite(*c44,sizeof(float),nz*nx,c44fp);
	efwrite(*c66,sizeof(float),nz*nx,c66fp);
	efwrite(*c15,sizeof(float),nz*nx,c15fp);
	efwrite(*c35,sizeof(float),nz*nx,c35fp);
	efwrite(*c55,sizeof(float),nz*nx,c55fp);
	efwrite(*rhooutput,sizeof(float),nz*nx,rhofp);

	if(verbose){
		  warn("Output file for c11 parameter plane : %s ",c11_file);
	  	  warn("Output file for c13 parameter plane : %s ",c13_file);
	          warn("Output file for c33 parameter plane : %s ",c33_file);
	          warn("Output file for c44 parameter plane : %s ",c44_file);
	          warn("Output file for c66 parameter plane : %s ",c66_file);
	          warn("Output file for c15 parameter plane : %s ",c15_file);
	          warn("Output file for c35 parameter plane : %s ",c35_file);
	          warn("Output file for c55 parameter plane : %s ",c55_file);
	}


	/* Free workspace */
	free1float(vp);
	free1float(vs);
	free1float(rho);
	free1float(epsilon);
	free1float(delta);
	free1float(gamma);
	free1float(dvpdx);
	free1float(dvpdz);
	free1float(dvsdx);
	free1float(dvsdz);
	free1float(dedx);
	free1float(dedz);
	free1float(dgdx);
	free1float(dgdz);
	free1float(dddx);
	free1float(dddz);
	free1float(drdx);
	free1float(drdz);
	free2float(c11);
	free2float(c13);
	free2float(c33);
	free2float(c44);
	free2float(c66);
	free2float(c15);
	free2float(c35);
	free2float(c55);

					
	return(CWP_Exit());
		
}