par_upweik.c

该程序是用vc开发的对动态数组进行管理的DLL

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

#include "par.h"

/*********************** self documentation **********************/
/*****************************************************************************
UPWEIK - Upwind Finite Difference Eikonal Solver

eikpex - Eikonal equation extrapolation of times and derivatives in 
         polar coordinates
ray_theoretic_sigma - difference equation extrapolation of "ray_theoretic_sigma" in polar coordinates
ray_theoretic_beta - difference equation extrapolation of "ray_theoretic_beta" in polar coordinates
eiktam - Compute traveltimes t(x,z) and  propagation angle a(x,z) via 
         eikonal equation, and ray_theoretic_sigma sig(x,z), incident angle bet(x,z) 
         via Crank-Nicolson Method
******************************************************************************
Function Prototypes:
void eikpex (int na, float da, float r, float dr, 
	float sc[], float uc[], float wc[], float tc[],
	float sn[], float un[], float wn[], float tn[]);
void ray_theoretic_sigma (int na, float da, float r, float dr, 
	float uc[], float wc[], float sc[],
	float un[], float wn[], float sn[]);
void ray_theoretic_beta (int na, float da, float r, float dr, 
	float uc[], float wc[], float bc[],
	float un[], float wn[], float bn[]);
void eiktam (float xs, float zs, 
	int nz, float dz, float fz, int nx, float dx, float fx, float **vel,
	float **time, float **angle, float **sig, float **bet)
******************************************************************************
eikpex:
Input:
na		number of a samples
da		a sampling interval
r		current radial distance r
dr		radial distance to extrapolate
sc		array[na] of slownesses at current r
uc		array[na] of dt/dr at current r
wc		array[na] of dt/da at current r
tc		array[na] of times t at current r
sn		array[na] of slownesses at next r

Output:
un		array[na] of dt/dr at next r (may be equivalenced to uc)
wn		array[na] of dt/da at next r (may be equivalenced to wc)
tn		array[na] of times t at next r (may be equivalenced to tc)

******************************************************************************
ray_theoretic_sigma:
Input:
na		number of a samples
da		a sampling interval
r		current radial distance r
dr		radial distance to extrapolate
uc		array[na] of dt/dr at current r
wc		array[na] of dt/da at current r
sc		array[na] of ray_theoretic_sigma  at current r
un		array[na] of dt/dr at next r
wn		array[na] of dt/da at next r

Output:
sn		array[na] of ray_theoretic_sigma at next r 
******************************************************************************
ray_theoretic_beta:
Input:
na		number of a samples
da		a sampling interval
r		current radial distance r
dr		radial distance to extrapolate
uc		array[na] of dt/dr at current r
wc		array[na] of dt/da at current r
bc		array[na] of ray_theoretic_beta  at current r
un		array[na] of dt/dr at next r
wn		array[na] of dt/da at next r

Output:
bn		array[na] of ray_theoretic_beta at next r 
******************************************************************************
eiktam:
Input:
xs		x coordinate of source (must be within x samples)
zs		z coordinate of source (must be within z samples)
nz		number of z samples
dz		z sampling interval
fz		first z sample
nx		number of x samples
dx		x sampling interval
fx		first x sample
vel		array[nx][nz] containing velocities

Output:
time		array[nx][nz] containing first-arrival times
angle		array[nx][nz] containing propagation angles
sig  		array[nx][nz] containing ray_theoretic_sigmas
bet		array[nx][nz] containing ray_theoretic_betas
******************************************************************************
Notes:
eikpex:
If na*da==2*PI, then the angular coordinate is wrapped around (periodic). 

This function implements the finite-difference method described by Bill
Symes (Rice University) and Jos van Trier (Stanford University) in a
(1990) preprint of a paper submitted to Geophysics.

ray_theoretic_sigma:
This routine implements the Crank-Nicolson finite-difference method with
boundary conditions dray_theoretic_sigma/da=0.

ray_theoretic_beta:
This function implements the Crank-Nicolson finite-difference 
method, with boundary conditions dray_theoretic_beta/da=1. 

eiktam:
The actual computation of times and ray_theoretic_sigmas is done in polar coordinates,
with bilinear interpolation used to map to/from rectangular coordinates.

******************************************************************************
Authors: CWP: Zhenuye Liu, based on code by Dave Hale, 1992.
******************************************************************************/

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

#define TINY 1.0e-3f	/* avoid divide by zero */
#define CFL 0.98f	/* Courant/Friedrichs/Lewy stability factor */

void eikpex (int na, float da, float r, float dr, 
	float sc[], float uc[], float wc[], float tc[],
	float sn[], float un[], float wn[], float tn[])
