/* expfit.c -- model functions for exponential + background */

struct data {
  size_t n;
  double * y;
};

int
expb_f (const gsl_vector * x, void *data, gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *y = ((struct data *)data)->y;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp (-lambda * t) + b;
      gsl_vector_set (f, i, Yi - y[i]);
    }

  return GSL_SUCCESS;
}

int
gaussian_f1peak (const gsl_vector * x, void *data, gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *y = ((struct data *)data)->y;
  // amplitude
  double x1 = gsl_vector_get(x, 0);
  // sigma
  double x2 = gsl_vector_get(x, 1);
  // peak
  double x3 = gsl_vector_get(x, 2);
  size_t i;

  for (i = 0; i < n; ++i)
    {
//      double ti = (7.0 - i) / 2.0;
      double ti = (x3 - i) / 2.0;
      double yi = y[i];
//      double term = ti - x3;
      double term = ti;
      double fi = x1 * exp(-x2*term*term/2.0) - yi;
      gsl_vector_set(f, i, fi);
    }

  return GSL_SUCCESS;
}

int
gaussian_f2peak (const gsl_vector * x, void *data, gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *y = ((struct data *)data)->y;
  // amplitude, sigma, peak
  double x11 = gsl_vector_get(x, 0);
  double x12 = gsl_vector_get(x, 1);
  double x13 = gsl_vector_get(x, 2);
  double x21 = gsl_vector_get(x, 3);
  double x22 = gsl_vector_get(x, 4);
  double x23 = gsl_vector_get(x, 5);
  size_t i;

  for (i = 0; i < n; ++i)
    {
      double ti = (x13 - i) / 2.0;
//	  double ti = i;
      double yi = y[i];
//      double term = ti - x3;
      double term = ti;
      double fi = x11 * exp(-x12*term*term/2.0);

      ti = (x23 - i) / 2.0;
      term = ti;
      fi += x21 * exp(-x22*term*term/2.0);

      fi -= yi;
      gsl_vector_set(f, i, fi);
    }

  return GSL_SUCCESS;
}

int
gaussian_f (const gsl_vector * x, void *data, gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *y = ((struct data *)data)->y;
  // amplitude, sigma, peak
  double x11 = gsl_vector_get(x, 0);
  double x12 = gsl_vector_get(x, 1);
  double x13 = gsl_vector_get(x, 2);
  double x21 = gsl_vector_get(x, 3);
  double x22 = gsl_vector_get(x, 4);
  double x23 = gsl_vector_get(x, 5);
  double x31 = gsl_vector_get(x, 6);
  double x32 = gsl_vector_get(x, 7);
  double x33 = gsl_vector_get(x, 8);
  size_t i;

  for (i = 0; i < n; ++i)
    {
      double ti = (x13 - i) / 2.0;
//	  double ti = i;
      double yi = y[i];
//      double term = ti - x3;
      double term = ti;
      double fi = x11 * exp(-x12*term*term/2.0);

      ti = (x23 - i) / 2.0;
      term = ti;
      fi += x21 * exp(-x22*term*term/2.0);

      ti = (x33 - i) / 2.0;
      term = ti;
      fi += x31 * exp(-x32*term*term/2.0);

      fi -= yi;
      gsl_vector_set(f, i, fi);
    }

  return GSL_SUCCESS;
}

int
expb_df (const gsl_vector * x, void *data, 
         gsl_matrix * J)
{
  size_t n = ((struct data *)data)->n;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double e = exp(-lambda * t);
      gsl_matrix_set (J, i, 0, e); 
      gsl_matrix_set (J, i, 1, -t * A * e);
      gsl_matrix_set (J, i, 2, 1.0);
    }
  return GSL_SUCCESS;
}