void unwrap_phase(int n, float w, float *phase)
/************************************************************************
unwrap_phase - unwrap the phase
*************************************************************************
Input:
n		number of samples
w		unwrapping flag; returns an error if w=0
phase		array[n] of input phase values

Output:
phase		array[n] of output phase values
*************************************************************************
Notes:
The phase is assumed to be continuously increasing. The strategy is
to look at the change in phase (dphase) with each time step. If it is larger
than PI/w, then use the previous value of dphase. No attempt is
made at smoothing the dphase curve.
*************************************************************************
Author: John Stockwell, CWP, 1994
************************************************************************/
{
	int i;
	float pibyw=0.0;
	float *dphase;
	float *temp;

	/* prevent division by zero in PI/w */
	if (w==0)  err("wrapping parameter is zero");
	else       pibyw = PI/w;

	/* allocate space */
	dphase = ealloc1float(n);
	temp = ealloc1float(n);

	/* initialize */
	temp[0]=phase[0];
	dphase[0]=0.0;

	/* compute unwrapped phase at each time step */
	for (i = 1; i < n; ++i) {

		/* compute jump in phase */
		dphase[i] = ABS(phase[i] - phase[i-1]);

		/* if dphase >= PI/w, use previous dphase value */
		if (ABS(dphase[i] - dphase[i-1]) >= pibyw )
			dphase[i] = dphase[i-1];

		/* sum up values in temporary vector */
		temp[i] = temp[i-1] + dphase[i];
	}

	/* assign values of temporary vector to phase[i] */
	for (i=0; i<n; ++i) phase[i] = temp[i];

	/* free space */
	free1float(temp);
	free1float(dphase);
}

void integ(float **mig,int nz,float dz,int nx,int m,float **migi) 
/* integration of a two-dimensional array	 
  input: 
    mig[nx][nz]		two-dimensional array
  output:
    migi[nx][nz+2*m] 	integrated array 
*/
{
	int nfft, nw, ix, iz, iw;
	float  *amp, dw, *rt;
	complex *ct;


        /* Set up FFT parameters */
        nfft = npfaro(nz+m, 2 * (nz+m));
        if (nfft >= SU_NFLTS || nfft >= 720720)
                err("Padded nt=%d -- too big", nfft);

        nw = nfft/2 + 1;
	dw = 2.0*PI/(nfft*dz);

	amp = ealloc1float(nw);
	for(iw=1; iw<nw; ++iw) 
		amp[iw] = 0.5/(nfft*(1-cos(iw*dw*dz)));
	amp[0] = amp[1];

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

	for(ix=0; ix<nx; ++ix) {
        	memcpy(rt, mig[ix], nz*FSIZE);
       	 	memset((void *) (rt + nz), 0, (nfft-nz)*FSIZE); 
        	pfarc(1, nfft, rt, ct);

        	/* Integrate traces   */
		for(iw=0; iw<nw; ++iw){
			ct[iw].i = ct[iw].i*amp[iw];
			ct[iw].r = ct[iw].r*amp[iw];
		}

        	pfacr(-1, nfft, ct, rt);

        	for (iz=0; iz<m; ++iz)  migi[ix][iz] = rt[nfft-m+iz];
        	for (iz=0; iz<nz+m; ++iz)  migi[ix][iz+m] = rt[iz];
	}

	free1float(amp);
	free1float(rt);
	free1complex(ct);
}

/* 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 makericker (float fpeak, float dt, Wavelet **w)
/*****************************************************************************
Make Ricker wavelet
******************************************************************************
Input:
fpeak		peak frequency of wavelet
dt		time sampling interval

Output:
w		Ricker wavelet
*****************************************************************************/
{
	int iw,lw,it,jt;
	float t,x,*wv;
	
	iw = -(1+1.0/(fpeak*dt));
	lw = 1-2*iw;
	wv = ealloc1float(lw);
	for (it=iw,jt=0,t=it*dt; jt<lw; ++it,++jt,t+=dt) {
		x = PI*fpeak*t;
		x = x*x;
		wv[jt] = exp(-x)*(1.0-2.0*x);
	}
	*w = ealloc1(1,sizeof(Wavelet));
	(*w)->lw = lw;
	(*w)->iw = iw;
	(*w)->wv = wv;
}

