unif.c

地震波正演和显示模块

#include <math.h>
#include <stdarg.h>
#include "../lib/getpars.c"
#include "../lib/atopkge.c"
#define O2 0.5000000
#define O6 0.1666667
#define EXIT_SUCCESS 0

/*********************** self documentation ******************************/
char *sdoc[] = {
" 									 ",
" UNIF2 - generate a 2-D UNIFormly sampled velocity profile from a layered",
"  	 model. In each layer, velocity is a linear function of position.",
" 									",
"  unif2 < infile > outfile [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 v00 is specified		",
" z0=0.0,0.0,..., 	depth z at which v00 is specified		",
" v00=1500,2000,2500...,	velocity at each x0,z0 (m/sec)		",
" dvdx=0.0		derivative of velocity with distance x (dv/dx)	",
" dvdz=0.0		derivative of velocity with depth z (dv/dz)	",
" 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 are 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",
" 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.	",
"									",
" Example using test input file generating feature:			",
" unif2 tfile=testfilename    produces a 5 interface demonstration model",
" unif2 < testfilename | psimage n1=100 n2=100 d1=10 d2=10 | ...	",
"									",
NULL};
/**************** end self doc *******************************************/

char *sample[] = {
"   0        0",
" 500        0",
"   1      -99999",
"   0       50",
" 500       50",
"   1      -99999",
"   0      200",
"  50      200",
" 200      280",
" 275      280",
" 500      180",
"   1      -99999",
"   0      400",
" 50       400",
" 200      375",
" 250      390",
" 275      390",
" 500      250",
"   1      -99999",
"   0      475",
" 50       475",
" 200      450",
" 250      400",
" 275      400.",
" 375      460",
" 500      465",
"   1      -99999",
"   0      500",
" 500      500",
"   1      -99999",
NULL};
/* pointer to that test dataset */
char **sampleptr=sample;

void warn(char *fmt, ...);
void err(char *fmt, ...);
void requestdoc(int flag);
void pagedoc(void);
void intcub (int ideriv, int nin, float xin[], float ydin[][4],int nout, float xout[], float yout[]);
void intlin (int nin, float xin[], float yin[], float yinl, float yinr,int nout, float xout[], float yout[]); 
void cmonot (int n, float x[], float y[], float yd[][4]);
void cakima (int n, float x[], float y[], float yd[][4]);
void strchop(char *s, char *t);
void xindex (int nx, float ax[], float x, int *index);
void csplin (int n, float x[], float y[], float yd[][4]);
void *alloc1 (size_t n1, size_t size);
void **alloc2 (size_t n1, size_t n2, size_t size);
float *alloc1float(size_t n1);
float **alloc2float(size_t n1, size_t n2);
void free1 (void *p);
void free2 (void **p);
void free1float(float *p);
int *alloc1int(size_t n1);
void free1int(int *p);
void free2float(float **p);
int xargc;  char **xargv;


