//
//  EXPSpline.m
//  Spline-Test
//
//  Created by ashley on 04/06/2008.
//  Copyright 2008 __MyCompanyName__. All rights reserved.
//

#import "EXPSpline.h"

@implementation EXPSpline

-(id)init
{
	self = [super init];

	if (self!=nil) {
		_x = NULL;
		_y = NULL;
		_b = NULL;
		_c = NULL;
		_d = NULL;
		_isRegular = NO;
		_last = 0;
	}

	return self;
}

-(BOOL)setPoints:(NSString *)interpType nPoints:(int)n xPoints:(double *)x yPoints:(double *)y
{
	if (n<=2) {
		return NO;
	}
	
	_isRegular = NO;
	
	int nm1, ib, i;
	double t;
		
	if (_b!=NULL) free(_b);
	if (_c!=NULL) free(_c);
	if (_d!=NULL) free(_d);
	if (_x!=NULL) free(_x);
	if (_y!=NULL) free(_y);
	[_interpType release];
	
	_x = (double *)calloc(n, sizeof(double));
	_y = (double *)calloc(n, sizeof(double));

// Check that c[i] are monotonic?
	_x[0] = x[0];
	_y[0] = y[0];
	BOOL increasing = x[1] > x[0];
	for(i=1; i<n; i++) {
		_x[i] = x[i];
		_y[i] = y[i];
		if ((x[i]>x[i-1])!=increasing) {
			[self release];
			return nil;
		}
	}
		
	_n = n;
	_interpType = [interpType retain];
	if (![interpType isEqualToString:@"spline"]) {
		return YES;
	}
		
	_b = (double *)calloc(n, sizeof(double));
	_c = (double *)calloc(n, sizeof(double));
	_d = (double *)calloc(n, sizeof(double));
		
	nm1 = n - 1;
	if (n<3) {
		_b[0] = (y[1] - y[0])/(x[1] - x[0]);
		_c[0] = 0.0;
		_d[0] = 0.0;
		_b[1] = _b[0];
		_c[1] = 0.0;
		_d[1] = 0.0;
		return YES;
	}
/*
	set up tridiagonal system

	b = diagonal, d = offdiagonal, c = right hand side.
*/
	_d[0] = x[1] - x[0];
	_c[1] = (y[1] - y[0])/_d[0];
	for(i=1; i<nm1; i++) {
		_d[i] = x[i+1] - x[i];
		_b[i] = 2.0*(_d[i-1] + _d[i]);
		_c[i+1] = (y[i+1] - y[i])/_d[i];
		_c[i] = _c[i+1] - _c[i];
//			printf("%d %g %g  %g\n", i, _b[i], _c[i], _d[i]);
	}
/*
	end conditions.  third derivatives at  x(1)  and  x(n)
	obtained from divided differences
*/
	_b[0] = -_d[0];
	_b[n-1] = -_d[n-2];
	_c[0] = 0.0;
	_c[n-1] = 0.0;
	if (n!=3) {
		_c[0] = _c[2]/(x[3] - x[1]) - _c[1]/(x[2] - x[0]);
		_c[n-1] = _c[n-2]/(x[n-1] - x[n-3]) - _c[n-3]/(x[n-2] - x[n-4]);
		_c[0] = _c[0]*pow(_d[0], 2)/(x[3] - x[0]);
		_c[n-1] = -_c[n-1]*pow(_d[n-2],2)/(x[n-1] - x[n-4]);
	}
/*
	forward elimination
*/
	for(i = 1; i<n; i++) {
		t = _d[i-1]/_b[i-1];
		_b[i] = _b[i] - t*_d[i-1];
		_c[i] = _c[i] - t*_c[i-1];
//			printf("%d %g %g  %g\n", i, t, _b[i], _c[i]);
	}
/*
	back substitution
*/
	_c[n-1] = _c[n-1]/_b[n-1];
	for(ib=0; ib<nm1; ib++) {
		i = n - ib - 2;
		_c[i] = (_c[i] - _d[i]*_c[i+1])/_b[i];
//			printf("%d %g  %g\n", i, _c[i], _c[i+1]);
	}
/*
	c(i) is now the sigma(i) of the text

	compute polynomial coefficients
*/
	_b[n-1] = (y[n-1] - y[n-2])/_d[n-2] + _d[n-2]*(_c[n-2] + 2.0*_c[n-1]);
	for(i=0; i<nm1; i++) {
		_b[i] = (y[i+1] - y[i])/_d[i] - _d[i]*(_c[i+1] + 2.0*_c[i]);
		_d[i] = (_c[i+1] - _c[i])/_d[i];
		_c[i] = 3.0*_c[i];
//			printf("%d %g %g  %g\n", i, _b[i], _c[i], _d[i]);
	}
	_c[n-1] = 3.0*_c[n-1];
	_d[n-1] = _d[n-2];
	return YES;
}

