fdmod.c

地震波正演和显示模块

                        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 cmonot (int n, float x[], float y[], float yd[][4])
/*****************************************************************************
compute cubic interpolation coefficients via the Fritsch-Carlson method,
which preserves monotonicity
******************************************************************************
Input:
n               number of samples
x               array[n] of monotonically increasing or decreasing abscissae
y               array[n] of ordinates

Output:
yd              array[n][4] of cubic interpolation coefficients (see notes)
******************************************************************************
Notes:
The computed cubic spline coefficients are as follows:
yd[i][0] = y(x[i])    (the value of y at x = x[i])
yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])
To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
use the computed coefficients as follows:
y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))

The Fritsch-Carlson method yields continuous 1st derivatives, but 2nd
and 3rd derivatives are discontinuous.  The method will yield a
monotonic interpolant for monotonic data.  1st derivatives are set
to zero wherever first divided differences change sign.

For more information, see Fritsch, F. N., and Carlson, R. E., 1980,
Monotone piecewise cubic interpolation:  SIAM J. Numer. Anal., v. 17,
n. 2, p. 238-246.

Also, see the book by Kahaner, D., Moler, C., and Nash, S., 1989,
Numerical Methods and Software, Prentice Hall.  This function was
derived from SUBROUTINE PCHEZ contained on the diskette that comes
with the book.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 09/30/89
Modified:  Dave Hale, Colorado School of Mines, 02/28/91
        changed to work for n=1.
Modified:  Dave Hale, Colorado School of Mines, 08/04/91
        fixed bug in computation of left end derivative
*****************************************************************************/
{
        int i;
        float h1,h2,del1,del2,dmin,dmax,hsum,hsum3,w1,w2,drat1,drat2,divdf3;

        /* copy ordinates into output array */
        for (i=0; i<n; i++)
                yd[i][0] = y[i];

        /* if n=1, then use constant interpolation */
        if (n==1) {
                yd[0][1] = 0.0;
                yd[0][2] = 0.0;
                yd[0][3] = 0.0;
                return;

        /* else, if n=2, then use linear interpolation */
        } else if (n==2) {
                yd[0][1] = yd[1][1] = (y[1]-y[0])/(x[1]-x[0]);
                yd[0][2] = yd[1][2] = 0.0;
                yd[0][3] = yd[1][3] = 0.0;
                return;
        }
 /* set left end derivative via shape-preserving 3-point formula */
        h1 = x[1]-x[0];
        h2 = x[2]-x[1];
        hsum = h1+h2;
        del1 = (y[1]-y[0])/h1;
        del2 = (y[2]-y[1])/h2;
        w1 = (h1+hsum)/hsum;
        w2 = -h1/hsum;
        yd[0][1] = w1*del1+w2*del2;
        if (yd[0][1]*del1<=0.0)
                yd[0][1] = 0.0;
        else if (del1*del2<0.0) {
                dmax = 3.0*del1;
                if (ABS(yd[0][1])>ABS(dmax)) yd[0][1] = dmax;
        }

 /* loop over interior points */
        for (i=1; i<n-1; i++) {

                /* compute intervals and slopes */
                h1 = x[i]-x[i-1];
                h2 = x[i+1]-x[i];
                hsum = h1+h2;
                del1 = (y[i]-y[i-1])/h1;
                del2 = (y[i+1]-y[i])/h2;

                /* if not strictly monotonic, zero derivative */
                if (del1*del2<=0.0) {
                        yd[i][1] = 0.0;

                /*
                 * else, if strictly monotonic, use Butland's formula:
                 *      3*(h1+h2)*del1*del2
                 * -------------------------------
 * ((2*h1+h2)*del1+(h1+2*h2)*del2)
                 * computed as follows to avoid roundoff error
                 */
                } else {
                        hsum3 = hsum+hsum+hsum;
                        w1 = (hsum+h1)/hsum3;
                        w2 = (hsum+h2)/hsum3;
                        dmin = MIN(ABS(del1),ABS(del2));
                        dmax = MAX(ABS(del1),ABS(del2));
                        drat1 = del1/dmax;
                        drat2 = del2/dmax;
                        yd[i][1] = dmin/(w1*drat1+w2*drat2);
                }
        }
 /* set right end derivative via shape-preserving 3-point formula */
        w1 = -h2/hsum;
        w2 = (h2+hsum)/hsum;
        yd[n-1][1] = w1*del1+w2*del2;
        if (yd[n-1][1]*del2<=0.0)
                yd[n-1][1] = 0.0;
        else if (del1*del2<0.0) {
                dmax = 3.0*del2;
                if (ABS(yd[n-1][1])>ABS(dmax)) yd[n-1][1] = dmax;
        }
 /* compute 2nd and 3rd derivatives of cubic polynomials */
        for (i=0; i<n-1; i++) {
                h2 = x[i+1]-x[i];
                del2 = (y[i+1]-y[i])/h2;
                divdf3 = yd[i][1]+yd[i+1][1]-2.0*del2;
                yd[i][2] = 2.0*(del2-yd[i][1]-divdf3)/h2;
                yd[i][3] = (divdf3/h2)*(6.0/h2);
        }
        yd[n-1][2] = yd[n-2][2]+(x[n-1]-x[n-2])*yd[n-2][3];
        yd[n-1][3] = yd[n-2][3];
}
/*4444444444444444444444444444444444444444444444*/


