# include "psphot.h"
# include <stdlib.h>
double drand48(void);

int main (int argc, char **argv) {

  int i, j, N;
  psPolynomial2D *poly;

  srand48(0);

  // *** test 2D fit without errors  
  { 
    psVector *x  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *y  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *z  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *dz = NULL;

    N = 0;
    for (i = 0; i < 50; i++) {
      for (j = 0; j < 20; j++, N++) {
	x->data.F64[N] = i;
	y->data.F64[N] = j;
	z->data.F64[N] = 3 + 5*i - 2*j;
      }
    }
  
    // use the code in psUtils.c
    poly = Polynomial2DAlloc (1, 1, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial2DOrd_EAM (poly, x, y, z, dz);

    for (i = 0; i < poly->nX + 1; i++) {
      for (j = 0; j < poly->nY + 1; j++) {
	fprintf (stderr, "x^%d y^%d : %6.3f +/- %6.3f\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]);
      }
    }
  }


  // *** test 2D fit with errors  
  { 
    psVector *x  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *y  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *z  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *dz = psVectorAlloc (1000, PS_TYPE_F64);

    N = 0;
    for (i = 0; i < 50; i++) {
      for (j = 0; j < 20; j++, N++) {
	x->data.F64[N] = i;
	y->data.F64[N] = j;
	dz->data.F64[N] = 1*(drand48() - 0.5);
	z->data.F64[N] = 3 + 5*i - 2*j + dz->data.F64[N];
      }
    }
  
    // use the code in psUtils.c
    poly = Polynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial2DOrd_EAM (poly, x, y, z, dz);

    for (i = 0; i < poly->nX; i++) {
      for (j = 0; j < poly->nY; j++) {
	fprintf (stderr, "x^%d y^%d : %f +/- %f\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]);
      }
    }
  }
}