-(BOOL)setPoints:(NSString *)interpType nPoints:(int)n xMin:(double)xMin xMax:(double)xMax yPoints:(double *)y
{
	_isRegular = YES;
	if (n<=2) {
		return NO;
	}
	
	if (_b!=NULL) free(_b);
	if (_c!=NULL) free(_c);
	if (_d!=NULL) free(_d);
	if (_x!=NULL) free(_x);
	if (_y!=NULL) free(_y);
	[_interpType release];
	
	_x = NULL;
	_y = (double *)calloc(n, sizeof(double));
	_b = (double *)calloc(n, sizeof(double));
	_c = (double *)calloc(n, sizeof(double));
	_d = (double *)calloc(n, sizeof(double));
		
	int i;
	double h = (xMax - xMin)/(n - 1);
	for(i=0; i<n; i++) {
		_y[i] = y[i];
	}
	
	_xmin = xMin;
	_xmax = xMax;
	_n = n;
	
	_interpType = [interpType retain];
	if (![interpType isEqualToString:@"spline"]) {
		return YES;
	}

	int nm1 = n - 1;
	if (n<3) {
		_b[0] = (y[1] - y[0])/h;
		_c[0] = 0.0;
		_d[0] = 0.0;
		_b[1] = _b[0];
		_c[1] = 0.0;
		_d[1] = 0.0;
		return YES;
	}
/*
	set up tridiagonal system

	b = diagonal, d = offdiagonal, c = right hand side.
*/
	_d[0] = h;
	_c[1] = (y[1] - y[0])/_d[0];
	for(i=1; i<nm1; i++) {
		_d[i] = h;
		_b[i] = 2.0*(_d[i-1] + _d[i]);
		_c[i+1] = (y[i+1] - y[i])/_d[i];
		_c[i] = _c[i+1] - _c[i];
//			printf("%d %g %g  %g\n", i, _b[i], _c[i], _d[i]);
	}
/*
	end conditions.  third derivatives at  x(1)  and  x(n)
	obtained from divided differences
*/
	_b[0] = -_d[0];
	_b[n-1] = -_d[n-2];
	_c[0] = 0.0;
	_c[n-1] = 0.0;
	if (n!=3) {
		_c[0] = (_c[2] - _c[1])/(2.0*h);
		_c[n-1] = (_c[n-2] - _c[n-3])/(2.0*h);
		_c[0] = _c[0]*pow(_d[0], 2)/(3.0*h);
		_c[n-1] = -_c[n-1]*pow(_d[n-2],2)/(3.0*h);
	}
/*
	forward elimination
*/
	for(i = 1; i<n; i++) {
		double t = _d[i-1]/_b[i-1];
		_b[i] = _b[i] - t*_d[i-1];
		_c[i] = _c[i] - t*_c[i-1];
//			printf("%d %g %g  %g\n", i, t, _b[i], _c[i]);
	}
/*
	back substitution
*/
	_c[n-1] = _c[n-1]/_b[n-1];
	int ib;
	for(ib=0; ib<nm1; ib++) {
		i = n - ib - 2;
		_c[i] = (_c[i] - _d[i]*_c[i+1])/_b[i];
//			printf("%d %g  %g\n", i, _c[i], _c[i+1]);
	}
/*
	c(i) is now the sigma(i) of the text

	compute polynomial coefficients
*/
	_b[n-1] = (y[n-1] - y[n-2])/_d[n-2] + _d[n-2]*(_c[n-2] + 2.0*_c[n-1]);
	for(i=0; i<nm1; i++) {
		_b[i] = (y[i+1] - y[i])/_d[i] - _d[i]*(_c[i+1] + 2.0*_c[i]);
		_d[i] = (_c[i+1] - _c[i])/_d[i];
		_c[i] = 3.0*_c[i];
//			printf("%d %g %g  %g\n", i, _b[i], _c[i], _d[i]);
	}
	_c[n-1] = 3.0*_c[n-1];
	_d[n-1] = _d[n-2];
	
	return YES;
}