static void mkvrms (int ndmo, float *tdmo, float *vdmo,
	int nt, float dt, float ft, float *vrms)
/*****************************************************************************
make uniformly sampled vrms(t) for DMO
******************************************************************************
Input:
ndmo		number of tdmo,vdmo pairs
tdmo		array[ndmo] of times
vdmo		array[ndmo] of rms velocities
nt		number of time samples
dt		time sampling interval
ft		first time sample

Output:
vrms		array[nt] of rms velocities
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 10/03/91
*****************************************************************************/
{
	int it;
	float t,(*vdmod)[4];

		vdmod = (float(*)[4])ealloc1float(ndmo*4);
	cmonot(ndmo,tdmo,vdmo,vdmod);
	for (it=0,t=ft; it<nt; ++it,t+=dt)
		intcub(0,ndmo,tdmo,vdmod,1,&t,&vrms[it]);
	free1float((float*)vdmod);
}

void twindow(int nt, int wtime, float *data)
/************************************************************
twindow - simple time gating
*************************************************************
Input:
nt	number of time samples
wtime	= n*dt   where n are integer ex=1,2,3,4,5,...
          wtime=3 as default
 used for Frequency Weighted and Thin-bed attributes
*************************************************************
Author:	UGM (Geophysics Students): Agung Wiyono, 2005
************************************************************/
{
	float val;
	float *temp;
	int i;
	float sum;
	int nwin;
	
	nwin=2*wtime+1;
	temp = ealloc1float(nt);
	sum=0.0;
	for (i = 0; i< wtime+1; ++i) {
		val = data[i];
		sum +=val;
	}
	/* weighted */
	temp[0] = sum/nwin;
	
	/* dt<wtime */
	for (i = 1; i < wtime; ++i) {
		val = data[i+wtime];
		sum+=val;
		++nwin;
		temp[i] = sum/nwin;
		}
	/*wtime<dt<dt-wtime */
	for (i = wtime ; i < nt-wtime; ++i) {
		val = data[i+wtime];
		sum += val;
		val = data[i-wtime];
		sum -=val;
		temp[i] = sum/nwin;
	}

	/*dt-wtime<dt*/
	for (i = nt - wtime; i < nt; ++i) {
		val = data[i-wtime];
		sum -= val;
		--nwin;
		temp[i] = sum/nwin;
	}
	
	
	for (i=0; i<nt; ++i) data[i] = temp[i];

	/* Memori free */
	free1float(temp);
}

int
main(int argc, char **argv)
{
	int j,nt,flag,ntout;
	float *buf,*ttn,dt,dtout=0.0,tmin,tmax;

	/* Initialize */
	initargs(argc, argv);
	requestdoc(1);

	/* Get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = (float) tr.dt/1000000.0;

	if (!getparfloat("tmin", &tmin)) tmin=0.1*nt*dt;
	if (!getparint("flag", &flag)) flag=1;
	if(flag==1) {
		dtout=tmin*2.*dt;
		tmax=nt*dt;
		ntout=1+tmax*tmax/dtout; CHECK_NT("ntout",ntout);
		ttn=ealloc1float(ntout);
		for(j=0;j<ntout;j++) ttn[j]=sqrt(j*dtout);
	}else{
		if (!getparfloat("dt", &dt)) dtout=0.004;
		ntout=1+sqrt(nt*dt)/dtout; CHECK_NT("ntout",ntout);
		ttn=ealloc1float(ntout);
		for(j=0;j<ntout;j++) ttn[j]=j*j*dtout*dtout;
	}
	buf = ealloc1float(nt);

	fprintf(stderr,"sutsq: ntin=%d dtin=%f ntout=%d dtout=%f\n",
		nt,dt,ntout,dtout);

	/* Main loop over traces */
	do {
		for(j=0;j<nt;j++) buf[j]=tr.data[j];
		tr.ns = ntout;
		tr.dt = dtout*1000000.;			
		ints8r(nt,dt,0.,buf,0.0,0.0,
			ntout,ttn,tr.data);
		puttr(&tr);
	} while (gettr(&tr));
	
	return(CWP_Exit());
}

void tabtrcoefs(int ninf, float *rho, float *v, float **theta,
					float *dip, float *trcoefs)
