package algorithmrepository;
/**
 *
 * Usage:
 * The Gauss-Legendre method for integration is based on finding
 * an approximation of the form
 *
 *      I = sum(w(i)*f(p(i))
 *
 * where w is a weight vector and p (poles) a vector of positions at
 * which the function should be evaluated. The weights and poles is
 * chosen in such a way as to make the approximation *exact* for any
 * polynomial of a degree < = 2n-1, where n is the size of the weight
 * and poles vectors (the order). The default integration interval
 * is [-1,1]. If another integration interval is chosen, the generated
 * weigths and poles have to be transformed, done in the setLimits()
 * method.
 *
 * Directly translated from Andreas Werners gaussLegendre.cpp
 *
 * See http://www.algorithmrepository.org
 *
 */
/** @deprecated Use static methods in algorithmrepository.Algorithms instead. */ 
public class GaussLegendre {

    static final int PRECISION = 15;
    
    int n;  
    double[] poles, weights, x, w;
    
        
    /**
     * Default constructor.
     */
    public GaussLegendre() {
      n = -1;
      poles = null;
      weights = null;
      x = null;
      w = null;
    }

    
    /**
     * Construct a Gauss-Legendre integrator of
     * the given order. Default limits are -1 to +1.
     * @param order The number of terms in the sum approximation
     * of the integral.
     */
    public GaussLegendre(int order) {
      n = order;
      poles = null;
      weights = null;
      x = null;
      w = null;
      initWeights();
    }
    
    /**
     * Constructs a Gauss-Legendre integrator for a given
     * order and lower, upper bounds of the integral.
     *
     * @param order Number of terms in the sum approximation
     * of the integral.
     * @param lower The lower integration limit.
     * @param upper The upper integration limit.
     */
    public GaussLegendre(int order, double lower, double upper) {
        this(order);
        setLimits(lower, upper);
    }

    
    /**
     * Convenience method to make a table of poles and weights
     * for given order and lower,upper limits.
     *
     * @param order The number of terms in the integral approximation.
     * @param minlim The lower bound of integration.
     * @param maxlim The upper bound of integration.
     */
    public static void printWeightsAndPoles(int order, double minlim, double maxlim) {
        GaussLegendre gl = new GaussLegendre(order);
        gl.setLimits(minlim,maxlim);
        
        double[] w = gl.getWeights();
        double[] p = gl.getPoles();
        
        System.out.println("Weights and Poles for order "+order+", interval "+minlim+"->"+maxlim);
        for(int i=0;i<w.length;++i) {
            System.out.println(i+": "+w[i]+", "+p[i]);
        }
        
    }
    
    /**
     * Sets the order of this integrator. This is
     * the number of terms in the sum that approximates
     * the integral.
     *
     * @param order The order of the Gauss-Legendre evaluation.
     */
    public void setOrder(int order)
    {
      if(n == order) return;

      n = order;
      poles = null;
      weights = null;
      x = null;
      w = null;
      initWeights();
    }

    /**
     * Changes the limits of the integral.
     *
     * @param l Lower limit of integral.
     * @param u Upper limit of integral.
     */
    public void setLimits(double l, double u) {
      double d = (u - l) / 2.0;

      x = new double[n];
      w = new double[n];

      for(int i = 0; i < n; i++) {
        x[i] = d * (poles[i] + 1.0) + l;
        w[i] = d * weights[i];
      }
    }
    
    /**
     * Returns the positions at which
     * the integrand should be evaluated.
     *
     * @return The positions at which the
     * integrand should be evaluated.
     */
    public double[] getPoles() {
      return x;
    }
    
    /**
     * Returns the weights with which the
     * integrand evaluated at the corresponding
     * poles should be multiplied.
     *
     * @return The weights with which the integrand
     * should be multiplied at the pole positions.
     */
    public double[] getWeights() {
      return w;
    }
    
    private void initWeights() {
      poles = new double[n];
      weights = new double[n];

            
     //Start calculation of poles and weights, C++ version first tries to read in file with coeffs     
      double[][] ja = new double[n][];
      double[][] jv = new double[n][];
      for(int i = 0; i < n; i++) ja[i] = new double[i + 1];
      for(int k = 0; k < n; k++) jv[k] = new double[n];
     
     for(int k = 0; k < n - 1; k++) {
       poles[k] = (k + 1) / Math.sqrt(4.0 * (k + 1) * (k + 1) - 1.0);
     }
     for(int k = 0; k < n; k++) {
       for(int i = 0; i <= k; i++) {
         if(i == k - 1)
           ja[k][i] = poles[i];
         else
           ja[k][i] = 0.0;
       }
     }
     jacobi(ja, jv);
     for(int k = 0; k < n; k++) {
       weights[k] = 2.0 * jv[0][k] * jv[0][k];
       poles[k] = ja[k][k];
     }

    }

   

    private void jacobi(double[][] ja, double[][] jv) {
      double delta = Math.pow(10.0, (-PRECISION));  // static in C++ version
      double eps = delta; // static in C++ version

      for(int i = 0; i < n; i++) {
        jv[i][i] = 1;
        for(int j = 0; j < n; j++)
          if(i != j)
            jv[i][j] = 0;
      }

      int k = 0;
      double residuum = Double.MAX_VALUE;
      double sum = 0;

      do {
        for(int p = 0; p < n - 1; p++) {
          for(int q = p + 1; q < n; q++) {
            if(Math.abs(ja[q][p]) > eps) {
              double teta = (ja[q][q] - ja[p][p]) / (2.0 * ja[q][p]);
              double t = 1.0;

              if(Math.abs(teta) > delta) {
                if(teta > 0.0)
                  t = 1.0 / (teta + Math.sqrt(teta * teta + 1.0));
                else
                  t = 1.0 / (teta - Math.sqrt(teta * teta + 1.0));
              }
              double c = 1.0 / Math.sqrt(1.0 + t * t);
              double s = c * t;
              double r = s / (1.0 + c);

              ja[p][p] -= t * ja[q][p];
              ja[q][q] += t * ja[q][p];
              ja[q][p] = 0;
              for(int j = 0; j < p; j++) {
                double g = ja[q][j] + r * ja[p][j];
                double h = ja[p][j] - r * ja[q][j];

                ja[p][j] -= s * g;
                ja[q][j] += s * h;
              }
              for(int i = p + 1; i < q; i++) {
                double g = ja[q][i] + r * ja[i][p];
                double h = ja[i][p] - r * ja[q][i];

                ja[i][p] -= s * g;
                ja[q][i] += s * h;
              }
              for(int i = q + 1; i < n; i++) {
                double g = ja[i][q] + r * ja[i][p];
                double h = ja[i][p] - r * ja[i][q];

                ja[i][p] -= s * g;
                ja[i][q] += s * h;
              }
              for(int i = 0; i < n; i++) {
                double g = jv[i][q] + r * jv[i][p];
                double h = jv[i][p] - r * jv[i][q];

                jv[i][p] -= s * g;
                jv[i][q] += s * h;
              }
            }
          }
        }
        for(int i = 1; i < n; i++) {
          for(int j = 0; j < i; j++) {
            sum += ja[i][j] * ja[i][j];
          }
        }
        if(Math.abs(residuum - (2.0 * sum - eps * eps)) < eps)
          k++;
        residuum = 2.0 * sum - eps * eps;
      } while((2.0 * sum) > (eps * eps) && k < 5);
    }
    
}