import java.awt.*;
import java.util.*;
import java.applet.Applet;

public abstract class LinearFit extends Pplotter { 

  int Ndata;   /* number of data values  */
     /*   the number of parameters is assumed to be 2  */
  double[] X0;    /* vector of reference values of data */
  double[][] A;  /*   matrix expressing linear dependence of X_i on p_k  */
  double[] data;  /* vector of data */  
  double[] error; /* vector of errors */
  double[] center; /* central value of fit */
  double[] centererr; /* standard error on each parameter p_k */
  double[]  sigma;  
           /* sigma_a for each principal direction in parameter space */
  double[][] vp;     /*   vp[k][a]  is the kth component 
                            of the ath principal (unit) vector */
  double chisq;   /* chisquared at the minimum    */
  int dof;        /* number of data points - number of fit parameters  */

  public LinearFit(int NDATA){
    Ndata = NDATA;
    X0 = new double[Ndata];
    fillX0();
    A = new double[Ndata][2];
    data = new double[Ndata];
    error = new double[Ndata];
    center = new double[2];
    centererr = new double[2];
    sigma = new double[2];
    vp = new double[2][2];
    dof = Ndata -2;
  }

  abstract void fillX0();

       /* use this method to initialize the matrix A */
  void setDependence(int i, double B0, double B1){
    A[i][0] = B0;
    A[i][1] = B1;
  }
      /* use this method to add data points, numbered from 0 to Ndata-1  */
  void addData(int i, double x, double err){
    data[i] = x;
    error[i] = err;
  }

  void fit(){
    int i, j, k;
    double[] datasiginverse = new double[Ndata];
    double[] Xdiff = new double[Ndata];
    for (i = 0; i < Ndata; i++){
      Xdiff[i] = data[i] - X0[i];
      datasiginverse[i] = 1.0/(error[i] * error[i]); 
    }
    double[][] Sigmasqinv = new double[2][2];
    double[][] Sigmasq = new double[2][2];
    for (i = 0; i < 2; i++){
      for (j = 0 ; j < 2; j++){
	for (k = 0; k < Ndata ; k++){
	  Sigmasqinv[i][j] += A[k][i] * datasiginverse[k] * A[k][j];
        }
      }
    }
    double detSigi = Sigmasqinv[0][0] *  Sigmasqinv[1][1]
                          - Sigmasqinv[0][1] *  Sigmasqinv[1][0];
    Sigmasq[0][0] = Sigmasqinv[1][1]/detSigi;
    Sigmasq[1][1] = Sigmasqinv[0][0]/detSigi;
    Sigmasq[0][1] = - Sigmasqinv[1][0]/detSigi;
    Sigmasq[1][0] = - Sigmasqinv[0][1]/detSigi;
    double[]    ASX = new double[2];
    for (j = 0 ; j < 2; j++){
	for (k = 0; k < Ndata ; k++){
	  ASX[j] += A[k][j] * datasiginverse[k] * Xdiff[k];
        }
    }
    for (i = 0; i < 2; i++){
        center[i] = 0.0;
        for (j = 0; j < 2; j++){
	  center[i] += Sigmasq[i][j] * ASX[j]; 
        }
    }
    centererr[0] = Math.sqrt(Sigmasqinv[1][1]/detSigi);
    centererr[1] = Math.sqrt(Sigmasqinv[0][0]/detSigi);
    chisq = 0.0;
    for (k = 0 ; k < Ndata; k++){
        chisq +=  datasiginverse[k] * Xdiff[k] * Xdiff[k];
    }
    for (i = 0; i < 2; i++){
        for (j = 0 ; j < 2 ; j++){ 
	    chisq -=  Sigmasq[i][j] * ASX[i] * ASX[j];
	}
    }
    double evone = (Sigmasqinv[0][0] + Sigmasqinv[1][1])/2.0;
    double evtwo = (Sigmasqinv[0][0] - Sigmasqinv[1][1])/2.0;
    double evthree = Math.sqrt( evtwo * evtwo 
				+ Sigmasqinv[0][1] *  Sigmasqinv[0][1]);
    double evplus = evone + evthree;
    double evminus = evone - evthree;
    sigma[0] = 1.0/Math.sqrt(evplus);
    sigma[1] = 1.0/Math.sqrt(evminus);
    double x2 = Sigmasqinv[0][1]/(evplus-Sigmasqinv[1][1]);
    double xn = Math.sqrt(1.0 + x2*x2);
    vp[0][0] = 1.0/xn;
    vp[1][0] = x2/xn;
    double x1 = Sigmasqinv[0][1]/(evminus-Sigmasqinv[0][0]);
    xn = Math.sqrt(1.0 + x1*x1);
    vp[0][1] = x1/xn;
    vp[1][1] = 1.0/xn;
  }

  static final double sixtyeight = 1.5096;
  static final double ninety = 2.146;   
  static final double ninetynine = 3.035; 
       /*  sigma for various 2-d confidence contours   */ 

  void plotContour(Graphics g){
      plotContour(g,center[0], center[1]);
  }

  void plotContour(Graphics g, double centerx, double centery){
      double level = sixtyeight;
      plotStream P = new plotStream();
      double theta;
      double dtheta = 0.1;
      double x,y;
      for (theta = 0.0; theta <= 2.0 * Math.PI + dtheta/2.0; theta += dtheta){
        x = centerx + level * (Math.cos(theta) * sigma[0] * vp[0][0] 
                            + Math.sin(theta) * sigma[1] * vp[0][1]);
        y = centery + level * (Math.cos(theta) * sigma[0] * vp[1][0] 
                            + Math.sin(theta) * sigma[1] * vp[1][1]);
        P.add(x,y);
      }
      setClip(g);
      plotCurve(g,P);
  }

  boolean contourOutside(){
          /* returns true if the contour falls outside the plot  */
     boolean outside = false;
     if ((center[0] - centererr[0] > xb)  || (center[0] + centererr[0] < xa)
       || (center[1] - centererr[1] > yb) || (center[1] + centererr[1] < ya)){
             outside = true;
     }
     return outside;
  }

  void plotConstraint(Graphics g, int i){
      double level = 1.0;
      double a0 = A[i][0];
      double absa0 = (a0 > 0.0) ?  a0 : -a0;
      double a1 = A[i][1];
      double absa1 = (a1 > 0.0) ?  a1 : -a1;
      double DX = data[i] - X0[i];
      double x1, y1, x2, y2; 
      setClip(g);
      if (absa0 > absa1){
         y1 = ya;
         x1 = (DX + error[i]*level - a1 * ya)/a0;
         y2 = yb;
         x2 = (DX + error[i]*level - a1 * yb)/a0;
         drawLine(g,x1,y1,x2,y2);
         y1 = ya;
         x1 = (DX - error[i]*level - a1 * ya)/a0;
         y2 = yb;
         x2 = (DX - error[i]*level - a1 * yb)/a0;
         drawLine(g,x1,y1,x2,y2);
      } else { 
         x1 = xa;
         y1 = (DX + error[i]*level - a0 * xa)/a1;
         x2 = xb;
         y2 = (DX + error[i]*level - a0 * xb)/a1;
         drawLine(g,x1,y1,x2,y2);
         x1 = xa;
         y1 = (DX - error[i]*level - a0 * xa)/a1;
         x2 = xb;
         y2 = (DX - error[i]*level - a0 * xb)/a1;
         drawLine(g,x1,y1,x2,y2);
      }
  }
}