#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#include "ajusteParab.h"

// Sum ( yi - a * xi * xi - b * xi - c ) ^ 2
//   minimo:    d/da, d/db, d/dc  =  0
// 
// (sumatorios)
//   /  xi^4  xi^3  xi^2 \  / a \     / y * x^2 ) 
//   |  xi^3  xi^2  xi   |  | b |  =  | y * x   )
//   \  xi^2  xi    1    /  \ c /     \ y       )
//   ... a =    b =    c =

void acerosdAjusteParab (tajub ajub)
{
  ajub->n    = 0;    ajub->sx   = 0.0;  ajub->sx2  = 0.0;
  ajub->sx3  = 0.0;  ajub->sx4  = 0.0;
  ajub->sy   = 0.0;  ajub->syx  = 0.0;  ajub->syx2 = 0.0;
  ajub->a    = 0.0;  ajub->b    = 0.0;  ajub->c    = 0.0;
}

void nuevodAjusteParab (tajub * najub)
{
  tajub ajub;
  ajub = (tajub) malloc (sizeof(struct teajub));
  acerosdAjusteParab (ajub);
  *najub = ajub ;
}

void borradAjusteParab (tajub ajub)
{
  free (ajub);
}

void imprimedAjusteParab (tajub ajub)
{
  printf ("n:    %5d\n", ajub->n);
  printf ("sx:   %9.3f\n", ajub->sx);
  printf ("sx2:  %9.3f\n", ajub->sx2);
  printf ("sx3:  %9.3f\n", ajub->sx3);
  printf ("sx4:  %9.3f\n", ajub->sx4);
  printf ("sy:   %9.3f\n", ajub->sy);
  printf ("syx:  %9.3f\n", ajub->syx);
  printf ("syx2: %9.3f\n", ajub->syx2);
  printf ("a:    %9.4f\n", ajub->a);
  printf ("b:    %9.4f\n", ajub->b);
  printf ("c:    %9.4f\n", ajub->c);
  printf ("\n");
}

void anyadedAAjusteParab (tajub ajub, int x, int y)
{
  ajub->n    += 1 ;
  ajub->sx   += x ;         ajub->sx2  += x * x ;
  ajub->sx3  += x * x * x;  ajub->sx4  += x * x * x * x;
  ajub->sy   += y ;        
  ajub->syx  += y * x ;     ajub->syx2 += y * x * x;  
}

void calculaAjusteParab (tajub ajub)
{
  double del;
  del  =        ajub->sx4 * ajub->sx2 * ajub->n
	  + 2 * ajub->sx3 * ajub->sx2 * ajub->sx
	  -     ajub->sx2 * ajub->sx2 * ajub->sx2
	  -     ajub->sx4 * ajub->sx  * ajub->sx 
	  -     ajub->sx3 * ajub->sx3 * ajub->n ;
  // printf ("del: %9.2f\n", del);
  ajub->a =   (   ajub->syx2 * (ajub->sx2 * ajub->n   - ajub->sx  * ajub->sx )
	        + ajub->syx  * (ajub->sx  * ajub->sx2 - ajub->sx3 * ajub->n  )
		+ ajub->sy   * (ajub->sx3 * ajub->sx  - ajub->sx2 * ajub->sx2)
	      ) / del ;

  ajub->b =   (   ajub->syx2 * (ajub->sx2 * ajub->sx  - ajub->sx3 * ajub->n  )
	        + ajub->syx  * (ajub->sx4 * ajub->n   - ajub->sx2 * ajub->sx2)
		+ ajub->sy   * (ajub->sx3 * ajub->sx2 - ajub->sx4 * ajub->sx )
	      ) / del ;

  ajub->c =   (   ajub->syx2 * (ajub->sx3 * ajub->sx  - ajub->sx2 * ajub->sx2)
	        + ajub->syx  * (ajub->sx2 * ajub->sx3 - ajub->sx4 * ajub->sx )
		+ ajub->sy   * (ajub->sx4 * ajub->sx2 - ajub->sx3 * ajub->sx3)
	      ) / del ;
}

double distanciaAParab (tajub ajub, int x, int y)
{
  return (y - (ajub->a*x*x + ajub->b*x + ajub->c)) ;
}