/*
 *	Karinto Library Project
 *
 *	This software is distributed under a zlib-style license.
 *	See license.txt for more information.
 */

using System;
using System.Collections.Generic;
using System.Text;

namespace Karinto
{
    /// <summary>
    ///    3次スプライン
    /// </summary>
    public class Spline
    {
        private int n;
        private double[] x;
        private double[] a; //0次係数
        private double[] b; //1次係数
        private double[] c; //2次係数
        private double[] d; //3次係数
        public int N
        {
            get { return n; }
        }

        public double[] XArray
        {
            get { return x; }
        }
        public double[] YArray
        {
            get { return a; }
        }

        /// <summary>
        ///     自然3次スプライン(両端の2階微分がゼロ)
        /// </summary>
        /// <param name="list">点列</param>
        public Spline(PointList list)
            : this(list, Double.NaN, Double.NaN)
        {
        }


        /// <summary>
        ///     締3次スプライン
        /// </summary>
        /// <param name="list">点列</param>
        /// <param name="firstGradient">最初のxにおける1階微分(NaNで自然境界)</param>
        /// <param name="lastGradient">最後のxにおける1階微分(NaNで自然境界)</param>
        public Spline(PointList list, double firstGradient, double lastGradient)
        {
            x = list.XArray;
            a = list.YArray;
            n = list.Count;
            if (n < 3) throw new ArgumentException("The point list requires more than 2 points.");
            b = new double[n];
            c = new double[n];
            d = new double[n];

            double[] h = new double[n - 1];
            double[] z = new double[n]; //A * c
            double[] m = new double[n]; //右側
            double l;
            for (int i = 0; i < n - 1; i++)
            {
                h[i] = x[i + 1] - x[i];
            }

            if (Double.IsNaN(firstGradient))
            {
                z[0] = 0;
                m[0] = 0;
            }
            else
            {
                z[0] = (3 * (a[1] - a[0]) / h[0] - 3 * firstGradient) / (2 * h[0]);
                m[0] = 0.5;
            }
            for (int i = 1; i < n - 1; i++)
            {
                z[i] = 3 * (a[i + 1] - a[i]) / h[i] - 3 * (a[i] - a[i - 1]) / h[i - 1];
            }

            for (int i = 1; i < n - 1; i++) //左下をゼロにする
            {
                l = 2 * (h[i] + h[i - 1]) - h[i - 1] * m[i - 1]; // 対角成分
                // lで割ることで対角成分を1にする
                m[i] = h[i] / l;
                z[i] = (z[i] - h[i - 1] * z[i - 1]) / l;
            }
            if (Double.IsNaN(lastGradient))
            {
                c[n - 1] = 0;
            }
            else
            {
                z[n - 1] = 3 * lastGradient - 3 * (a[n - 1] - a[n - 2]) / h[n - 2];
                l = 2 * h[n - 2] - h[n - 2] * m[n - 2]; // 対角成分
                c[n - 1] = (z[n - 1] - h[n - 2] * z[n - 2]) / l;
            }
            for (int i = n - 2; i >= 0; i--)
            {
                c[i] = z[i] - m[i] * c[i + 1];
                b[i] = (a[i + 1] - a[i]) / h[i] - h[i] * (c[i + 1] + 2 * c[i]) / 3;
                d[i] = (c[i + 1] - c[i]) / (3 * h[i]);
            }
            b[n - 1] = b[n - 2] + (2 * c[n - 2] + 3 * d[n - 2] * h[n - 2]) * h[n - 2];
        }

        public override string ToString()
        {
            string s = "i\tx\ta\tb\tc\td\r\n";
            int i;
            for (i = 0; i < x.Length - 1; i++)
            {
                s += "" + i + "\t" + x[i] + "\t" + a[i] + "\t" + b[i] + "\t" + c[i] + "\t" + d[i] + "\r\n";
            }
            s += "" + i + "\t" + x[i] + "\t" + a[i] + "\t\t\t\r\n";
            return s;
        }


        public double GetY(double t)
        {
            int i;
            if (x[0] > t) return a[0] - b[0] * (x[0] - t);//線形補外
            if (x[n - 1] < t) return a[n - 1] + b[n - 1] * (t - x[n - 1]);//線形補外
            for (i = 0; i < n - 2; i++)
            {
                if (x[i + 1] >= t) break;
            }
            double h = t - x[i];
            return a[i] + (b[i] + (c[i] + d[i] * h) * h) * h;
        }

        public double[] GetCoefficients(int order)
        {
            switch (order)
            {
                case 0:
                    return a;
                case 1:
                    return b;
                case 2:
                    return c;
                case 3:
                    return d;
            }
            throw new ArgumentException("Invalid order");
        }

        public PointList GetFrozenList(double step)
        {
            return GetFrozenList(new Range(x[0], x[n - 1]), step);
        }

        public PointList GetFrozenList(Range range, double step)
        {
            PointList list = new PointList();
            if (Double.Epsilon > step) throw new ArgumentException("step is too small.");
            foreach (double t in range.Step(step))
            {
                list.Add(new Point(t, GetY(t)));
            }
            return list;
        }
    }
}