/*****************************************************************************
table transmission coefficients
******************************************************************************
Input:
ninf		x coordinate of source (must be within x samples)
xxxx

Output:
trcoefs		array[1..ninf] containing transmission coefficients
******************************************************************************
Notes:
The parameters are exactly the same as in the main program.
This is a subroutine so that the temporary arrays may be disposed
of after exit. This routine is only called once.
******************************************************************************
Author:  Brian Sumner, Colorado School of Mines, 1985
******************************************************************************/
{
/**
* Local variables:
*    TEMP, TMP - temporary arrays
*    I, J - loop variables
*    T1, T2 - temporaries
*    R, P - more temporaries
**/

	int i,j;
	float t1, t2, r, p, *temp, **tmp;

	temp = ealloc1float(ninf+1);
	tmp = ealloc2float(ninf+1, ninf+1);
	
	for (i = 0; i <= ninf; ++i)  temp[i] = rho[i]*v[i];

	for (j = 1; j <= ninf; ++j) {
		for (i = 1; i < j; ++i) {
			t1 = temp[i]*cos(theta[j][i-1]+dip[i]);
			t2 = temp[i-1]*cos(theta[j][i]+dip[i]);
			r = (t1 - t2)/(t1 + t2);
			tmp[j][i] = 1.0 - r*r;
		}
	}

	for (j = 1; j <= ninf; ++j) {
		t1 = temp[j];
		t2 = temp[j-1];
		p = (t1 - t2)/(t1 + t2);
		for (i = 1; i < j; ++i)  p *= tmp[j][i];
		trcoefs[j] = p;
	}
	
	/* free space */
	free1float(temp);
	free2float(tmp);
}

void remove_fb(float *yp,float *ym,int n,short *scaler,short *shft)
/* Find Scale and timeshift 

	Yp = data trace
	Ym = wavelet
	
	F=min(Yp(x) - Ym(x)*a)^2 is a function to minimize for scaler
	
	find shift first with xcorrelation then
	
	solve for a - scaler 
*/



{

#define RAMPR 5

	void find_p(float *ym,float *yp,float *a,int *b,int n);
	
	int it,ir;
	
	float a=1.0;
	int b=0;
	int ramps;
	float *a_ramp;
	
	ramps=NINT(n-n/RAMPR);
	
	find_p(ym,yp,&a,&b,n);
	
	a_ramp = ealloc1float(n);
	
	for(it=0;it<ramps;it++)
		a_ramp[it]=1.0;
	
	for(it=ramps,ir=0;it<n;it++,ir++) {
		a_ramp[it]=1.0-(float)ir/(float)(n-ramps-1);
	}
	
	
	
	if (b<0) {
		for(it=0;it<n+b;it++) 
			yp[it-b] -=ym[it]*a*a_ramp[it-b];
	} else {
		for(it=0;it<n-b;it++) 
			yp[it] -=ym[it-b]*a*a_ramp[it+b];
	}
	
	*scaler = NINT((1.0-a)*100.0);
	*shft = b;
	
	free1float(a_ramp);
}

void find_p(float *ym,float *yp,float *a,int *b,int n)
{
#define NP 1

	float *y,**x;
	int n_s;
	
	int k;
	float *res;
	int jpvt[NP];
	float qraux[NP];
	float work[NP]; 
	
														
	x = ealloc2float(n,1);
	y = ealloc1float(n);
	res = ealloc1float(n);
	
	memcpy((void *) &x[0][0], (const void *) &ym[0], n*FSIZE);
	memset((void *) y, (int) '\0', n*FSIZE);	
	
	/* Solve for shift */
 	xcor (n,0,ym,n,0,yp,n,-n/2,y);
	/* pick the maximum */
	*b = -max_index(n,y,1)+n/2;

	n_s = n-abs(*b);
 	
	if (*b < 0) {
		memcpy((void *) &x[0][0], (const void *) &ym[*b], n_s*FSIZE);
	} else {
		memcpy((void *) &x[0][*b], (const void *) &ym[0], n_s*FSIZE);
													
	}
	/* Solve for scaler */
	sqrst(x, n_s, 1,yp,0.0,a,res,&k,&jpvt[0],&qraux[0],&work[0]);
	
	
	
	free2float(x);
	free1float(res);
	free1float(y);
}

