unif2aniso.c

su 的源代码库

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

/* UNIF2ANISO: $Revision: 1.8 $ ; $Date: 2005/10/25 17:59:27 $	*/

#include "par.h"
#include "anisotropy.h"

/*********************** self documentation ******************************/
char *sdoc[] = {
" 									 ",
" UNIF2ANISO - generate a 2-D UNIFormly sampled profile of elastic	",
"	constants from a layered model.					",
" 									",
"  unif2aniso < infile [Parameters]					",
" 									",
" Required Parameters:							",
" none 									",
" 									",
" Optional Parameters:							",
" ninf=5	number of interfaces					",
" nx=100	number of x samples (2nd dimension)			",
" nz=100	number of z samples (1st dimension)			",
" dx=10		x sampling interval					",
" dz=10		z sampling interval					",
" 									",
" npmax=201	maximum number of points on interfaces			",
" 									",
" fx=0.0	first x sample						",
" fz=0.0	first z sample						",
" 									",
" 									",
" x0=0.0,0.0,..., 	distance x at which vp00 and vs00 are specified	",
" z0=0.0,0.0,..., 	depth z at which vp00 and vs00 are specified	",
" 									",
" vp00=1500,2000,...,	P-velocity at each x0,z0 (m/sec)		",
" vs00=866,1155...,	S-velocity at each x0,z0 (m/sec)		",
" rho00=1000,1100,...,	density at each x0,z0 (kg/m^3)			",
" 									",
" eps00=0,0,0...,	Thomsen or Sayers epsilon			",
" delta00=0,0,0...,	Thomsen or Sayers delta				",
" gamma00=0,0,0...,	Thomsen or Sayers gamma				",
" 									",
" drdx=0.0,0.0,...,	x-derivative of density (d rho/dx)		",
" drdz=0.0,0.0,...,	z-derivative of density (d rho/dz)		",
" 									",
" dvpdx=0.0,0.0,...,	x-derivative of P-velocity (dvp/dx)		",
" dvpdz=0.0,0.0,...,	z-derivative of P-velocity (dvs/dz)		",
" 									",
" dvsdx=0.0,0.0,...,	x-derivative of S-velocity (dvs/dx)		",
" dvsdz=0.0,0.0,...,	z-derivative of S-velocity (dvs/dz)		",
" 									",
" dedx=0.0,0.0,...,	x-derivative of epsilon (de/dx)			",
" dedz=0.0,0.0,...,	z-derivative of epsilon with depth z (de/dz)	",
" 									",
" dddx=0.0,0.0,...,	x-derivative of delta (dd/dx)			",
" dddz=0.0,0.0,...,	z-derivative of delta (dd/dz)			",
" 									",
" dgdz=0.0,0.0,...,	x-derivative of gamma (dg/dz)			",
" dgdx=0.0,0.0,...,	z-derivative of gamma (dg/dx)			",
" 									",
" phi00=0,0,..., 	rotation angle(s) in each layer			",
" 									",
" ...output filenames 							",
" c11_file=c11_file	output filename for c11 values	 		",
" c13_file=c13_file	output filename for c13 values	 		",
" c15_file=c15_file	output filename for c15 values	 		",
" c33_file=c33_file	output filename for c33 values	 		",
" c35_file=c35_file	output filename for c35 values	 		",
" c44_file=c44_file	output filename for c44 values	 		",
" c55_file=c55_file	output filename for c55 values	 		",
" c66_file=c66_file	output filename for c66 values	 		",
" 									",
" rho_file=rho_file	output filename for density values 		",
" 									",
" paramtype=1   =1 Thomsen parameters, =0 Sayers parameters(see below)	",
" method=linear		for linear interpolation of interface		",
" 			=mono for monotonic cubic interpolation of interface",
"			=akima for Akima's cubic interpolation of interface",
"			=spline for cubic spline interpolation of interface",
" 									",
" tfile=		=testfilename  if set, a sample input dataset is",
" 			 output to \"testfilename\".			",
" 			 						",
" Notes:								",
" The input file is an ASCII file containing x z values representing a	",
" piecewise continuous velocity model with a flat surface on top.	",
" 			 						",
" The surface and each successive boundary between media is represented ",
" by a list of selected x z pairs written column form. The first and	",
" last x values must be the same for all boundaries. Use the entry	",
" 1.0  -99999  to separate the entries for successive boundaries. No	",
" boundary may cross another. Note that the choice of the method of	",
" interpolation may cause boundaries to cross that do not appear to	",
" cross in the input data file.						",
" 			 						",
" The number of interfaces is specified by the parameter \"ninf\". This ",
" number does not include the top surface of the model. The input data	",
" format is the same as a CSHOT model file with all comments removed.	",
" 			 						",
" The algorithm works by transforming the P-wavespeed , S-wavespeed,	",
" density and the Thomsen or Sayers parameters epsilon, delta, and gamma",
" into elastic stiffness coefficients. Furthermore, the	user can specify",
" rotations, phi, to the elasticity tensor in each layer.		",
"  									",
" Common ranges of Thomsen parameters are				",
"  epsilon:  0.0 -> 0.5							",
"  delta:   -0.2 -> 0.4							",
"  gamma:    0.0 -> 0.4							",
" 									",
" If only P-wave, S-wave velocities and density is given as input,	",
" the model is, by definition,  isotropic.				",
" 									",
" If files containing Thomsen/Sayers parameters are given, the model	",
" will be assumed to have VTI symmetry.		 			",
"									",
" Example using test input file generating feature:			",
" unif2aniso tfile=testfilename  produces a 5 interface demonstration model",
" unif2aniso < testfilename 						",
" ximage < c11_file n1=100 n2=100					",
" ximage < c13_file n1=100 n2=100					",
" ximage < c15_file n1=100 n2=100					",
" ximage < c33_file n1=100 n2=100					",
" ximage < c35_file n1=100 n2=100					",
" ximage < c44_file n1=100 n2=100					",
" ximage < c55_file n1=100 n2=100					",
" ximage < c66_file n1=100 n2=100					",
" ximage < rho_file n1=100 n2=100					",
"									",
NULL};

