/* img_filter_2007.c
 *
 * Filter an img file, using cosine transform in Y direction and
 * cosine blending of windows.  Use fftpack routines cost and rfft.
Built upon img_filter_99, with the following modifications:

Add -Y switch to support extended latitude range files 
Add feature to use B-spline filters with knot # -K
Add flexibility to define A parameter in W2 filter
WHFS, 1 April 2006

In 2007 version (February 8) WHFS added Gaussian and spherical harmonic
taper filters to try to reproduce all the old features of img_filter.f
 *
 */
 
#include <stdio.h>
#include <math.h>
#include <stdlib.h>		/* apparently atof() wasn't  working without this */
#include <string.h>		/* for memmove()  */

struct FPARAMS {
	double	depth;
	double	crust;
	double	flex_lambda;
	double	gauss_sigma;
	double	w2a;
	int		do_flex;
	int		do_w2;
	int		do_gauss;
	int		do_harm;
	int		n1;
	int		n2;
	int		knot;
};

int	main (int argc, char **argv) {

	struct FPARAMS fparams;
	double	dxm = -1.0;
	double	circum = 40075.01229;	/* Equatorial circumference, km */
	double	f = 1.0/298.257;
	double	esq, d2r, radius, phi, width, height, sinphi, wf;
	float	*blender, *z, *wrkx, *wrky, *tpr, *kx, *ky, *transp;
	short int	*s;
	int	tracks = 0, verbose = 0, error = 0, tall = 0;
	int	nximg, nyimg, nyfft = -1, noverflow = 0;
	int	i, iswath, nswaths, nyimghalf, nyffthalf, nyfftp1;
	size_t	kbyte;
	int	kbytesum;
	char	*infile, *outfile;
	
	void	load_wavenumbers(), costi_(), rffti_(), load_data(), blend_and_write();
	void	rfft2dg1c2(), filter();
	FILE	*fpin, *fpout;
	
	infile = (char *)NULL;
	outfile = (char *)NULL;
	
	fparams.depth = 0.0;	/* Depth in km, positive to downward continue */
	fparams.crust = 7.0;	/* Crust thickness in km, for W1 filter  */
	fparams.flex_lambda = 135.0;	/* Wavelength of flexure in km for W1 filter */
	fparams.gauss_sigma = 30.0;		/* (wavelength in km)/5.34 for Gaussian filter */
	fparams.w2a = 9500.0;	/* Constant A in W2 filter equation, in km^4  */
	fparams.do_flex = 0;
	fparams.do_w2 = 0;
	fparams.do_gauss = 0;
	fparams.do_harm = 0;
	fparams.n1 = 50;
	fparams.n2 = 70;
	fparams.knot = 0;
	

	/* Parse arg list:  */
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'D':
					fparams.depth = atof(&argv[i][2]);
					break;
				case 'F':
					fparams.do_flex = (argv[i][2] == 'h') ? 1 : -1;
					break;
				case 'G':
					fparams.do_gauss = (argv[i][2] == 'h') ? 1 : -1;
					fparams.gauss_sigma = atof (&argv[i][3]);
					break;
				case 'H':
					fparams.do_harm = (argv[i][2] == 'h') ? 1 : -1;
					sscanf(&argv[i][3], "%d/%d", &fparams.n1, &fparams.n2);
					break;
				case 'K':
					fparams.knot = atoi(&argv[i][2]);
					break;
				case 'M':
					dxm = atof(&argv[i][2]);
					break;
				case 'N':
					nyfft = atoi(&argv[i][2]);
					break;
				case 'O':
					outfile = &argv[i][2];
					break;
				case 'T':
					tracks = 1;
					break;
				case 'V':
					verbose = 1;
					break;
				case 'W':
					fparams.do_w2 = 1;
					if ((sscanf(&argv[i][2], "%lf", &phi)) == 1) fparams.w2a = phi;
					break;
				case 'Y':
					tall = 1;
					break;
				default:
					fprintf (stderr, "img_filter_2007:  Unsupported option:  %s\n", argv[i]);
					error ++;
					break;
			}
		}
		else {
			infile = argv[i];
		}
	}
	if (dxm <= 0.0) {
		fprintf (stderr, "img_filter_2007:  You must supply -M<pixel width in minutes>\n");
		error ++;
	}
	if (nyfft <= 0) {
		fprintf (stderr, "img_filter_2007:  You must supply -N<nyfft>, the height of an FFT swath, in pixels.\n");
		error ++;
	}
	if (fparams.knot < 0) {
		fprintf (stderr, "img_filter_2007:  You must supply a positive knot to use option -K<knot>.\n");
		error ++;
	}
	if (nyfft%2 != 0) {
		fprintf (stderr, "img_filter_2007:  <nyfft> must be even.\n");
		error ++;
	}
	if (infile == (char *)NULL) {
		fprintf (stderr, "img_filter_2007:  You must supply an input file.\n");
		error ++;
	}
	if (outfile == (char *)NULL) {
		fprintf (stderr, "img_filter_2007:  You must supply -O<output file>.\n");
		error ++;
	}
	if (fparams.w2a <= 0.0) {
		fprintf (stderr, "img_filter_2007:  -W<a> must supply a positive a parameter.\n");
		error ++;
	}
	
	if (error) {
		fprintf (stderr, "usage:  img_filter_2007 <infile> -M<dxm> -N<nyfft> -O<outfile> [-D<d>] [-F<lorh>] [-G<lorh><sigma>] [-H<lorh><n1>/<n2>] [-K<k>] [-T] [-V] [-W[<a>]] [-Y]\n\n");
		fprintf (stderr, "\t<infile> and <outfile> are img files with pixel width <dxm> minutes.\n");
		fprintf (stderr, "\t<nyfft> is the height of each FFT swath, in pixels.  Must be even.\n");
		fprintf (stderr, "\tOption -D<d> downward continues to <d> km depth.\n");
		fprintf (stderr, "\tOption -Fh highpasses the flexure band (W1 by Eqn 9 of S&S 1994).\n");
		fprintf (stderr, "\tOption -Fl  lowpasses the flexure band (1-W1 by Eqn 9 of S&S 1994).\n");
		fprintf (stderr, "\tOption -Gh<sigma> does a Gaussian  hipass filter with sigma in km.\n");
		fprintf (stderr, "\tOption -Gl<sigma> does a Gaussian lowpass filter with sigma in km.\n");
		fprintf (stderr, "\t\tThe Gaussian sigma = (full wavelength at 1/2 amplitude) / 5.34\n");
		fprintf (stderr, "\tOption -Hl<n1>/<n2> lowpasses with spherical harmonic taper between degree n1 and n2\n");
		fprintf (stderr, "\tOption -Hh<n1>/<n2> hipasses  with spherical harmonic taper between degree n1 and n2\n");
		fprintf (stderr, "\t\tThese options expect n1 >= 0 and n2 >= n1.\n");
		fprintf (stderr, "\tOption -K<k> bandpasses with the k'th Bspline filter.\n");
		fprintf (stderr, "\tOption -T removes track coding from input file.\n");
		fprintf (stderr, "\tOption -V reports verbose operations to stderr.\n");
		fprintf (stderr, "\tOption -W applies W2 filter, including d if -D used (Eqn 11 of S&S 1994).\n");
		fprintf (stderr, "\tOption -W<a> as above but changes value of filter constant a [9500].\n");
		fprintf (stderr, "\tOption -Y means files have extended latitude range.\n");
		exit (-1);
	}


	/* Set up some stuff  */	
	if (dxm == 1.0) {
		nximg = 21600;
		nyimg = (tall) ? 17280 : 12672;
	}
	else if (dxm == 2.0) {
		nximg = 10800;
		nyimg = (tall) ? 8640 : 6336;
	}
	else {
		fprintf (stderr, "img_filter_2007:  ERROR.  Expecting -M<dxm> equal 1.0 or 2.0, not %.8lg\n", dxm);
		exit (-1);
	}
	nyimghalf = nyimg / 2;
	d2r = M_PI / 180.0;
	esq = f * (2.0 - f);
	radius = nximg/(2.0 * M_PI);
	if ( (2 * nyimg)%nyfft != 0) {
		fprintf (stderr, "img_filter_2007:  nyfft = %d  does not fit nyimg = %d\n",
			nyfft, nyimg);
		exit (-1);
	}
	nswaths = (2 * nyimg) / nyfft;
	nyfftp1 = nyfft + 1;
	nyffthalf = nyfft / 2;


	/* Malloc stuff:  */
	kbyte = nximg*nyfft/2;
	if ( (blender = (float *)calloc(kbyte, (size_t)4) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  calloc failure for blender array.\n");
		exit (-1);
	}
	kbytesum = kbyte * 4;
	
	kbyte = nximg*nyfftp1*4;
	if ( (z = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for z array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = (2*nximg + 16) * 4;
	if ( (wrkx = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for wrkx array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = (3*nyfftp1 + 16) * 4;
	if ( (wrky = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for wrky array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = nyffthalf*4;
	if ( (tpr = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for tpr array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = nximg*4;
	if ( (kx = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for kx array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = nyfftp1*4;
	if ( (ky = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for ky array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = nyfftp1*4;
	if ( (transp = (float *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for transp array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	kbyte = nximg*nyfftp1*2;
	if ( (s = (short int *)malloc(kbyte) ) == NULL) {
		fprintf (stderr, "img_filter_2007:  malloc failure for s array.\n");
		exit (-1);
	}
	kbytesum += kbyte;

	if (verbose) printf("img_filter_2007.c  %3.1lf files are %d by %d; nyfft is %d; %6.1lf Mbytes used.\n",
		dxm, nximg, nyimg, nyfft, kbytesum/1048576.0);
		
	if (verbose && fparams.knot) printf ("img_filter_2007:  Knot = %d filter centers on wavelength %lg km.\n", fparams.knot, pow(2.0, 9.0-fparams.knot));



	/* Initialize blending taper:  */
	wf = M_PI/nyffthalf;
	for (i = 0; i < nyffthalf; i++) tpr[i] = (float) (0.5 * (1.0 + cos(wf*(i-nyffthalf+0.5))));
	
	
	/* Initialize FT routines  */
	costi_(&nyfftp1, wrky);
	rffti_(&nximg, wrkx);

	/* Open files:  */
	if ( (fpin = fopen(infile, "r")) == NULL) {
		fprintf (stderr, "img_filter_2007:  Error.  Cannot open r %s\n", infile);
		exit (-1);
	}
	if ( (fpout = fopen(outfile, "w")) == NULL) {
		fprintf (stderr, "img_filter_2007:  Error.  Cannot open w %s\n", outfile);
		exit (-1);
	}


	/* Do swaths  */
	
	for (iswath = 0; iswath <= nswaths; iswath++) {
		phi = 2.0 * atan(exp( (nyimghalf - iswath * (nyfft/2) ) / radius) ) - M_PI_2;
		sinphi = sin(phi);
		width = circum * cos(phi) / sqrt(1.0 - esq * sinphi * sinphi);
		height = (2.0 * nyfft) * (width / nximg);	/* Cos FT doubles Y length  */
		if (verbose) printf ("Strip %2.2d of %d at latitude %8.4lf is %8.2lf by %7.2lf km.\n",
			iswath, nswaths, phi/d2r, width, height);
		
		load_data (iswath, nswaths, nximg, nyfftp1, tracks, s, z, fpin);
		load_wavenumbers (width, height, kx, ky, nximg, nyfftp1);
		rfft2dg1c2 (z, nximg, nyfftp1, wrkx, wrky, transp, -1);
		filter (z, kx, ky, nximg, nyfftp1, &fparams);
		rfft2dg1c2 (z, nximg, nyfftp1, wrkx, wrky, transp, 1);
		blend_and_write (iswath, nswaths, nximg, nyfft, &noverflow, blender, tpr, z, fpout);
	}
	/* close files and free mallocs  */
	if (noverflow) fprintf (stderr, "img_filter_2007:  WARNING.  %d pixels overflowed out of %d total\n", noverflow, nximg*nyimg);
	fclose (fpout);
	fclose (fpin);
	free ( (void *)s);
	free ( (void *)transp);
	free ( (void *)ky);
	free ( (void *)kx);
	free ( (void *)wrky);
	free ( (void *)wrkx);
	free ( (void *)z);
	free ( (void *)blender);
	exit (0);
}

void	load_wavenumbers (width, height, kx, ky, nximg, nyfftp1)
double	width, height;
float	*kx, *ky;
int	nximg, nyfftp1;
{
	int	i, j;
	double	f;

	/* Define kx,ky as 1.0/wavelength.
	   Set up kx for rfft, ky for cost.
	   Put kx-squared in kx[]; ky-squared in ky[].  */
	 
	kx[0] = 0.0;
	f = 1.0 / width;
	for (i = 1; i < nximg; i+= 2) {
		j = (i + 1)/2;
		kx[i] = (float) (j * f);
		kx[i+1] = kx[i];
	}
	if (nximg%2 == 0) kx[nximg-1] = (float) ( (nximg/2) * f);
	
	for (i = 1; i < nximg; i++) kx[i] = (kx[i] * kx[i]);

	ky[0] = 0.0;
	f = 1.0 / height;
	for (i = 1; i < nyfftp1; i++) ky[i] = (float) (i * f);

	for (i = 1; i < nyfftp1; i++) ky[i] = (ky[i] * ky[i]);
	
	return;
}

void load_data (iswath, nswaths, nximg, nyfftp1, tracks, s, z, fp)
int	iswath, nswaths, nximg, nyfftp1, tracks;
short int *s;
float	*z;
FILE	*fp;
{
	/* Load short int array s[] with nyfftp1 rows of data.
		Shift and shuffle data as needed, reading from fpin
		to load new data into s after shuffling.
		If (tracks), subtract 1 from odd values after read.
		Then copy contents of s to z.
	*/

	int	k, kstart, kstop, kread, kmove, nyffthalf, nxy;

	nyffthalf = (nyfftp1 - 1)/2;
	nxy = nyfftp1 * nximg;
	
	if (iswath == 1 || iswath == nswaths) {
		/* Re-use previous data; no need to shuffle or read.  */
		for (k = 0; k < nxy; k++) z[k] = (float)s[k];
		return;
	}
	
	if (iswath == 0) {
		/* Need to load entire amount:  */
		kread = nximg * nyfftp1;
		if ( (fread( (void *)s, (size_t)2, kread, fp) ) != kread) {
			fprintf (stderr, "img_filter_2007:  Read error at iswath = %d\n", iswath);
			exit (-1);
		}
		if (tracks) {
			/* Need to fix entire amount:  */
			for (k = 0; k < kread; k++) {
				if ( (abs((int)s[k]))%2 == 1) s[k] --;
			}
		}
		for (k = 0; k < nxy; k++) z[k] = (float)s[k];
		return;
	}

	if (iswath == nswaths - 1) {
		/* Discard, shuffle, and load one row less than usual.  */
		kstart = (nyffthalf - 1) * nximg;
		kmove = (nyffthalf + 2) * nximg * 2;
	}
	else {
		/* Need to discard first nyhalf rows, shuffle remainder, and load new on end.  */
		kstart = nyffthalf * nximg;
		kmove = (nyfftp1 - nyffthalf) * nximg * 2;
	}
	memmove ( (void *)s, (void *)&s[kstart], (size_t)kmove);
	kread = kstart;
	kstart = kmove / 2;
	if ( (fread( (void *)&s[kstart], (size_t)2, kread, fp) ) != kread) {
		fprintf (stderr, "img_filter_2007:  Read error at iswath = %d\n", iswath);
		exit (-1);
	}
	if (tracks) {
		kstop = nyfftp1 * nximg;
		for (k = kstart; k < kstop; k++) {
			if ( (abs((int)s[k]))%2 == 1) s[k] --;
		}
	}
	for (k = 0; k < nxy; k++) z[k] = (float)s[k];
	return;
}

void blend_and_write (iswath, nswaths, nximg, nyfft, noverflow, blender, tpr, z, fpout)
int	iswath, nswaths, nximg, nyfft, *noverflow;
float	blender[], tpr[], z[];
FILE	*fpout;
{
	int	nyffthalf, i, j, j1, j2, k, k1, k2;
	double	x;
	float	tpr1;
	short int h;

	nyffthalf = nyfft/2;
	
	if (iswath == 0) {
		/* Simply copy the first half of z to blender and return  */
		k2 = nyffthalf * nximg;
		for (k = 0; k < k2; k++) blender[k] = z[k];
		return;
	}
	
	
	j1 = (iswath >= nswaths - 1) ? 1 : 0;	/* Start of first half  */
	j2 = j1 + nyffthalf;			/* Start of 2nd half  */
	if (iswath == nswaths) {
		j1 = j2;	/* Need to blend 2nd half instead of 1st half  */
	}
	
	k1 = j1 * nximg;
	k2 = j2 * nximg;
	
	for (j = 0, k = 0; j < nyffthalf; j++) {
		tpr1 = 1.0 - tpr[j];
		for (i = 0; i < nximg; i++, k++) {
			x = rint(tpr[j] * z[k1+k] + tpr1 * blender[k]);
			if (fabs(x) > 32767.0) {
			/*	fprintf (stderr, "img_filter_2007:  Overflow problem.\n");  */
				*noverflow++;
				x = copysign(32767.0, x);
			}
			h = (short int)x;
			if ( (fwrite( (void*)&h, (size_t)2, (size_t)1, fpout) ) != 1) {
				fprintf (stderr, "img_filter_2007:  Write error.\n");
				exit (-1);
			}
			blender[k] = z[k2+k];
		}
	}
	return;
}


void	rfft2dg1c2 (a, n1, n2, wsave1, wsave2, transp, isgn)
float	*a, *wsave1, *wsave2, *transp;
int	n1, n2, isgn;
{
	/* A 2-d real array is stored contiguously in a[], with the n1
		dimension looped first (i.e., like fortran a(n1,n2) ).
		The n1 dimension is given a real fourier transform.
		The n2 dimension is given a real cosine  transform.
		If isgn = -1, the forward transforms are calculated;
		If isgn = +1, the inverse transforms are calculated,
		and the data are then normalized by the appropriate factor.

		The fftpack routines rfftf_, rfftb_, cost_ are used.
		It is assumed that wsave1 and wsave2 have been initialized
		by calling rffti_(&n1, wsave1) and costi_(&n2, wsave2), 
		respectively.
		
		transp[] is a work array of size n2 used to compute FTs in
		the n2 direction ("transpose").  */
		

	int	i, j, k, n;
	double	r;
	void	rfftf_(), rfftb_(), cost_();

	if (isgn < 0) {	

		/* Forward transform.  
			First do real ft of n1 dimension in place:  */
		n = n1;
		for (j = 0, k = 0; j < n2; j++, k+=n1) {
			rfftf_(&n, &a[k], wsave1);
		}

		/* Now do cosft of n2 dimension by column vectors in transp:  */
		n = n2;
		for (i = 0; i < n1; i++) {
			for (j = 0, k = i; j < n2; j++, k+=n1) {
				transp[j] = a[k];
			}
			cost_(&n, transp, wsave2);
			for (j = 0, k = i; j < n2; j++, k+=n1) {
				a[k] = transp[j];
			}
		}
	}
	else {
		/* Inverse transform.  Get normalization factor:  */
		n = 2*n1*(n2-1);
		r = 1.0/n;
		
		/* Do cosft of n2 dimension by column vectors in transp,
			and multiply by r to normalize:  */
		n = n2;
		for (i = 0; i < n1; i++) {
			for (j = 0, k = i; j < n2; j++, k+=n1) {
				transp[j] = (float) (r * a[k]);
			}
			cost_(&n, transp, wsave2);
			for (j = 0, k = i; j < n2; j++, k+=n1) {
				a[k] = transp[j];
			}
		}
		
		/* Now, do real ft of n1 dimension in place:  */
			
		n = n1;
		for (j = 0, k = 0; j < n2; j++, k+=n1) {
			rfftb_(&n, &a[k], wsave1);
		}

	}
	return;
}

void	filter (z, kx, ky, nx, ny, fparams)
float	*z, *kx, *ky;
int	nx, ny;
struct FPARAMS *fparams;
{
	int	i, j, k;
	double	pow_flex_lam_4, two_pi_c, four_pi_d, two_pi_d, gauss_factor, gf;
	double	f, dky2, dkr2, dkr4, dkr, w2, w1, phi, bf, hf, qh1, qh2; 
	double	cbspl_basis_fn (double x);
	double  log2(double x);


	two_pi_d = 2.0 * M_PI * fparams->depth;
	four_pi_d = 2.0 * two_pi_d;
	two_pi_c = 2.0 * M_PI * fparams->crust;
	pow_flex_lam_4 = pow(fparams->flex_lambda, 4.0);
	f = 5.34 * fparams->gauss_sigma;
	gauss_factor = -M_LN2 * f * f;
	qh1 = fparams->n1 / 40000.0;
	qh2 = fparams->n2 / 40000.0;	/* Old routine uses simple 40k scaling  */

	for (j = 0, k = 0; j < ny; j++) {
		dky2 = (double)ky[j];
		for (i = 0; i < nx; i++, k++) {
			dkr2 = dky2 + (double)kx[i];
			dkr4 = dkr2 * dkr2;
			dkr = sqrt(dkr2);
			f = (fparams->depth == 0.0) ? 1.0 : exp(two_pi_d * dkr);	/* Downward continue.  */
			if (fparams->do_w2) {
				w2 = 1.0 / (1.0 + fparams->w2a*dkr4*exp(dkr*four_pi_d) );
				f *= w2;
			}
			if (fparams->do_flex) {
				phi = 1.0 / (1.0 + pow_flex_lam_4 * dkr4);
				w1 = (fparams->do_flex < 0) ? phi * exp(-two_pi_c * dkr) : 1.0 - phi * exp(-two_pi_c * dkr);
				f *= w1;
			}
			if (fparams->knot) {
				bf = (dkr > 0.0) ? cbspl_basis_fn (log2(dkr) + 9.0 - fparams->knot) : 0.0;
				f *= bf;
 				f *= 1.45;  /* Try a 14.5 m/mGal scale factor for now; meaning rho_c = 2.67 */
			}
			if (fparams->do_gauss) {
				gf = (fparams->do_gauss < 0) ? exp(gauss_factor * dkr2) : 1.0 - exp(gauss_factor * dkr2);
				f *= gf;
			}
			if (fparams->do_harm) {
				hf = 0.0;
				if (fparams->do_harm > 0) {
					/* high pass */
					if (dkr >= qh2) {
						hf = 1.0;
					}
					else if (dkr > qh1) {
						hf = cos(M_PI_2 * (dkr-qh1)/(qh2-qh1));
						hf = 1.0 - hf * hf;
					}
				} else {
					/* low pass */
					if (dkr <= qh1) {
						hf = 1.0;
					}
					else if (dkr < qh2) {
						hf = cos(M_PI_2 * (dkr-qh1)/(qh2-qh1));
						hf = hf * hf;
					}
				}
				f *= hf;
			}
			z[k] = (float) (z[k] * f);
		}
	}
	return;
}

double  cbspl_basis_fn (double x) {

	/* Given x, return the value of the cubic B-spline
	basis function B(x) on equally spaced knots of width
	1 unit of x and centered at x = 0.
	
	B(x), has the following properties:
	B(x) is piecewise cubic and has continuous 2nd
	derivatives.
	B(x) has compact support; B(x) is identically
	zero for |x| >= 2 and B' and B'' are zero
	at |x| = 2.
	Integral from -oo to +oo of B(x) equals 1.
	B(0) = 2/3; B(-1) = B(1) = 1/6.
	B'(0) = 0; B'(1) = -B'(-1) = -1/2.
	B''(0) = -2; B''(-1) = B''(1) = 1.
	
	Using Bracewell's definition of equivalent
	width, the width of B(x) is 3/2.  It turns
	out that B(x) is a very good approximation
	of a Gaussian function, but with the 
	advantage of having finite support.  
	The gaussian with the same equivalent width 
	would be exp[-pi * (2*x/3)^2].
	
	This function can be used to build a cubic spline
	on equally spaced knots.  One needs a linear 
	combination of these B(x) each translated by 1
	unit of x along the x axis.  In this combination
	there will be only four translated B functions
	which are non-zero at any x.  */
	
	
	double  t, a, a2;
	
	t = fabs(x);
	
	if (t >= 2.0) return (0.0);
	
	if (t > 1.0) {
		a = 2 - t;
		a2 = a * a;
		return (a2 * a / 6.0);
	}
	
	a2 = t * t;
	return ((2.0/3.0) + a2*(0.5*t - 1));
}