#define enso_N       168

#define enso_P       9

static double enso_x0a[enso_P] = { 11.0, 3.0, 0.5, 40.0, -0.7, -1.3, 25.0, -0.3, 1.4 };

static double enso_x0b[enso_P] = { 10.0, 3.0, 0.5, 44.0, -1.5, 0.5, 26.0, -0.1, 1.5 };

static double enso_epsrel = 1.0e-3;

static double enso_sigma[enso_P] = {

 1.7488832467E-01, 2.4310052139E-01, 2.4354686618E-01, 9.4408025976E-01, 2.8078369611E-01,

 4.8073701119E-01, 4.1612939130E-01, 5.1460022911E-01, 2.5434468893E-01

};

static double enso_F[enso_N] = {

    12.90000, 11.30000, 10.60000, 11.20000, 10.90000, 

    7.500000, 7.700000, 11.70000, 12.90000, 14.30000, 

    10.90000, 13.70000, 17.10000, 14.00000, 15.30000, 

    8.500000, 5.700000, 5.500000, 7.600000, 8.600000, 

    7.300000, 7.600000, 12.70000, 11.00000, 12.70000, 

    12.90000, 13.00000, 10.90000, 10.400000, 10.200000, 

    8.000000, 10.90000, 13.60000, 10.500000, 9.200000, 

    12.40000, 12.70000, 13.30000, 10.100000, 7.800000, 

    4.800000, 3.000000, 2.500000, 6.300000, 9.700000, 

    11.60000, 8.600000, 12.40000, 10.500000, 13.30000, 

   10.400000, 8.100000, 3.700000, 10.70000, 5.100000, 

   10.400000, 10.90000, 11.70000, 11.40000, 13.70000, 

    14.10000, 14.00000, 12.50000, 6.300000, 9.600000, 

    11.70000, 5.000000, 10.80000, 12.70000, 10.80000, 

    11.80000, 12.60000, 15.70000, 12.60000, 14.80000, 

    7.800000, 7.100000, 11.20000, 8.100000, 6.400000, 

    5.200000, 12.00000, 10.200000, 12.70000, 10.200000, 

    14.70000, 12.20000, 7.100000, 5.700000, 6.700000, 

    3.900000, 8.500000, 8.300000, 10.80000, 16.70000, 

    12.60000, 12.50000, 12.50000, 9.800000, 7.200000, 

    4.100000, 10.60000, 10.100000, 10.100000, 11.90000, 

    13.60000, 16.30000, 17.60000, 15.50000, 16.00000, 

    15.20000, 11.20000, 14.30000, 14.50000, 8.500000, 

    12.00000, 12.70000, 11.30000, 14.50000, 15.10000, 

   10.400000, 11.50000, 13.40000, 7.500000, 0.6000000, 

   0.3000000, 5.500000, 5.000000, 4.600000, 8.200000, 

    9.900000, 9.200000, 12.50000, 10.90000, 9.900000, 

    8.900000, 7.600000, 9.500000, 8.400000, 10.70000, 

    13.60000, 13.70000, 13.70000, 16.50000, 16.80000, 

    17.10000, 15.40000, 9.500000, 6.100000, 10.100000, 

    9.300000, 5.300000, 11.20000, 16.60000, 15.60000, 

    12.00000, 11.50000, 8.600000, 13.80000, 8.700000, 

    8.600000, 8.600000, 8.700000, 12.80000, 13.20000, 

    14.00000, 13.40000, 14.80000

};

static void

enso_checksol(const double x[], const double sumsq,

              const double epsrel, const char *sname,

              const char *pname)

{

  size_t i;

  const double sumsq_exact = 7.8853978668E+02;

  const double enso_x[enso_P] = {

    1.0510749193E+01, 3.0762128085E+00, 5.3280138227E-01,

    4.4311088700E+01, -1.6231428586E+00, 5.2554493756E-01,

    2.6887614440E+01, 2.1232288488E-01, 1.4966870418E+00 };

  gsl_test_rel(sumsq, sumsq_exact, epsrel, "%s/%s sumsq",

               sname, pname);

  for (i = 0; i < enso_P; ++i)

    {

      gsl_test_rel(x[i], enso_x[i], epsrel, "%s/%s i=%zu",

                   sname, pname, i);

    }

}

static int

enso_f (const gsl_vector * x, void *params, gsl_vector * f)

{

  double b[enso_P];

  size_t i;

  for (i = 0; i < enso_P; i++)

    {

      b[i] = gsl_vector_get(x, i);

    }

  for (i = 0; i < enso_N; i++)

    {

      double t = (i + 1.0);

      double y;

      y = b[0];

      y += b[1] * cos(2*M_PI*t/12);

      y += b[2] * sin(2*M_PI*t/12);

      y += b[4] * cos(2*M_PI*t/b[3]);

      y += b[5] * sin(2*M_PI*t/b[3]);

      y += b[7] * cos(2*M_PI*t/b[6]);

      y += b[8] * sin(2*M_PI*t/b[6]);

      gsl_vector_set (f, i, enso_F[i] - y);

    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;

}

static int

enso_df (const gsl_vector * x, void *params, gsl_matrix * df)

{

  double b[enso_P];

  size_t i;

  for (i = 0; i < enso_P; i++)

    {

      b[i] = gsl_vector_get(x, i);

    }

  for (i = 0; i < enso_N; i++)

    {

      double t = (i + 1.0);

      gsl_matrix_set (df, i, 0, -1.0);

      gsl_matrix_set (df, i, 1, -cos(2*M_PI*t/12));

      gsl_matrix_set (df, i, 2, -sin(2*M_PI*t/12));

      gsl_matrix_set (df, i, 3, 

                      -b[4]*(2*M_PI*t/(b[3]*b[3]))*sin(2*M_PI*t/b[3])

                      +b[5]*(2*M_PI*t/(b[3]*b[3]))*cos(2*M_PI*t/b[3]));

      gsl_matrix_set (df, i, 4, -cos(2*M_PI*t/b[3]));

      gsl_matrix_set (df, i, 5, -sin(2*M_PI*t/b[3]));

      gsl_matrix_set (df, i, 6, 

                     -b[7] * (2*M_PI*t/(b[6]*b[6])) * sin(2*M_PI*t/b[6])

                     +b[8] * (2*M_PI*t/(b[6]*b[6])) * cos(2*M_PI*t/b[6]));

      gsl_matrix_set (df, i, 7, -cos(2*M_PI*t/b[6]));

      gsl_matrix_set (df, i, 8, -sin(2*M_PI*t/b[6]));

    }

  (void)params; /* avoid unused parameter warning */

  return GSL_SUCCESS;

}

static gsl_multifit_nlinear_fdf enso_func =

{

  enso_f,

  enso_df,

  NULL, /* analytic expression too complex */

  enso_N,

  enso_P,

  NULL,

  0,

  0,

  0

};

static test_fdf_problem ensoa_problem =

{

  "nist-ENSOa",

  enso_x0a,

  NULL,

  enso_sigma,

  &enso_epsrel,

  &enso_checksol,

  &enso_func

};

static test_fdf_problem ensob_problem =

{

  "nist-ENSOb",

  enso_x0b,

  NULL,

  enso_sigma,

  &enso_epsrel,

  &enso_checksol,

  &enso_func

};