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

#include "ajustePlano.h"

void acerosdAjustePlano (tajup ajup)
{
  ajup->n    = 0;
  ajup->sx   = 0.0;  ajup->sy   = 0.0;  ajup->sz   = 0.0;
  ajup->sxy  = 0.0;  ajup->sxz  = 0.0;  ajup->syz  = 0.0;
  ajup->sx2  = 0.0;  ajup->sy2  = 0.0;  ajup->sz2  = 0.0;
  ajup->az   = 0.0;  ajup->bz   = 0.0;  ajup->cz   = 0.0;
}

void nuevodAjustePlano (tajup * najup)
{
  tajup ajup;
  ajup = (tajup) malloc (sizeof(struct teajup));
  acerosdAjustePlano (ajup);
  *najup = ajup ;
}

void borradAjustePlano (tajup ajup)
{
  free (ajup);
}

void imprimedAjustePlano (tajup ajup)
{
  printf ("n:   %5d\n", ajup->n);
  printf ("sx:  %9.3f\n", ajup->sx);
  printf ("sy:  %9.3f\n", ajup->sy);
  printf ("sz:  %9.3f\n", ajup->sz);
  printf ("sxy: %9.3f\n", ajup->sxy);
  printf ("sxz: %9.3f\n", ajup->sxz);
  printf ("syz: %9.3f\n", ajup->syz);
  printf ("sx2: %9.3f\n", ajup->sx2);
  printf ("sy2: %9.3f\n", ajup->sy2);
  printf ("sz2: %9.3f\n", ajup->sz2);
  printf ("az:  %9.4f\n", ajup->az);
  printf ("bz:  %9.4f\n", ajup->bz);
  printf ("cz:  %9.4f\n", ajup->cz);
  printf ("\n");
}

void anyadedAAjustePlano (tajup ajup, int x, int y, int z)
{
  ajup->n    += 1 ;
  ajup->sx   += x ;      ajup->sy   += y ;      ajup->sz   += z ;
  ajup->sxy  += x * y ;  ajup->sxz  += x * z ;  ajup->syz  += y * z ;
  ajup->sx2  += x * x ;  ajup->sy2  += y * y ;  ajup->sz2  += z * z ;
}

void calculaAjustePlano (tajup ajup)
{
  double delz;
  delz =    ajup->sx2 * (ajup->sy  * ajup->sy - ajup->n * ajup->sy2)
	  + ajup->sx *   ajup->sx  * ajup->sy2
      - 2 * ajup->sx *   ajup->sy  * ajup->sxy
          + ajup->n *    ajup->sxy * ajup->sxy ;
  ajup->az = - ( ajup->sx *  (ajup->sy * ajup->syz - ajup->sy2 * ajup->sz)
                + ajup->sxy * (ajup->sy * ajup->sz -  ajup->n *   ajup->syz)
	       + (ajup->n *    ajup->sy2 - ajup->sy * ajup->sy) * ajup->sxz
	          ) / delz ;
  ajup->bz =   ( ajup->sx2 * (ajup->sy * ajup->sz -  ajup->n *   ajup->syz)
                + ajup->sx *   ajup->sx *  ajup->syz
               + (ajup->n *    ajup->sxy - ajup->sx * ajup->sy) * ajup->sxz
		- ajup->sx *   ajup->sxy * ajup->sz
		  ) / delz ;
  ajup->cz = - ( ajup->sx2 * (ajup->sy2 * ajup->sz - ajup->sy * ajup->syz)
	        + ajup->sx *   ajup->sxy * ajup->syz
               + (ajup->sxy *  ajup->sy -  ajup->sx * ajup->sy2) * ajup->sxz
	        - ajup->sxy *  ajup->sxy * ajup->sz
		  ) / delz ;
}

double distanciaAPlano (tajup ajup, int x, int y, int z)
{
  return (z - (ajup->az*x + ajup->bz*y + ajup->cz)) ;
}