/*
 * Credits:
 *      CWP: John Stockwell, April 2005. 
 * 	CWP: based on program unif2 by Zhenyue Liu, 1994 
 */

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

/* sample input data set */
char *sample[] = {
"   0         0",
"1000         0",
"   1    -99999",
"   0       200",
"1000       200",
"   1    -99999",
"   0       400",
" 100.      400",
" 400       560",
" 550.      560",
"1000.      360",
"   1    -99999",
"   0       800",
" 100.      800",
" 400.      750.",
" 500.      780",
" 550.      780",
"1000.      500",
"   1.   -99999",
"   0       950",
" 100       950",
" 400       900.",
" 500       800.",
" 550       800.",
" 750       920",
"1000       930",
"   1    -99999",
"   0      1000",
"1000      1000",
"   1    -99999",
NULL};

/* pointer to that test dataset */
char **sampleptr=sample;

Stiff2D spar;		/* stiffness tensor			*/

int
main(int argc, char **argv)
{
	int nx;		/* size of 2nd (slow) dimension of data */
	int nz;		/* size of 1st (fast) dimension of data */
	int i,j,ix,iz;	/* counters				*/
	int ninf;	/* number of interfaces			*/
	int npmax;	/* max number of points on an interface	*/
	int npt;	/* number of points on an interface	*/

	int paramtype=0;/* =0 Thomsen, =1 Sayers 		*/
	int verbose=0;	/* =0 silent, =1 chatty 		*/

	float z;	/* fast dimension coordinate		*/
	float fx;	/* first value in x			*/
	float fz;	/* first value in z			*/
	float dx;	/* x sampling interval			*/
	float dz;	/* z sampling interval				*/


	float *vp=NULL;		/* P-velocity function calculated	*/
	float *vs=NULL;		/* S-velocity function calculated	*/
	float *rho=NULL;	/* density function calculated		*/
	float *epsilon=NULL;	/* epsilon function  calculated		*/
	float *delta=NULL;	/* delta function calculated		*/
	float *gamma=NULL;	/* delta function calculated		*/
	double *phi=NULL;	/* angle to rotate elasticity tensor */


	float *x0=NULL;		/* initial x position values		*/
	float *z0=NULL;		/* initial z position values		*/
	float *vp00=NULL;	/* P-wavespeeds at x0,z0 positions	*/
	float *vs00=NULL;	/* S-wavespeeds at x0,z0 positions	*/
	float *rho00=NULL;	/* densities at x0,z0 positions		*/
	float *eps00=NULL;	/* epsilon at x0,z0 positions		*/
	float *delta00=NULL;	/* delta at x0,z0 positions		*/
	float *gamma00=NULL;	/* gamma at x0,z0 positions		*/
	double *phi00=NULL;	/* angle to rotate elasticity tensor */

	float *dvpdx=NULL;	/* x rates of change of V-P by layer	*/ 
	float *dvpdz=NULL;	/* z rates of change of V-P by layer	*/
	float *dvsdx=NULL;	/* x rates of change of V-S by layer	*/ 
	float *dvsdz=NULL;	/* z rates of change of V-S by layer	*/
	float *drdx=NULL;	/* x rates of change of density by layer*/ 
	float *drdz=NULL;	/* z rates of change of density by layer*/
	float *dedx=NULL;	/* x rates of change of epsilon by layer*/ 
	float *dedz=NULL;	/* z rates of change of epsilon by layer*/
	float *dddx=NULL;	/* x rates of change of delta by layer*/ 
	float *dddz=NULL;	/* z rates of change of delta by layer*/
	float *dgdx=NULL;	/* x rates of change of gamma by layer*/ 
	float *dgdz=NULL;	/* z rates of change of gamma by layer	*/


	float **zs=NULL;	/* array of z values for interpolator	*/
	float *xs=NULL;		/* array of x values for interpolator	*/
	float *xint=NULL;	/* x values defining an interface	*/
	float *zint=NULL;	/* z values defining an interface	*/

	float (*zind)[4]=NULL;	/* array of interpolation coefficients	*/ 

   	char *method="linear";	/* interpolation method		*/
	char *tfile="";		/* test filename		*/
	FILE *tfilefp;		/* ...its file pointer		*/	

	/* .. stiffness values */
	float **c11=NULL;	/* stiffness c11		*/
	float **c13=NULL;	/* stiffness c13		*/
	float **c33=NULL;	/* stiffness c33		*/
	float **c44=NULL;	/* stiffness c44		*/
	float **c66=NULL;	/* stiffness c66		*/
	float **c15=NULL;	/* stiffness c15		*/
	float **c35=NULL;	/* stiffness c15		*/
	float **c55=NULL;	/* stiffness c15		*/

	float **rhooutput=NULL;	/* density values		*/

	/* ... output files */
	char *c11_file;		/* c11 filename		*/
	FILE *c11fp=NULL;	/* ... its file pointer */
	char *c13_file;		/* c13 filename		*/
	FILE *c13fp=NULL;	/* ... its file pointer */
	char *c33_file;		/* c33 filename		*/
	FILE *c33fp=NULL;	/* ... its file pointer */
	char *c44_file;		/* c44 filename		*/
	FILE *c44fp=NULL;	/* ... its file pointer */
	char *c66_file;		/* c66 filename		*/
	FILE *c66fp=NULL;	/* ... its file pointer */
	char *c15_file;		/* c15 filename		*/
	FILE *c15fp=NULL;	/* ... its file pointer */
	char *c35_file;		/* c35 filename		*/
	FILE *c35fp=NULL;	/* ... its file pointer */
	char *c55_file;		/* c55 filename		*/
	FILE *c55fp=NULL;	/* ... its file pointer */
	char *rho_file;		/* rho filename		*/
	FILE *rhofp=NULL;	/* ... its file pointer */

	/* Hook up getpar */
	initargs(argc, argv);
	requestdoc(1);
	
	/* If the user specifies the name of a test datafile */
	/* output test dataset to this file */
	getparstring("tfile", &tfile);
	if (*tfile!='\0') {
		if((tfilefp=efopen(tfile,"w"))==NULL)
			err("cannot open tfile=%s\n",tfile);

		while(*sampleptr) fprintf(tfilefp,"%s\n", *sampleptr++);
		efclose(tfilefp);
		exit(EXIT_SUCCESS);
	}

	/* Get parameters */
	if (!getparint("ninf",&ninf))	ninf=5;
	if (!getparint("nx",&nx))	nx=100;
	if (!getparint("nz",&nz))	nz=100;
	if (!getparfloat("dx",&dx))	dx=10.0;
	if (!getparfloat("dz",&dz))	dz=10.0;
	if (!getparfloat("fx",&fx)) fx = 0.;
	if (!getparfloat("fz",&fz)) fz = 0.;
	if (!getparint("npmax",&npmax)) npmax=201;
    	getparstring("method",&method);
	if (!getparint("paramtype",&paramtype))	paramtype=0;
	if (!getparint("verbose",&paramtype))	verbose=0;

	/* Open output files */
	if( !getparstring("c11_file",&c11_file)) c11_file = "c11_file";
	c11fp = efopen(c11_file,"w");
	if( !getparstring("c13_file",&c13_file)) c13_file = "c13_file";
	c13fp = efopen(c13_file,"w");
	if( !getparstring("c33_file",&c33_file)) c33_file = "c33_file";
	c33fp = efopen(c33_file,"w");
	if( !getparstring("c44_file",&c44_file)) c44_file = "c44_file";
	c44fp = efopen(c44_file,"w");
	if( !getparstring("c66_file",&c66_file)) c66_file = "c66_file";
	c66fp = efopen(c66_file,"w");
	if( !getparstring("c15_file",&c15_file)) c15_file = "c15_file";
	c15fp = efopen(c15_file,"w");
	if( !getparstring("c35_file",&c35_file)) c35_file = "c35_file";
	c35fp = efopen(c35_file,"w");
	if( !getparstring("c55_file",&c55_file)) c55_file = "c55_file";
	c55fp = efopen(c55_file,"w");
	if( !getparstring("rho_file",&rho_file)) rho_file = "rho_file";
	rhofp = efopen(rho_file,"w");
	
	/* Allocate space */