-(double)evaluate:(double)u primes:(int)nPrimes
{
/*
	this subroutine evaluates the cubic spline function

	seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3

	where  x(i) .lt. u .lt. x(i+1), using horner's rule

	f  u .lt. x(1) then  i = 1  is used.
	if  u .ge. x(n) then  i = n  is used.

	input..

	n = the number of data points
	u = the abscissa at which the spline is to be evaluated
	x,y = the arrays of data abscissas and ordinates
	b,c,d = arrays of spline coefficients computed by spline

	if  u  is not in the same interval as the previous call, then a
	binary search is performed to determine the proper interval.
*/
	int i = _last;
	double dx;
	double h;
	if (i>=_n-1) i = 0;
/*
	binary search
*/
	if (![self isRegular]) {
		if (u<_x[0]) {
			return _y[0];
		} else if (u>_x[_n-1]) {
			return _y[_n-1];
		} 
		int j, k;
		if((u>_x[i+1]) || (u<_x[i])) {
			i = 0;
			j = _n;
			do {
				k = (i+j)/2;
				if (u<_x[k]) j = k;
				if (u>=_x[k]) i = k;
			} while (j>(i+1));
		}
		dx = u - _x[i];
		h = _x[i+1] - _x[i];
	} else {
		if (u<_xmin) {
			return _y[0];
		} else if (u>_xmax) {
			return _y[_n-1];
		} 
		h = (_xmax - _xmin)/(_n - 1);
		i = floor(u/h);
		dx = u - i*h;
	}

	double seval = 0.0;
	NSString *interpType = [self interpType];
/*
	evaluate spline
*/
	if ([interpType isEqualToString:@"stepwise"]) {
		if (nPrimes==0) {
			seval = _y[i];
		}
	} else if ([interpType isEqualToString:@"linear"]) {
		if (nPrimes==0) {
			seval = _y[i] + dx*(_y[i+1] - _y[i])/h;
		} else if (nPrimes==1) {
			seval = (_y[i+1] - _y[i])/h;
		}
	} else if ([interpType isEqualToString:@"spline"]) {
		if (nPrimes==0) {
			seval = _y[i] + dx*(_b[i] + dx*(_c[i] + dx*_d[i]));
		} else if (nPrimes==1) {
		} else if (nPrimes==2) {
		} else if (nPrimes==3) {
		}
	} else {
	}
	
	return seval;
}

- (NSString *) description
{
	NSMutableString *desc = [[NSMutableString alloc] init];
	[desc appendString:@"EXPSpline:\n"];
	int i;
	for(i=0; i<_n; i++) {
		NSString *str = [[NSString alloc] initWithFormat:@"%d %15.6f %15.6f %15.6f\n", i, _b[i], _c[i], _d[i]];
		[desc appendString:str];
		[str release];
	}
	
	return desc;
}

-(NSString *)interpType
{
	return _interpType;
}

- (BOOL)isRegular
{
	return _isRegular;
}

- (void) dealloc
{
	if (_b!=NULL) free(_b);
	if (_c!=NULL) free(_c);
	if (_d!=NULL) free(_d);
	if (_x!=NULL) free(_x);
	if (_y!=NULL) free(_y);
	[_interpType release];
	[super dealloc];
}

@end