void strchop(char *s, char *t)
/***********************************************************************
strchop - chop off the tail end of a string "s" after a "," returning
          the front part of "s" as "t".
************************************************************************
Notes:
Based on strcpy in Kernighan and Ritchie's C [ANSI C] book, p. 106.
************************************************************************
Author: CWP: John Stockwell and Jack K. Cohen, July 1995
***********************************************************************/
{

        while ( (*s != ',') && (*s != '\0') ) {
                 *t++ = *s++;
        }
        *t='\0';
}


/*****************************************************************************
XINDEX - determine index of x with respect to an array of x values

xindex          determine index of x with respect to an array of x values

******************************************************************************
Input:
nx              number of x values in array ax
ax              array[nx] of monotonically increasing or decreasing x values
x               the value for which index is to be determined
index           index determined previously (used to begin search)

Output:
index           for monotonically increasing ax values, the largest index
                for which ax[index]<=x, except index=0 if ax[0]>x;
                for monotonically decreasing ax values, the largest index
                for which ax[index]>=x, except index=0 if ax[0]<x
******************************************************************************
Notes:
This function is designed to be particularly efficient when called
repeatedly for slightly changing x values; in such cases, the index
returned from one call should be used in the next.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 12/25/89
*****************************************************************************/
/**************** end self doc ********************************/
void
xindex (int nx, float ax[], float x, int *index)
/*****************************************************************************
determine index of x with respect to an array of x values
******************************************************************************Input:
nx              number of x values in array ax
ax              array[nx] of monotonically increasing or decreasing x values
x               the value for which index is to be determined
index           index determined previously (used to begin search)

Output:
index           for monotonically increasing ax values, the largest index
                for which ax[index]<=x, except index=0 if ax[0]>x;
                for monotonically decreasing ax values, the largest index
                for which ax[index]>=x, except index=0 if ax[0]<x
******************************************************************************
Notes:
This function is designed to be particularly efficient when called
repeatedly for slightly changing x values; in such cases, the index
returned from one call should be used in the next.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 12/25/89
*****************************************************************************/
{
        int lower,upper,middle,step;

        /* initialize lower and upper indices and step */
        lower = *index;
        if (lower<0) lower = 0;
        if (lower>=nx) lower = nx-1;
        upper = lower+1;
        step = 1;
  /* if x values increasing */
        if (ax[nx-1]>ax[0]) {

                /* find indices such that ax[lower] <= x < ax[upper] */
                while (lower>0 && ax[lower]>x) {
                        upper = lower;
                        lower -= step;
                        step += step;
                }
                if (lower<0) lower = 0;
                while (upper<nx && ax[upper]<=x) {
                        lower = upper;
                        upper += step;
                        step += step;
                }
                if (upper>nx) upper = nx;
   /* find index via bisection */
                while ((middle=(lower+upper)>>1)!=lower) {
                        if (x>=ax[middle])
                                lower = middle;
                        else
                                upper = middle;
                }

        /* else, if not increasing */
        } else {
 /* find indices such that ax[lower] >= x > ax[upper] */
                while (lower>0 && ax[lower]<x) {
                        upper = lower;
                        lower -= step;
                        step += step;
                }
                if (lower<0) lower = 0;
                while (upper<nx && ax[upper]>=x) {
                        lower = upper;
                        upper += step;
                        step += step;
                }
                if (upper>nx) upper = nx;
  /* find index via bisection */
                while ((middle=(lower+upper)>>1)!=lower) {
                        if (x<=ax[middle])
                                lower = middle;
                        else
                                upper = middle;
                }
        }

        /* return lower index */
        *index = lower;
}


void warn(char *fmt, ...)
{
        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 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);
}

void *alloc1 (size_t n1, size_t size)
{
	void *p;

	if ((p=malloc(n1*size))==NULL)
		return NULL;
	return p;
}
void **alloc2 (size_t n1, size_t n2, size_t size)
{
	size_t i2;
	void **p;

	if ((p=(void**)malloc(n2*sizeof(void*)))==NULL) 
		return NULL;
	if ((p[0]=(void*)malloc(n2*n1*size))==NULL) {
		free(p);
		return NULL;
	}
	for (i2=0; i2<n2; i2++)
		p[i2] = (char*)p[0]+size*n1*i2;
	return p;
}

float *alloc1float(size_t n1) 
{
	return (float*)alloc1(n1,sizeof(float)); 
} 

float **alloc2float(size_t n1, size_t n2)
{
	return (float**)alloc2(n1,n2,sizeof(float));
}

void free1 (void *p) 
{ 
	free(p);
}

void free2 (void **p)
{
	free(p[0]);
	free(p);
}

void free1float(float *p) 
{ 
	free1(p);
} 

int *alloc1int(size_t n1)
{ 
	return (int*)alloc1(n1,sizeof(int));
}

void free1int(int *p) 
{ 
	free1(p); 
}

void free2float(float **p)
{
	free2((void**)p);
}