int
main(int argc, char **argv)
{    
	int nx, nz, i,j,ix,iz,ninf,npmax,npt;
	float z,fx,fz,dx,dz;
	float *v,*x0,*z0,*v00,*dvdx,*dvdz,**zs,*xs,*xint,*zint,(*zind)[4];
   	char *method="linear";
	FILE *out_file=stdout;
	char *tfile="";
	FILE *tfilefp;
	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(0);

	/* if the user specifies the name of a test datafile */
	/* output test dataset to this file */
	getparstring("tfile", &tfile);
      	
	
	if (*tfile!='\0') {
		if((tfilefp=fopen(tfile,"w"))==NULL)
			err("cannot open tfile=%s\n",tfile);

		while(*sampleptr) fprintf(tfilefp,"%s\n", *sampleptr++);
		fclose(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);
    	
    	warn("ninf=%d\n",ninf);
    	warn("nx=%d\n",nx);
    	warn("nz=%d\n",nz);
    	warn("dx=%f\n",dx);
    	warn("dz=%f\n",dz);
    	warn("fx=%f\n",fx);
    	warn("fz=%f\n",fz);
    	warn("npmax=%d\n",npmax);
	warn("method=%s\n",method);
	/* allocate space */
        v = alloc1float(nz);
        x0 = alloc1float(ninf+1);
        z0 = alloc1float(ninf+1);
        v00 = alloc1float(ninf+1);
        dvdx = alloc1float(ninf+1);
        dvdz = alloc1float(ninf+1);
        xint =  alloc1float(npmax);
        zint =  alloc1float(npmax);
        xs =  alloc1float(nx);
        zs =  alloc2float(nx,ninf+1);

	/* 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);exit(0);}
			scanf("%f%f", &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])alloc1float(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])alloc1float(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])alloc1float(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 and velocity 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("v00",v00))
 		for(i=0;i<=ninf;++i) v00[i] = 1500.+ 500*i;
	if (!getparfloat("dvdx",dvdx)) 
		for(i=0;i<=ninf;++i) dvdx[i] = 0.;
	if (!getparfloat("dvdz",dvdz)) 
		for(i=0;i<=ninf;++i) dvdz[i] = 0.;

	/* compute linear velocity */
	for(ix=0; ix<nx; ++ix){ 
		j = 1;
		for(iz=0,z=fz; iz<nz; ++iz,z+=dz){
			if(z>=zs[ninf][ix]) 
				v[iz] = v00[ninf]+(xs[ix]-x0[ninf])*dvdx[ninf]
					+(z-z0[ninf])*dvdz[ninf];
			else if(z<=zs[1][ix]) 
				v[iz]=v00[0]+(xs[ix]-x0[0])*dvdx[0]
					+(z-z0[0])*dvdz[0];
			else{
					for(j=j; z>zs[j][ix]; ++j);
				j -= 1;
				v[iz] =	v00[j]+(xs[ix]-x0[j])*dvdx[j]
					+(z-z0[j])*dvdz[j];
			}
		}
		fwrite(v,sizeof(float),nz,out_file);
	}
	
				
	return(EXIT_SUCCESS);
		
}

void warn(char *fmt, ...)
/***************************************************************************
print warning on application program error
****************************************************************************/
{
	va_list args;

	if (EOF == fflush(stdout)) {
		fprintf(stderr, "\nwarn: fflush failed on stdout");
	}
	fprintf(stderr, "\n%s: ", xargv[0]);
	va_start(args,fmt);
	vfprintf(stderr, fmt, args);
	va_end(args);
	fprintf(stderr, "\n");
	return;
}


void err(char *fmt, ...)
/***************************************************************************
print warning on application program error and die
****************************************************************************/
{
	va_list args;

 
	if (EOF == fflush(stdout)) {
		fprintf(stderr, "\nerr: fflush failed on stdout");
	}
	fprintf(stderr, "\n%s: ", xargv[0]);
	va_start(args,fmt);
	vfprintf(stderr, fmt, args);
	va_end(args);
	fprintf(stderr, "\n");
	exit(EXIT_FAILURE);
}


void requestdoc(int flag)
/*************************************************************************** 
print selfdocumentation as directed by the user-specified flag
**************************************************************************** 

The flag argument distinguishes these cases:
            flag = 0; fully defaulted, no stdin
            flag = 1; usual case
            flag = n > 1; no stdin and n extra args required

pagedoc:
Intended to be called by pagedoc(), but conceivably could be
used directly as in:
      if (xargc != 3) selfdoc();*/


{
        switch(flag) {
        case 1:
                if (xargc == 1 && isatty(STDIN)) pagedoc();
        break;
        case 0:
                if (xargc == 1 && isatty(STDIN) && isatty(STDOUT)) pagedoc();
        break;
        default:
                if (xargc <= flag) pagedoc();
        break;
        }
        return;
}

void pagedoc(void)
{
        extern char *sdoc[];
		char **p = sdoc;
        FILE *fp;
		char *pager;
		char cmd[32];

		if ((pager=getenv("PAGER")) != (char *)NULL)
			sprintf(cmd,"%s 1>&2", pager);
		else 
			sprintf(cmd,"more 1>&2");


        fflush(stdout);
       /*  fp = popen("more -22 1>&2", "w"); */
       /*  fp = popen("more  1>&2", "w"); */
        fp = popen(cmd, "w");
	while(*p) (void)fprintf(fp, "%s\n", *p++);
        pclose(fp);

        exit(EXIT_FAILURE);
}


/* sample input data set */

void intcub (int ideriv, int nin, float xin[], float ydin[][4], 
	int nout, float xout[], float yout[])
/*****************************************************************************
evaluate y(x), y'(x), y''(x), ... via piecewise cubic interpolation
******************************************************************************
Input:
ideriv		=0 if y(x) desired; =1 if y'(x) desired, ...
nin		length of xin and ydin arrays
xin		array[nin] of monotonically increasing or decreasing x values
ydin		array[nin][4] of y(x), y'(x), y''(x), and y'''(x)
nout		length of xout and yout arrays
xout		array[nout] of x values at which to evaluate y(x), y'(x), ...

Output:
yout		array[nout] of y(x), y'(x), ... values
******************************************************************************
Notes:
xin values must be monotonically increasing or decreasing.

Extrapolation of the function y(x) for xout values outside the range
spanned by the xin values is performed using the derivatives in 
ydin[0][0:3] or ydin[nin-1][0:3], depending on whether xout is closest
to xin[0] or xin[nin-1], respectively.
*/
{
	static int idx;
	int iout;
	float delx;

	/* y(x) is desired, then */
	if (ideriv==0) {
		for (iout=0; iout<nout; iout++) {
			xindex(nin,xin,xout[iout],&idx);
			delx = xout[iout]-xin[idx];
			yout[iout] = (ydin[idx][0]+delx*
				(ydin[idx][1]+delx*
				(ydin[idx][2]*O2+delx*
				(ydin[idx][3]*O6))));
		}
		
	/* else, if y'(x) is desired, then */
	} else if (ideriv==1) {
		for (iout=0; iout<nout; iout++) {
			xindex(nin,xin,xout[iout],&idx);
			delx = xout[iout]-xin[idx];
			yout[iout] = (ydin[idx][1]+delx*
				(ydin[idx][2]+delx*
				(ydin[idx][3]*O2)));
		}
		
	/* else, if y''(x) is desired, then */
	} else if (ideriv==2) {
		for (iout=0; iout<nout; iout++) {
			xindex(nin,xin,xout[iout],&idx);
			delx = xout[iout]-xin[idx];
			yout[iout] = (ydin[idx][2]+delx*
				(ydin[idx][3]));
		}
		
	/* else, if y'''(x) is desired, then */
	} else if (ideriv==3) {
		for (iout=0; iout<nout; iout++) {
			xindex(nin,xin,xout[iout],&idx);
			delx = xout[iout]-xin[idx];
			yout[iout] = (ydin[idx][3]);
		}
		
	/* else, if any other derivative is desired, then */
	} else {
		for (iout=0; iout<nout; iout++)
			yout[iout] = 0.0;
	}
}

void
intlin (int nin, float xin[], float yin[], float yinl, float yinr, 
	int nout, float xout[], float yout[])
/*****************************************************************************
evaluate y(x) via linear interpolation of y(x[0]), y(x[1]), ...
******************************************************************************
Input:
nin		length of xin and yin arrays
xin		array[nin] of monotonically increasing or decreasing x values
yin		array[nin] of input y(x) values
yinl		value used to extraplate y(x) to left of input yin values
yinr		value used to extraplate y(x) to right of input yin values
nout		length of xout and yout arrays
xout		array[nout] of x values at which to evaluate y(x)

Output:
yout		array[nout] of linearly interpolated y(x) values
******************************************************************************
Notes:
xin values must be monotonically increasing or decreasing.

Extrapolation of the function y(x) for xout values outside the range
spanned by the xin values in performed as follows:

	For monotonically increasing xin values,
		yout=yinl if xout<xin[0], and yout=yinr if xout>xin[nin-1].

	For monotonically decreasing xin values, 
		yout=yinl if xout>xin[0], and yout=yinr if xout<xin[nin-1].

If nin==1, then the monotonically increasing case is used.
*/
{
	static int idx;
	int jout;
	float x;

	/* if input x values are monotonically increasing, then */
	if(xin[0]<=xin[nin-1]) 
	{
		for (jout=0; jout<nout; jout++) {
			x = xout[jout];
			if (x<xin[0])
				yout[jout] = yinl;
			else if (x>xin[nin-1])
				yout[jout] = yinr;
			else if (x==xin[nin-1] || nin==1)
				yout[jout] = yin[nin-1];
			else {   
				xindex(nin,xin,x,&idx);
				yout[jout] = yin[idx]+(x-xin[idx])
					*(yin[idx+1]-yin[idx])
					/(xin[idx+1]-xin[idx]);
			}
		}
	
	/* else, if input x values are monotonically decreasing, then */
	} else {
		for (jout=0; jout<nout; jout++) {
			x = xout[jout];
			if (x>xin[0])
				yout[jout] = yinl;
			else if (x<xin[nin-1])
				yout[jout] = yinr;
			else if (x==xin[nin-1] || nin==1)
				yout[jout] = yin[nin-1];
			else {
				xindex(nin,xin,x,&idx);
				yout[jout] = yin[idx]+(x-xin[idx])
					*(yin[idx+1]-yin[idx])
					/(xin[idx+1]-xin[idx]);
			}
		}
	}
}