int main( int argc, char *argv[] )
{
	int n1;		/* number of x y values */
	int stinc;	/* x increment  */
	int f;		/* filter length */
	int m;		/* filter method flag */
	float *x;	/* array of x index values */
	float *y; 	/* array of y values */

        /* Initialize */
        initargs(argc, argv);
        requestdoc(1);

    
	MUSTGETPARINT("n1",&n1);
	if( !getparint("stinc",&stinc)) stinc=1;
	if( !getparint("f",&f)) f=5;
	if( !getparint("m",&m)) m=1;
	
	/* allocate arrays */
	x = ealloc1float(n1);
	y = ealloc1float(n1);
	
	/* Read data into the arrays */
	{ int i;
		for(i=0;i<n1;i++) {
			fscanf(stdin," %f %f\n",&x[i],&y[i]);
		}
	}
	/* smooth */
	sm_st(x,y,n1,f,stinc,m);
	
	/* Write out */
	{ int i;
		for(i=0;i<n1;i++) {
			fprintf(stdout," %10.3f %10.3f\n",x[i],y[i]);
		}
	}
	
	free1float(x);
	free1float(y);
   	return EXIT_SUCCESS;
}

int
main(int argc, char **argv)
{
	int nt;				/* number of samples on input */
	int ntout;			/* number of samples on output */
	int it;				/* counter */
	int istart;			/* beginning sample */
	int izero;			/* - istart */
	int norm;			/* user defined normalization value */
	int sym;			/* symmetric plot? */
	float scale;			/* scale factor computed from norm */
	float *temp=NULL;		/* temporary array */
	float dt;			/* time sampling interval (sec) */

	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(1);

	/* get information from the first header */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = tr.dt/1000000.0;

	/* get parameters */
	if (!getparint("ntout",&ntout)) ntout=101;
	if (!getparint("norm",&norm)) norm = 1;
	if (!getparint("sym",&sym)) sym = 1;
        checkpars();
	
	/* allocate workspace */
	temp = ealloc1float(ntout);
	
	/* index of first sample */
	if (sym == 0) istart = 0;
	else istart = -(ntout-1)/2;

	/* index of sample at time zero */
	izero = -istart;
	
	/* loop over traces */
	do {
		xcor(nt,0,tr.data,nt,0,tr.data,ntout,istart,temp);
		if (norm) {
			scale = 1.0/(temp[izero]==0.0?1.0:temp[izero]);
			for (it=0; it<ntout; ++it)  temp[it] *= scale;
		}
		memcpy((void *) tr.data, (const void *) temp, ntout*FSIZE);
		tr.ns = ntout;
		tr.f1 = -dt*ntout/2.0;
		tr.delrt = 0;
		puttr(&tr);
	} while(gettr(&tr));

	return(CWP_Exit());
}

int main( int argc, char *argv[] )
{
	/* Segy data constans */
	int nt;                 /* number of time samples               */
        int ntr=0;              /* number of traces                     */
	
	float *filter;
	int fnl,fnr;
	int fnp;
	int fld;
	int fm;
        float dt;               /* sample interval in secs              */
	float prw;		/* pre-withening */
	
	initargs(argc, argv);
   	requestdoc(1);
	
        /* get information from the first header */
        if (!gettr(&tr)) err("can't get first trace");
        nt = tr.ns;

        if (!getparfloat("dt", &dt)) dt = ((double) tr.dt)/1000000.0;
        if (!dt) {
                dt = .002;
                warn("dt not set, assumed to be .002");
        }
	
	if(!getparint ("fnl", &fnl)) fnl=15;
	fnr=fnl;
	if(!getparint ("fnp", &fnp)) fnp=fnr+fnl+fnr/2;
	if(!getparfloat ("prw", &prw)) prw=1.0;
		
	if(fnl!=0) {
		fld=0; fm=0; fnr=fnl;
		filter = ealloc1float(fnp);
		SG_smoothing_filter(fnp,fnl,fnr,fld,fm,filter); 
/*		rwa_smoothing_filter(1,fnl,fnr,filter); */ 
	} else {
		filter= NULL;
	}
	
	do {
		do_minphdec(tr.data,nt,filter,fnl,fnr,prw);
		
		tr.ns=nt;
		ntr++;		
		puttr(&tr);
	} while(gettr(&tr));
	
   return EXIT_SUCCESS;
}

