/* 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 clonesdecodeReflectors	parse reflectors parameter stringdecodeReflector		decode a particular reflectorbreakReflectors		break up reflectors duplicating interior (x,z) pointsmakeref			make piecewise cubic reflectors		raylv2		 	Trace ray between two points,				for linear velocity v = v00+dvdx*x+dvdz*zaddsinc		Add sinc wavelet to trace at specified time and			with specified amplitudemakericker	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 modelingtypedef 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 reflectorsaPtr		pointer to a specifying reflector amplitudesnxzPtr		pointer to nxz specifying number of (x,z) pairs defining the		reflectorsxPtr		pointer to array[x][nr] of x values for the entire modelzPtr		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 reflectorOutput:aPtr		pointer to amplitude valuenxzPtr		pointer to number of x,z pairsxPtr		array of x values for one reflectorzPtr		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 segmentnr		number of reflectorsar		array[nr] of reflector amplitudesnu		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, sincethey are only used to construct the reflectors.This function is meant to be called only once, so it need not be veryefficient.  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) */