package jp.hasc.hasctool.core.data;

//* スプライン補間(１次元) */
public class Spline {

	/** 補助データの作成 */
	private double[] x, y, diff2;
	private int N;

	// コンストラクタ
	// x.length == y.length であること
	// x[0] < x[1] < ... < x[N-1] であること
	public Spline(double[] x, double[] y) {
		this.x = x;
		this.y = y;
		N = x.length;
		diff2 = new double[N];
		maketable();
	}

	private void maketable() {
		double[] h = new double[N];
		double[] diff1 = new double[N];

		// 両端点での y''(x) / 6 ... (自然スプライン)
		diff2[0] = 0;  diff2[N - 1] = 0;
		for(int i=0; i<N-1; i++) {
			h[i    ] =  x[i + 1] - x[i];
			diff1[i + 1] = (y[i + 1] - y[i]) / h[i];
		}
		diff2[1] = diff1[2] - diff1[1] - h[0] * diff2[0];
		diff1[1] = 2 * (x[2] - x[0]);
		for(int i = 1;i < N - 2;i++) {
			double t = h[i] / diff1[i];
			diff2[i + 1] = diff1[i + 2] - diff1[i + 1] - diff2[i] * t;
			diff1[i + 1] = 2 * (x[i + 2] - x[i]) - h[i] * t;
		}
		diff2[N - 2] -= h[N - 2] * diff2[N - 1];
		for(int i=N-2; i>0; i--) {
			diff2[i] = (diff2[i] - h[i] * diff2[i + 1]) / diff1[i];
		}
	}

	public double map(double t) {
		int i, j, k;
		double d, h;

		i = 0;  j = N - 1;
		while (i < j) {
			k = (i + j) / 2;
			if (x[k] < t) i = k + 1;  
			else j = k;
		}
		if (i > 0) i--;
		h = x[i + 1] - x[i];  
		d = t - x[i];
		return (((diff2[i + 1] - diff2[i]) * d / h + diff2[i] * 3) * d + ((y[i + 1] - y[i]) / h - (diff2[i] * 2 + diff2[i + 1]) * h)) * d + y[i];
	}

	public double getSpline(double t){
		int i;
		double temp_x,temp_y = 0,yy0,yy1,yy2,yy3;
		double[] h=new double[N];     //間隔
		double[] dif1=new double[N];  //一次微分
		double[] dif2=new double[N];  //二次微分

		h[0]=0.0;
		dif2[0]=0.0;
		dif2[N-1]=0.0; 

		for(i=1;i<N;i++){
			h[i]=x[i]-x[i-1];             //間隔を計算
			dif1[i]=(y[i]-y[i-1])/h[i];   //一次微分を計算
		}

		for(i=1;i<N - 1;i++)
			//二次微分を計算
			dif2[i]=(dif1[i+1]-dif1[i])/(x[i+1]-x[i-1]);

		i=1;
		for(temp_x=x[0];temp_x<x[N-1];temp_x+=0.01){
			if(temp_x<x[i]){
				yy0=dif2[i-1]/(6*h[i])*(x[i]-temp_x)*(x[i]-temp_x)*(x[i]-temp_x);       //第１項
				yy1=dif2[i]/(6*h[i])*(temp_x-x[i-1])*(temp_x-x[i-1])*(temp_x-x[i-1]);   //第２項
				yy2=(y[i-1]/h[i]-h[i]*dif2[i-1]/6)*(x[i]-temp_x);                    //第３項
				yy3=(y[i]/h[i]-h[i]*dif2[i]/6)*(temp_x-x[i-1]);                   //第４項
				temp_y=yy0+yy1+yy2+yy3;
				return temp_y;
			}
			else i++;
		}
		return temp_y;
	}
}