int
main(int argc, char **argv)
{
	char *outpar;		/* name of file holding output parfile	*/
	FILE *outparfp;		/* ... its file pointer			*/
	int n1;			/* number of floats per line		*/
	size_t n1read;		/* number of items read			*/
	size_t n2 = 0;		/* number of lines in input file 	*/
	float *z;		/* binary floats			*/

	/* Hook up getpar */
	initargs(argc, argv);
	requestdoc(1);

	/* Get parameters and do set up */
	if (!getparstring("outpar", &outpar))	outpar = "/dev/tty" ;
	outparfp = efopen(outpar, "w");
	MUSTGETPARINT("n1",&n1);

	z = ealloc1float(n1);

	/* Loop over data converting to ascii */
	while ((n1read = efread(z, FSIZE, n1, stdin))) {
		register int i1;

		if (n1read != n1)
			err("out of data in forming line #%d", n2+1);
		for (i1 = 0; i1 < n1; ++i1)
/* z2xyz.c:70: warning: format ‘%d’ expects type ‘int’, but argument 2 has type ‘size_t’ */
/* 			printf("%d %d %11.4e \n",n2,i1,z[i1]); */
#if __WORDSIZE == 64
			printf("%lu %d %11.4e \n",n2,i1,z[i1]);
#else
			printf("%u %d %11.4e \n",n2,i1,z[i1]);
#endif
		++n2;
	}


	/* Make par file */
/* z2xyz.c:76: warning: format ‘%d’ expects type ‘int’, but argument 3 has type ‘size_t’ */
/* 	fprintf(outparfp, "n2=%d\n", n2); */
#if __WORDSIZE == 64
	fprintf(outparfp, "n2=%lu\n", n2);
#else
	fprintf(outparfp, "n2=%u\n", n2);
#endif

	return(CWP_Exit());
}

int
main (int argc, char **argv)
{
	int n1,n2,i2;
	float f1,f2,d1,d2,*x;
	char *label2="Trace",label[256];
	FILE *infp=stdin,*outfp=stdout;

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);

	/* get optional parameters */
	if (!getparint("n1",&n1)) {
		if (efseeko(infp,(off_t) 0,SEEK_END)==-1)
			err("input file size is unknown; specify n1!\n");
		if ((n1=((int) (eftello(infp)/((off_t) sizeof(float)))))<=0)
			err("input file size is unknown; specify n1!\n");
		efseeko(infp,(off_t) 0,SEEK_SET);
	}

	if (!getparfloat("d1",&d1)) d1 = 1.0;
	if (!getparfloat("f1",&f1)) f1 = d1;
	if (!getparint("n2",&n2)) n2 = -1;
	if (!getparfloat("d2",&d2)) d2 = 1.0;
	if (!getparfloat("f2",&f2)) f2 = d2;
	getparstring("label2",&label2);

	/* allocate space */
	x = ealloc1float(n1);

	/* loop over 2nd dimension */
	for (i2=0; i2<n2 || n2<0; i2++) {

		/* read input array, watching for end of file */
		if (efread(x,sizeof(float),n1,infp)!=n1) break;
			
		/* make plot label */
		sprintf(label,"%s %0.4g",label2,f2+i2*d2);

		/* plot the array */
		prp1d(outfp,label,n1,d1,f1,x);
	}
	
	return(CWP_Exit());
}

void differentiate(int n, float h, float *f)
/************************************************************************
differentiate - compute the 1st derivative of a function f[]
************************************************************************
Input:
n		number of samples
h		sample rate
f		array[n] of input values

Output:
f		array[n], the derivative of f
************************************************************************
Notes:
This is a simple 2 point centered-difference differentiator.
The derivatives at the endpoints are computed via 2 point leading and
lagging differences. 
************************************************************************
Author: John Stockwell, CWP, 1994
************************************************************************/
{
	int i;	
	float *temp;
	float h2=2*h;

	/* allocate space in temporary vector */
	temp = ealloc1float(n);

	/* do first as a leading difference */
	temp[0] = (f[1] - f[0])/h;

	/* do the middle values as a centered difference */
	for (i=1; i<n-1; ++i) temp[i] = (f[i+1] - f[i-1])/h2;

	/* do last value as a lagging difference */
	temp[n-1] = (f[n-1] - f[n-2])/h;

	for (i=0 ; i < n ; ++i) f[i] = temp[i];

	free1float(temp);
}

