package visad.util;

import Jama.Matrix;
import Jama.LUDecomposition;
import java.util.Arrays;

public class CubicInterpolator {

      LUDecomposition solver;

      double[][] solution = null;

      double x0 = 0;
      double x1 = 0;
      double x2 = 0;
      double x0_last = 0;
      double x0_save;
      
      float[] values0 = null;
      float[] values1 = null;
      float[] values2 = null;
      float[] values0_last = null;
      float[] values0_save = null;

      int numSpatialPts = 1;

      boolean doIntrp = true;
      
      boolean[] needed = null;
      boolean[] computed = null;
      
      public CubicInterpolator(boolean doIntrp, int numSpatialPts) {
         this.doIntrp = doIntrp;
         this.numSpatialPts = numSpatialPts;
         this.solution = new double[4][numSpatialPts];
         this.needed = new boolean[numSpatialPts];
         this.computed = new boolean[numSpatialPts];
         Arrays.fill(needed, false);
         Arrays.fill(computed, false);
      }

      private void buildSolver() {
         double x0_p3 = x0*x0*x0;
         double x1_p3 = x1*x1*x1;

         double x0_p2 = x0*x0;
         double x1_p2 = x1*x1;

         Matrix coeffs = new Matrix(new double[][]
              { {x0_p3, x0_p2, x0, 1},
                {x1_p3, x1_p2, x1, 1},
                {3*x0_p2, 2*x0, 1, 0},
                {3*x1_p2, 2*x1, 1, 0}}, 4, 4);

         solver = new LUDecomposition(coeffs);         
      }

      public void interpolate(double xt, float[] interpValues) {
         if (!doIntrp) {
            if (xt == x0) {
               System.arraycopy(values0, 0, interpValues, 0, numSpatialPts);
            }
            else if (xt == x1) {
               System.arraycopy(values1, 0, interpValues, 0, numSpatialPts);
            }
            return;
         }
         java.util.Arrays.fill(interpValues, Float.NaN);
         
         for (int k=0; k