/*****************************************************************************
eikpex - Eikonal equation extrapolation of times and derivatives in 
         polar coordinates
******************************************************************************
Input:
na		number of a samples
da		a sampling interval
r		current radial distance r
dr		radial distance to extrapolate
sc		array[na] of slownesses at current r
uc		array[na] of dt/dr at current r
wc		array[na] of dt/da at current r
tc		array[na] of times t at current r
sn		array[na] of slownesses at next r

Output:
un		array[na] of dt/dr at next r (may be equivalenced to uc)
wn		array[na] of dt/da at next r (may be equivalenced to wc)
tn		array[na] of times t at next r (may be equivalenced to tc)
******************************************************************************
Notes:
If na*da==2*PI, then the angular coordinate is wrapped around (periodic). 

This function implements the finite-difference method described by Bill
Symes (Rice University) and Jos van Trier (Stanford University) in a
(1990) preprint of a paper submitted to Geophysics.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 07/16/90
******************************************************************************/
{
	int i,wrap;
	float drleft,drorig,frac,cmax,umaxl,uminr,uminm,umaxm,
		uu,unew,uold,ueol,ueor,wor,or,*wtemp,*s;
	
	/* allocate workspace */
	wtemp = alloc1float(na);
	s = alloc1float(na);
	
	/* remember the step size */
	drleft = drorig = dr;
	
	/* initialize slownesses to values at current r */
	for (i=0; i<na; ++i)
		s[i] = sc[i];
	
	/* copy inputs to output */
	for (i=0; i<na; ++i) {
		un[i] = uc[i];
		wn[i] = wc[i];
		tn[i] = tc[i];
	}
	
	/* determine if angular coordinate wraps around */
	wrap = ABS(na*da-2.0*PI)<0.01*ABS(da);
	
	/* loop over intermediate steps with adaptive stepsize */
	while (drleft>0.0) {
		
		/* determine adaptive step size according to CFL condition */
		for (i=0,cmax=TINY; i<na; ++i) {
			if (r*ABS(un[i])<TINY*ABS(wn[i]))
				cmax = 1.0f/TINY;
			else
				cmax = MAX(cmax,ABS(wn[i]/(r*un[i])));
		}
		dr = MIN(drleft,CFL/cmax*r*da);
		
		/* if angles wrap around */
		if (wrap) {
			umaxl = (wn[na-1]>0.0 ? un[na-1] : s[0]);
			if (wn[0]>0.0) {
				uminm = s[0];
				umaxm = un[0];
			} else {
				uminm = un[0];
				umaxm = s[0];
			}
			uminr = (wn[1]>0.0 ? s[0] : un[1]);
			ueol = uminm+umaxl;
			ueor = uminr+umaxm;
			wtemp[0] = wn[0]+dr*(ueor-ueol)/da;
			umaxl = (wn[na-2]>0.0 ? un[na-2] : s[na-1]);
			if (wn[na-1]>0.0) {
				uminm = s[na-1];
				umaxm = un[na-1];
			} else {
				uminm = un[na-1];
				umaxm = s[na-1];
			}
			uminr = (wn[0]>0.0 ? s[na-1] : un[0]);
			ueol = uminm+umaxl;
			ueor = uminr+umaxm;
			wtemp[na-1] = wn[na-1]+dr*(ueor-ueol)/da;
		
		/* else, if angles do not wrap around */
		} else {
			if (wn[0]<=0.0)
				wtemp[0] = wn[0] + 
					dr*(un[1]-un[0])/da; 
			else
				wtemp[0] = 0.0;
			if (wn[na-1]>=0.0) 
				wtemp[na-1] = wn[na-1] +
					dr*(un[na-1]-un[na-2])/da;
			else
				wtemp[na-1] = 0.0;
		}
		
		/* update interior w values via Enquist/Osher scheme */
		for (i=1; i<na-1; ++i) {
			umaxl = (wn[i-1]>0.0 ? un[i-1] : s[i]);
			if (wn[i]>0.0) {
				uminm = s[i];
				umaxm = un[i];
			} else {
				uminm = un[i];
				umaxm = s[i];
			}
			uminr = (wn[i+1]>0.0 ? s[i] : un[i+1]);
			ueol = uminm+umaxl;
			ueor = uminr+umaxm;
			wtemp[i] = wn[i]+dr*(ueor-ueol)/da;
		}
		
		/* decrement the size of step left to do */
		drleft -= dr;
		
		/* update radial coordinate and its inverse */
		r += dr;
		or = 1.0f/r;
		
		/* linearly interpolate slowness for new r */
		frac = drleft/drorig;
		for (i=0; i<na; ++i)
			s[i] = frac*sc[i]+(1.0f-frac)*sn[i];
		
		/* update w and u; integrate u to get t */
		for (i=0; i<na; i++) {
			wn[i] = wtemp[i];
			wor = wn[i]*or;
			uu = (s[i]-wor)*(s[i]+wor);
			if(uu<=0) printf("\tRaypath has a too large curvature!\n\t A smoother velocity is required. \n");
 			unew = (float)sqrt(uu); 
			uold = un[i];
			un[i] = unew;
			tn[i] += 0.5f*dr*(unew+uold);
		}
	}
	
	/* free workspace */
	free1float(wtemp);
	free1float(s);
}

void ray_theoretic_sigma (int na, float da, float r, float dr,