static void makezt (int nz, float dz, float fz, float v[], 
	int nt, float dt, float ft, float z[])
/************************************************************************
makezt - compute z(t) from v(z)
*************************************************************************
Input:
nz	number of z values (output)
dz	depth sampling interval (output)
fz	first depth value (output)
v[]	array of velocities as a function of time t
nt	number of time samples
dt	time sampling interval
ft	first time sample
Output:
z[]	array of z values as a function of t
*************************************************************************
Author: CWP: based on maketz by Dave Hale (c. 1992)
*************************************************************************/
{
	int iz;			/* counter */
	float vfz;		/* velocity at the first depth sample */
	float vlz;		/* velocity at the last depth sample */
	float *t=NULL;		/* array of time values as a function of z */

	/* allocate space */
	t = ealloc1float(nz);

	/* calculate t(z) from v(z) */
	t[0] = 2.0*fz/v[0];
	for (iz=1; iz<nz; ++iz)
		t[iz] = t[iz-1]+2.0*dz/v[iz-1];
	vfz = v[0];
	vlz = v[nz-1];

	/* compute z(t) from t(z) */
	tzzt(nz,dz,fz,t,vfz,vlz,nt,dt,ft,z);
	free1float(t);
}

void bandpass(float *data, int nt, int nfft, int nfreq,
		float *filterj, float *ftracej)
{
	float *rt;
	complex *ct;
	int i;

	rt  = ealloc1float(nfft);
	ct  = ealloc1complex(nfreq);
        /* Load trace into rt (zero-padded) */
        memcpy((char*) rt, (char*) data, nt*FSIZE);
        bzero(rt + nt, (nfft-nt)*FSIZE);

        /* FFT, filter, inverse FFT */
        pfarc(1, nfft, rt, ct);
        for (i = 0; i < nfreq; ++i)  ct[i] = crmul(ct[i], filterj[i]);
        pfacr(-1, nfft, ct, rt);

        /* Load traces back in, recall filter had nfft factor */
        for (i = 0; i < nt; ++i)  ftracej[i] = rt[i]; /* ftracej = rt ?? */
	free(rt);
	free(ct);
}

void do_smooth(float *data, int nt, int isl)
/**********************************************************************
do_smooth - smooth data in a window of length isl samples
**********************************************************************
Input:
data[]		array of floats of size nt
nt		size of array
isl		integerized window length
Output:
returns smoothed data.

**********************************************************************
Author: Nils Maercklin,
 	 GeoForschungsZentrum (GFZ) Potsdam, Germany, 2001.
 	 E-mail: [email protected]
**********************************************************************/
{
	register int it,jt;
	float *tmpdata, sval;

	tmpdata=ealloc1float(nt);
	for (it=0;it<nt;it++) {
	sval=0.0;
	if ( (it >= isl/2) && (it < nt-isl/2) ) {
		for (jt=it-isl/2;jt<it+isl/2;jt++) {
			sval += data[jt];
		}
		tmpdata[it] = sval / (float) isl;
		} else {
			tmpdata[it] = data[it];
		}
	}

	memcpy((void *) data, (const void *) tmpdata, nt*FSIZE);

	free1float(tmpdata);
}

void do_smooth(float *data, int nt, int isl)
{
    int it,jt;
    float *tmpdata, sval;
    
    tmpdata=ealloc1float(nt);

    for (it=0;it<nt;it++) {
          sval=0.0;
          if ( (it >= isl/2) && (it < nt-isl/2) ) {
            for (jt=it-isl/2;jt<it+isl/2;jt++) {
                  sval += data[jt];
            }
            tmpdata[it] = sval / (float) isl;
      }
      else {
            tmpdata[it] = 0.0;
      }
    }
    for (it=0;it<nt;it++) {
          data[it] = tmpdata[it];
    }
    free1float(tmpdata);
}