modeling.c

su 的源代码库

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

/* MODELING: $Revision: 1.4 $ ; $Date: 1996/09/06 18:06:51 $	*/

/*********************** self documentation **********************/
/*************************************************************************
MODELING - Seismic Modeling Subroutines for SUSYNLV and clones

decodeReflectors	parse reflectors parameter string
decodeReflector		decode a particular reflector
breakReflectors		break up reflectors duplicating interior (x,z) points
makeref			make piecewise cubic reflectors		
raylv2		 	Trace ray between two points,
				for linear velocity v = v00+dvdx*x+dvdz*z

addsinc		Add sinc wavelet to trace at specified time and
			with specified amplitude

makericker	make a Ricker wavelet

**************************************************************************
Function Prototypes:
void decodeReflectors (int *nrPtr,
	float **aPtr, int **nxzPtr, float ***xPtr, float ***zPtr);
int decodeReflector (char *string,
	float *aPtr, int *nxzPtr, float **xPtr, float **zPtr);
void breakReflectors (int *nr, float **ar, 
	int **nu, float ***xu, float ***zu);
void makeref (float dsmax, int nr, float *ar, 
	int *nu, float **xu, float **zu, Reflector **r);
void raylv2 (float v00, float dvdx, float dvdz,
	float x0, float z0, float x, float z,
	float *c, float *s, float *t, float *q);
void addsinc (float time, float amp,
	int nt, float dt, float ft, float *trace);
void makericker (float fpeak, float dt, Wavelet **w);
	
**************************************************************************
Notes:
Typedefs used by Hale's modeling
typedef struct ReflectorSegmentStruct {
	float x;	* x coordinate of segment midpoint *
	float z;	* z coordinate of segment midpoint *
	float s;	* x component of unit-normal-vector *
	float c;	* z component of unit-normal-vector *
} ReflectorSegment;
typedef struct ReflectorStruct {
	int ns;			* number of reflector segments *
	float ds;		* segment length *
	float a;		* amplitude of reflector *
	ReflectorSegment *rs;	* array[ns] of reflector segments *
} Reflector;
typedef struct WaveletStruct {
	int lw;			* length of wavelet *
	int iw;			* index of first wavelet sample *
	float *wv;		* wavelet sample values *
} Wavelet;

These are items used in SUSYNLV, SUSYNVXZ, SUSYNLVCW.

**************************************************************************
Author: Dave Hale, Colorado School of Mines, 09/17/91
**************************************************************************/
/**************** end self doc ********************************/

#include "par.h"

/*
 * Credits: CWP Dave Hale, 09/17/91,  Colorado School of Mines
 */


void decodeReflectors (int *nrPtr,
	float **aPtr, int **nxzPtr, float ***xPtr, float ***zPtr)
/*************************************************************************
decodeReflectors - parse reflectors parameter string
**************************************************************************
Output:
nrPtr		pointer to nr an int specifying number of reflectors
aPtr		pointer to a specifying reflector amplitudes
nxzPtr		pointer to nxz specifying number of (x,z) pairs defining the
		reflectors
xPtr		pointer to array[x][nr] of x values for the entire model
zPtr		array[z][nr] of z values for the entire model

***************************************************************************
Author: Dave Hale, Colorado School of Mines, 09/17/91
**************************************************************************/
{
	int nr,*nxz,ir;
	float *a,**x,**z;
	char t[1024],*s;

	/* count reflectors */
	nr = countparname("ref");
	if (nr==0) nr = 1;
	
	/* allocate space */
	a = ealloc1(nr,sizeof(float));
	nxz = ealloc1(nr,sizeof(int));
	x = ealloc1(nr,sizeof(float*));
	z = ealloc1(nr,sizeof(float*));

	/* get reflectors */
	for (ir=0; ir<nr; ++ir) {
		if (!getnparstring(ir+1,"ref",&s)) s = "1:1,2;4,2";
		strcpy(t,s);
		if (!decodeReflector(t,&a[ir],&nxz[ir],&x[ir],&z[ir]))
			err("Reflector number %d specified "
				"incorrectly!\n",ir+1);
	}

	/* set output parameters before returning */
	*nrPtr = nr;
	*aPtr = a;
	*nxzPtr = nxz;
	*xPtr = x;
	*zPtr = z;
}

/* parse one reflector specification; return 1 if valid, 0 otherwise */
int decodeReflector (char *string,
	float *aPtr, int *nxzPtr, float **xPtr, float **zPtr)
