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

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

  int i;
  psPolynomial1D *poly;

  // *** test 1D fit without errors  
  { 
    psVector *x  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *y  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *dy = NULL;
    psVector *dy2 = NULL;

    srand48(0);

    for (i = 0; i < 1000; i++) {
      x->data.F64[i] = i;
      y->data.F64[i] = 3 + 5*i - 2*i*i +0.4*i*i*i;
    }
  
    // use the code in psUtils.c
    poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial1DOrd_EAM (poly, x, y, dy);

    for (i = 0; i < poly->n; i++) {
      fprintf (stderr, "x^%d : %f +/- %f\n", i, poly->coeff[i], poly->coeffErr[i]);
    }
  
    // use the code in psMinimize.c
    poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial1DOrd (poly, x, y, dy2);

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

  // *** test 1D fit with errors  
  {
    psVector *x  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *y  = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *dy = psVectorAlloc (1000, PS_TYPE_F64);
    psVector *dy2 = psVectorAlloc (1000, PS_TYPE_F64);

    srand48(0);

    for (i = 0; i < 1000; i++) {
      x->data.F64[i] = i;
      dy->data.F64[i] = 10*(drand48() - 0.5);
      dy2->data.F64[i] = PS_SQR(dy->data.F64[i]);
      y->data.F64[i] = 3 + 5*i - 2*i*i +0.4*i*i*i + dy->data.F64[i];
    }
  
    // use the code in psUtils.c
    poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial1DOrd_EAM (poly, x, y, dy);

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

    // use the code in psMinimize.c
    poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD);
    poly = VectorFitPolynomial1DOrd (poly, x, y, dy2);

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