htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/expfit.c

   1 /* expfit.c -- model functions for exponential + background */
   2
   3 struct data {
   4   size_t n;
   5   double * y;
   6   double * sigma;
   7 };
   8
   9 int
  10 expb_f (const gsl_vector * x, void *data,
  11         gsl_vector * f)
  12 {
  13   size_t n = ((struct data *)data)->n;
  14   double *y = ((struct data *)data)->y;
  15   double *sigma = ((struct data *) data)->sigma;
  16
  17   double A = gsl_vector_get (x, 0);
  18   double lambda = gsl_vector_get (x, 1);
  19   double b = gsl_vector_get (x, 2);
  20
  21   size_t i;
  22
  23   for (i = 0; i < n; i++)
  24     {
  25       /* Model Yi = A * exp(-lambda * i) + b */
  26       double t = i;
  27       double Yi = A * exp (-lambda * t) + b;
  28       gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
  29     }
  30
  31   return GSL_SUCCESS;
  32 }
  33
  34 int
  35 expb_df (const gsl_vector * x, void *data,
  36          gsl_matrix * J)
  37 {
  38   size_t n = ((struct data *)data)->n;
  39   double *sigma = ((struct data *) data)->sigma;
  40
  41   double A = gsl_vector_get (x, 0);
  42   double lambda = gsl_vector_get (x, 1);
  43
  44   size_t i;
  45
  46   for (i = 0; i < n; i++)
  47     {
  48       /* Jacobian matrix J(i,j) = dfi / dxj, */
  49       /* where fi = (Yi - yi)/sigma[i],      */
  50       /*       Yi = A * exp(-lambda * i) + b  */
  51       /* and the xj are the parameters (A,lambda,b) */
  52       double t = i;
  53       double s = sigma[i];
  54       double e = exp(-lambda * t);
  55       gsl_matrix_set (J, i, 0, e/s);
  56       gsl_matrix_set (J, i, 1, -t * A * e/s);
  57       gsl_matrix_set (J, i, 2, 1/s);
  58     }
  59   return GSL_SUCCESS;
  60 }
  61
  62 int
  63 expb_fdf (const gsl_vector * x, void *data,
  64           gsl_vector * f, gsl_matrix * J)
  65 {
  66   expb_f (x, data, f);
  67   expb_df (x, data, J);
  68
  69   return GSL_SUCCESS;
  70 }