/**************************************************************************
decodeReflector - parse one reflector specification
***************************************************************************
Input:
string		string representing reflector

Output:
aPtr		pointer to amplitude value
nxzPtr		pointer to number of x,z pairs
xPtr		array of x values for one reflector
zPtr		array of z values for one reflector
***************************************************************************
Author: Dave Hale, Colorado School of Mines, 09/17/91
**************************************************************************/
{
	int nxz,ixz;
	float a,*x,*z;
	char *s,*t;

	/* if specified, get reflector amplitude; else, use default */
	s = string;
	if (strchr(s,':')==NULL) {
		a = 1.0;
		s = strtok(s,",;\0");
	} else {
		/* if (strcspn(s,":")>=strcspn(s,",;\0")) return 0;*/
		if (((unsigned int) strcspn(s,":")) >=
		    ((unsigned int) strcspn(s,",;\0"))) return 0;
		a = atof(s=strtok(s,":"));
		s = strtok(NULL,",;\0");
	}

	/* count x and z values, while splitting string into tokens */
	for (t=s,nxz=0; t!=NULL; ++nxz)
		t = strtok(NULL,",;\0");
	
	/* total number of values must be even */
	if (nxz%2) return 0;

	/* number of (x,z) pairs */
	nxz /= 2;

	/* 2 or more (x,z) pairs are required */
	if (nxz<2) return 0;

	/* allocate space */
	x = ealloc1(nxz,sizeof(float));
	z = ealloc1(nxz,sizeof(float));

	/* convert (x,z) values */
	for (ixz=0; ixz<nxz; ++ixz) {
		x[ixz] = atof(s);
		s += strlen(s)+1;
		z[ixz] = atof(s);
		s += strlen(s)+1;
	}

	/* set output parameters before returning */
	*aPtr = a;
	*nxzPtr = nxz;
	*xPtr = x;
	*zPtr = z;
	return 1;
}

/* Break up reflectors by duplicating interior (x,z) points */
void breakReflectors (int *nr, float **ar, 
	int **nu, float ***xu, float ***zu)
{
	int nri,nro,*nui,*nuo,ir,jr,iu;
	float *ari,*aro,**xui,**zui,**xuo,**zuo;

	/* input reflectors */
	nri = *nr;
	ari = *ar;
	nui = *nu;
	xui = *xu;
	zui = *zu;

	/* number of output reflectors */
	for (ir=0,nro=0; ir<nri; ++ir)
		nro += nui[ir]-1;

	/* make output reflectors and free space for input reflectors */
	aro = ealloc1float(nro);
	nuo = ealloc1int(nro);
	xuo = ealloc1(nro,sizeof(float*));
	zuo = ealloc1(nro,sizeof(float*));
	for (ir=0,jr=0; ir<nri; ++ir) {
		for (iu=0; iu<nui[ir]-1; ++iu,++jr) {
			aro[jr] = ari[ir];
			nuo[jr] = 2;
			xuo[jr] = ealloc1float(2);
			zuo[jr] = ealloc1float(2);
			xuo[jr][0] = xui[ir][iu];
			zuo[jr][0] = zui[ir][iu];
			xuo[jr][1] = xui[ir][iu+1];
			zuo[jr][1] = zui[ir][iu+1];
		}
		free1float(xui[ir]);
		free1float(zui[ir]);
	}
	free1float(ari);
	free1int(nui);
	free1(xui);
	free1(zui);

	/* output reflectors */
	*nr = nro;
	*ar = aro;
	*nu = nuo;
	*xu = xuo;
	*zu = zuo;
}

void makeref (float dsmax, int nr, float *ar, 
	int *nu, float **xu, float **zu, Reflector **r)
/*****************************************************************************
Make piecewise cubic reflectors
******************************************************************************
Input:
dsmax		maximum length of reflector segment
nr		number of reflectors
ar		array[nr] of reflector amplitudes
nu		array[nr] of numbers of (x,z) pairs; u = 0, 1, ..., nu[ir]
xu		array[nr][nu[ir]] of reflector x coordinates x(u)
zu		array[nr][nu[ir]] of reflector z coordinates z(u)

Output:
r		array[nr] of reflectors
******************************************************************************
Notes:
Space for the ar, nu, xu, and zu arrays is freed by this function, since
they are only used to construct the reflectors.

This function is meant to be called only once, so it need not be very
efficient.  Once made, the reflectors are likely to be used many times, 
so the cost of making them is insignificant.
*****************************************************************************/
{
	int ir,iu,nuu,iuu,ns,is;
	float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
		*u,*s,(*xud)[4],(*zud)[4],*us;
	ReflectorSegment *rs;
	Reflector *rr;
	
	/* allocate space for reflectors */
	*r = rr = ealloc1(nr,sizeof(Reflector));

	/* loop over reflectors */
	for (ir=0; ir<nr; ++ir) {

		/* compute cubic spline coefficients for uniformly sampled u */
		u = ealloc1float(nu[ir]);
		for (iu=0; iu<nu[ir]; ++iu)
			u[iu] = iu;
		xud = (float(*)[4])ealloc1float(4*nu[ir]);
		csplin(nu[ir],u,xu[ir],xud);
		zud = (float(*)[4])ealloc1float(4*nu[ir]);
		csplin(nu[ir],u,zu[ir],zud);

		/* finely sample x(u) and z(u) and compute length s(u) */