htswanalysis/MACS/lib/gsl/gsl-1.11/interpolation/test.c

   1 /* interpolation/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20
  21 /* Author:  G. Jungman
  22  */
  23 #include <config.h>
  24 #include <stddef.h>
  25 #include <stdlib.h>
  26 #include <stdio.h>
  27 #include <math.h>
  28 #include <gsl/gsl_test.h>
  29 #include <gsl/gsl_errno.h>
  30 #include <gsl/gsl_interp.h>
  31 #include <gsl/gsl_ieee_utils.h>
  32
  33 int
  34 test_bsearch(void)
  35 {
  36   double x_array[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  37   size_t index_result;
  38   int status = 0;
  39   int s;
  40
  41   /* check an interior point */
  42   index_result = gsl_interp_bsearch(x_array, 1.5, 0, 4);
  43   s = (index_result != 1);
  44   status += s;
  45   gsl_test (s, "simple bsearch");
  46
  47   /* check that we get the last interval if x == last value */
  48   index_result = gsl_interp_bsearch(x_array, 4.0, 0, 4);
  49   s = (index_result != 3);
  50   status += s;
  51   gsl_test (s, "upper endpoint bsearch");
  52
  53   /* check that we get the first interval if x == first value */
  54   index_result = gsl_interp_bsearch(x_array, 0.0, 0, 4);
  55   s = (index_result != 0);
  56   status += s;
  57   gsl_test (s, "lower endpoint bsearch");
  58
  59   /* check that we get correct interior boundary behaviour */
  60   index_result = gsl_interp_bsearch(x_array, 2.0, 0, 4);
  61   s = (index_result != 2);
  62   status += s;
  63   gsl_test (s, "degenerate bsearch");
  64
  65   /* check out of bounds above */
  66   index_result = gsl_interp_bsearch(x_array, 10.0, 0, 4);
  67   s = (index_result != 3);
  68   status += s;
  69   gsl_test (s, "out of bounds bsearch +");
  70
  71   /* check out of bounds below */
  72   index_result = gsl_interp_bsearch(x_array, -10.0, 0, 4);
  73   s = (index_result != 0);
  74   status += s;
  75   gsl_test (s, "out of bounds bsearch -");
  76
  77   return status;
  78 }
  79
  80
  81
  82 typedef double TEST_FUNC (double);
  83 typedef struct _xy_table xy_table;
  84
  85 struct _xy_table
  86   {
  87     double * x;
  88     double * y;
  89     size_t n;
  90   };
  91
  92 xy_table make_xy_table (double x[], double y[], size_t n);
  93
  94 xy_table
  95 make_xy_table (double x[], double y[], size_t n)
  96 {
  97   xy_table t;
  98   t.x = x;
  99   t.y = y;
 100   t.n = n;
 101   return t;
 102 }
 103
 104 static int
 105 test_interp (
 106   const xy_table * data_table,
 107   const gsl_interp_type * T,
 108   xy_table * test_table,
 109   xy_table * test_d_table,
 110   xy_table * test_i_table
 111   )
 112 {
 113   int status = 0, s1, s2, s3;
 114   size_t i;
 115
 116   gsl_interp_accel *a = gsl_interp_accel_alloc ();
 117   gsl_interp *interp = gsl_interp_alloc (T, data_table->n);
 118
 119   gsl_interp_init (interp, data_table->x, data_table->y, data_table->n);
 120
 121   for (i = 0; i < test_table->n; i++)
 122     {
 123       double x = test_table->x[i];
 124       double y;
 125       double deriv;
 126       double integ;
 127       double diff_y, diff_deriv, diff_integ;
 128       s1 = gsl_interp_eval_e (interp, data_table->x, data_table->y, x, a, &y);
 129       s2 = gsl_interp_eval_deriv_e (interp, data_table->x, data_table->y, x, a, &deriv);
 130       s3 = gsl_interp_eval_integ_e (interp, data_table->x, data_table->y, test_table->x[0], x, a, &integ);
 131
 132       gsl_test (s1, "gsl_interp_eval_e %s", gsl_interp_name(interp));
 133       gsl_test (s2, "gsl_interp_eval_deriv_e %s", gsl_interp_name(interp));
 134       gsl_test (s3, "gsl_interp_eval_integ_e %s", gsl_interp_name(interp));
 135
 136       gsl_test_abs (y, test_table->y[i], 1e-10, "%s %d", gsl_interp_name(interp), i);
 137       gsl_test_abs (deriv, test_d_table->y[i], 1e-10, "%s deriv %d", gsl_interp_name(interp), i);
 138       gsl_test_abs (integ, test_i_table->y[i], 1e-10, "%s integ %d", gsl_interp_name(interp), i);
 139
 140       diff_y = y - test_table->y[i];
 141       diff_deriv = deriv - test_d_table->y[i];
 142       diff_integ = integ - test_i_table->y[i];
 143       if (fabs (diff_y) > 1.e-10 || fabs(diff_deriv) > 1.0e-10 || fabs(diff_integ) > 1.0e-10) {
 144         status++;
 145       }
 146     }
 147
 148   gsl_interp_accel_free (a);
 149   gsl_interp_free (interp);
 150
 151   return status;
 152 }
 153
 154 static int
 155 test_linear (void)
 156 {
 157   int s;
 158
 159   double data_x[4] = { 0.0, 1.0, 2.0, 3.0 };
 160   double data_y[4] = { 0.0, 1.0, 2.0, 3.0 };
 161   double test_x[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
 162   double test_y[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
 163   double test_dy[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
 164   double test_iy[6] = { 0.0, 0.125, 0.5, 9.0/8.0, 25.0/8.0, 9.0/2.0 };
 165
 166   xy_table data_table = make_xy_table(data_x, data_y, 4);
 167   xy_table test_table = make_xy_table(test_x, test_y, 6);
 168   xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
 169   xy_table test_i_table = make_xy_table(test_x, test_iy, 6);
 170
 171   s = test_interp (&data_table, gsl_interp_linear, &test_table, &test_d_table, &test_i_table);
 172   gsl_test (s, "linear interpolation");
 173   return s;
 174 }
 175
 176 static int
 177 test_polynomial (void)
 178 {
 179   int s;
 180
 181   double data_x[4] = { 0.0, 1.0, 2.0, 3.0 };
 182   double data_y[4] = { 0.0, 1.0, 2.0, 3.0 };
 183   double test_x[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
 184   double test_y[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
 185   double test_dy[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
 186   double test_iy[6] = { 0.0, 0.125, 0.5, 9.0/8.0, 25.0/8.0, 9.0/2.0 };
 187
 188   xy_table data_table = make_xy_table(data_x, data_y, 4);
 189   xy_table test_table = make_xy_table(test_x, test_y, 6);
 190   xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
 191   xy_table test_i_table = make_xy_table(test_x, test_iy, 6);
 192
 193   s = test_interp (&data_table, gsl_interp_polynomial, &test_table, &test_d_table, &test_i_table);
 194   gsl_test (s, "polynomial interpolation");
 195   return s;
 196 }
 197
 198
 199 static int
 200 test_cspline (void)
 201 {
 202   int s;
 203
 204   double data_x[3] = { 0.0, 1.0, 2.0 };
 205   double data_y[3] = { 0.0, 1.0, 2.0 };
 206   double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
 207   double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
 208   double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
 209   double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };
 210
 211   xy_table data_table = make_xy_table(data_x, data_y, 3);
 212   xy_table test_table = make_xy_table(test_x, test_y, 4);
 213   xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
 214   xy_table test_i_table = make_xy_table(test_x, test_iy, 4);
 215
 216   s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
 217   gsl_test (s, "cspline interpolation");
 218   return s;
 219 }
 220
 221 static int
 222 test_cspline2 (void)
 223 {
 224   /* Test taken from Young & Gregory, A Survey of Numerical
 225      Mathematics, Vol 1 Chapter 6.8 */
 226
 227   int s;
 228
 229   double data_x[6] = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };
 230
 231   double data_y[6] = { 1.0,
 232                        0.961538461538461, 0.862068965517241,
 233                        0.735294117647059, 0.609756097560976,
 234                        0.500000000000000 } ;
 235
 236   double test_x[50] = {
 237     0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18,
 238     0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38,
 239     0.40, 0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58,
 240     0.60, 0.62, 0.64, 0.66, 0.68, 0.70, 0.72, 0.74, 0.76, 0.78,
 241     0.80, 0.82, 0.84, 0.86, 0.88, 0.90, 0.92, 0.94, 0.96, 0.98 };
 242
 243   double test_y[50] = {
 244     1.000000000000000, 0.997583282975581, 0.995079933416512,
 245     0.992403318788142, 0.989466806555819, 0.986183764184894,
 246     0.982467559140716, 0.978231558888635, 0.973389130893999,
 247     0.967853642622158, 0.961538461538461, 0.954382579685350,
 248     0.946427487413627, 0.937740299651188, 0.928388131325928,
 249     0.918438097365742, 0.907957312698524, 0.897012892252170,
 250     0.885671950954575, 0.874001603733634, 0.862068965517241,
 251     0.849933363488199, 0.837622973848936, 0.825158185056786,
 252     0.812559385569085, 0.799846963843167, 0.787041308336369,
 253     0.774162807506023, 0.761231849809467, 0.748268823704033,
 254     0.735294117647059, 0.722328486073082, 0.709394147325463,
 255     0.696513685724764, 0.683709685591549, 0.671004731246381,
 256     0.658421407009825, 0.645982297202442, 0.633709986144797,
 257     0.621627058157454, 0.609756097560976, 0.598112015427308,
 258     0.586679029833925, 0.575433685609685, 0.564352527583445,
 259     0.553412100584061, 0.542588949440392, 0.531859618981294,
 260     0.521200654035625, 0.510588599432241};
 261
 262   double test_dy[50] = {
 263     -0.120113913432180, -0.122279726798445, -0.128777166897241,
 264     -0.139606233728568, -0.154766927292426, -0.174259247588814,
 265     -0.198083194617734, -0.226238768379184, -0.258725968873165,
 266     -0.295544796099676, -0.336695250058719, -0.378333644186652,
 267     -0.416616291919835, -0.451543193258270, -0.483114348201955,
 268     -0.511329756750890, -0.536189418905076, -0.557693334664512,
 269     -0.575841504029200, -0.590633926999137, -0.602070603574326,
 270     -0.611319695518765, -0.619549364596455, -0.626759610807396,
 271     -0.632950434151589, -0.638121834629033, -0.642273812239728,
 272     -0.645406366983674, -0.647519498860871, -0.648613207871319,
 273     -0.648687494015019, -0.647687460711257, -0.645558211379322,
 274     -0.642299746019212, -0.637912064630930, -0.632395167214473,
 275     -0.625749053769843, -0.617973724297039, -0.609069178796061,
 276     -0.599035417266910, -0.587872439709585, -0.576731233416743,
 277     -0.566762785681043, -0.557967096502484, -0.550344165881066,
 278     -0.543893993816790, -0.538616580309654, -0.534511925359660,
 279     -0.531580028966807, -0.529820891131095};
 280
 281   double test_iy[50] = {
 282     0.000000000000000, 0.019975905023535, 0.039902753768792,
 283     0.059777947259733, 0.079597153869625, 0.099354309321042,
 284     0.119041616685866, 0.138649546385285, 0.158166836189794,
 285     0.177580491219196, 0.196875783942601, 0.216036382301310,
 286     0.235045759060558, 0.253888601161251, 0.272550937842853,
 287     0.291020140643388, 0.309284923399436, 0.327335342246135,
 288     0.345162795617181, 0.362760024244829, 0.380121111159890,
 289     0.397241442753010, 0.414117280448683, 0.430745332379281,
 290     0.447122714446318, 0.463246950320456, 0.479115971441505,
 291     0.494728117018421, 0.510082134029305, 0.525177177221407,
 292     0.540012809111123, 0.554589001813881, 0.568906157172889,
 293     0.582965126887879, 0.596767214344995, 0.610314174616794,
 294     0.623608214462242, 0.636651992326715, 0.649448618342004,
 295     0.662001654326309, 0.674315113784241, 0.686393423540581,
 296     0.698241001711602, 0.709861835676399, 0.721259443710643,
 297     0.732436874986582, 0.743396709573044, 0.754141058435429,
 298     0.764671563435718, 0.774989397332469 };
 299
 300   xy_table data_table = make_xy_table(data_x, data_y, 6);
 301   xy_table test_table = make_xy_table(test_x, test_y, 50);
 302   xy_table test_d_table = make_xy_table(test_x, test_dy, 50);
 303   xy_table test_i_table = make_xy_table(test_x, test_iy, 50);
 304
 305   s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
 306   gsl_test (s, "cspline 1/(1+x^2) interpolation");
 307   return s;
 308 }
 309
 310 static int
 311 test_cspline3 (void)
 312 {
 313   /* This data has been chosen to be random (uneven spacing in x and y)
 314      and the exact cubic spine solution computed using Octave */
 315
 316   int s;
 317
 318   double data_x[7] = {
 319     -1.2139767065644265,
 320     -0.792590494453907,
 321     -0.250954683125019,
 322     0.665867809951305,
 323     0.735655088722706,
 324     0.827622053027153,
 325     1.426592227816582
 326   };
 327
 328   double data_y[7] = {
 329     -0.00453877449035645,
 330     0.49763182550668716,
 331     0.17805472016334534,
 332     0.40514493733644485,
 333     -0.21595209836959839,
 334     0.47405586764216423,
 335     0.46561462432146072
 336   } ;
 337
 338   double test_x[19] = {
 339     -1.2139767065644265,
 340     -1.0735146358609200,
 341     -0.9330525651574135,
 342     -0.7925904944539071,
 343     -0.6120452240109444,
 344     -0.4314999535679818,
 345     -0.2509546831250191,
 346     0.0546528145670890,
 347     0.3602603122591972,
 348     0.6658678099513053,
 349     0.6891302362084388,
 350     0.7123926624655723,
 351     0.7356550887227058,
 352     0.7663107434908548,
 353     0.7969663982590039,
 354     0.8276220530271530,
 355     1.0272787779569625,
 356     1.2269355028867721,
 357     1.4265922278165817,
 358   };
 359
 360   double test_y[19] = {
 361     -0.00453877449035645,
 362     0.25816917628390590,
 363     0.44938881397673230,
 364     0.49763182550668716,
 365     0.31389980410075147,
 366     0.09948951681196887,
 367     0.17805472016334534,
 368     1.27633142487980233,
 369     2.04936553432792001,
 370     0.40514493733644485,
 371     0.13322324792901385,
 372     -0.09656315924697809,
 373     -0.21595209836959839,
 374     -0.13551147728045118,
 375     0.13466779030061801,
 376     0.47405586764216423,
 377     1.68064089899304370,
 378     1.43594739539458649,
 379     0.46561462432146072
 380   };
 381
 382   double test_dy[19] = {
 383     1.955137555965937,
 384     1.700662049790549,
 385     0.937235531264386,
 386     -0.335141999612553,
 387     -1.401385073563169,
 388     -0.674982149482761,
 389     1.844066772628670,
 390     4.202528085784793,
 391     -0.284432022227558,
 392     -11.616813551408383,
 393     -11.272731243226174,
 394     -7.994209291156876,
 395     -1.781247695200491,
 396     6.373970868827501,
 397     10.597456848997197,
 398     10.889210245308570,
 399     1.803124267866902,
 400     -3.648527318598099,
 401     -5.465744514086432,
 402   };
 403
 404   double test_iy[19] = {
 405     0.000000000000000,
 406     0.018231117234914,
 407     0.069178822023139,
 408     0.137781019634897,
 409     0.213936442847744,
 410     0.249280997744777,
 411     0.267492946016120,
 412     0.471372708120518,
 413     1.014473660088477,
 414     1.477731933018837,
 415     1.483978291717981,
 416     1.484256847945450,
 417     1.480341742628893,
 418     1.474315901028282,
 419     1.473972210647307,
 420     1.483279773396950,
 421     1.728562698140330,
 422     2.057796448999396,
 423     2.253662886537457,
 424   };
 425
 426   xy_table data_table = make_xy_table(data_x, data_y, 7);
 427   xy_table test_table = make_xy_table(test_x, test_y, 19);
 428   xy_table test_d_table = make_xy_table(test_x, test_dy, 19);
 429   xy_table test_i_table = make_xy_table(test_x, test_iy, 19);
 430
 431   s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
 432   gsl_test (s, "cspline arbitrary data interpolation");
 433   return s;
 434 }
 435
 436 static int
 437 test_csplinep (void)
 438 {
 439   /* This data has been chosen to be random (uneven spacing in x and y)
 440      and the exact cubic spine solution computed using Octave */
 441
 442   int s;
 443
 444   double data_x[11] = {
 445     0.000000000000000,
 446     0.130153674349869,
 447     0.164545962312740,
 448     0.227375461261537,
 449     0.256465324353657,
 450     0.372545206874658,
 451     0.520820016781720,
 452     0.647654717733075,
 453     0.753429306654340,
 454     0.900873984827658,
 455     1.000000000000000,
 456   };
 457
 458   double data_y[11] = {
 459     0.000000000000000,
 460     0.729629261832041,
 461     0.859286331568207,
 462     0.989913099419008,
 463     0.999175006262120,
 464     0.717928599519215,
 465     -0.130443237213363,
 466     -0.800267961158980,
 467     -0.999767873040527,
 468     -0.583333769240853,
 469     -0.000000000000000,
 470   } ;
 471
 472   double test_x[31] = {
 473     0.000000000000000,
 474     0.043384558116623,
 475     0.086769116233246,
 476     0.130153674349869,
 477     0.141617770337492,
 478     0.153081866325116,
 479     0.164545962312740,
 480     0.185489128629005,
 481     0.206432294945271,
 482     0.227375461261537,
 483     0.237072082292243,
 484     0.246768703322951,
 485     0.256465324353657,
 486     0.295158618527324,
 487     0.333851912700991,
 488     0.372545206874658,
 489     0.421970143510346,
 490     0.471395080146033,
 491     0.520820016781720,
 492     0.563098250432172,
 493     0.605376484082623,
 494     0.647654717733075,
 495     0.682912914040164,
 496     0.718171110347252,
 497     0.753429306654340,
 498     0.802577532712113,
 499     0.851725758769885,
 500     0.900873984827658,
 501     0.933915989885105,
 502     0.966957994942553,
 503     1.000000000000000
 504   };
 505
 506   double test_y[31] = {
 507     0.000000000000000,
 508     0.268657574670719,
 509     0.517940878523929,
 510     0.729629261832041,
 511     0.777012551497867,
 512     0.820298314554859,
 513     0.859286331568207,
 514     0.918833991960315,
 515     0.962624749226346,
 516     0.989913099419008,
 517     0.996756196601349,
 518     0.999858105635752,
 519     0.999175006262120,
 520     0.959248551766306,
 521     0.863713527741856,
 522     0.717928599519215,
 523     0.470065187871106,
 524     0.177694938589523,
 525     -0.130443237213363,
 526     -0.385093922365765,
 527     -0.613840011545983,
 528     -0.800267961158980,
 529     -0.912498361131651,
 530     -0.980219217412290,
 531     -0.999767873040528,
 532     -0.943635958253643,
 533     -0.800314354800596,
 534     -0.583333769240853,
 535     -0.403689914131666,
 536     -0.206151346799382,
 537     -0.000000000000000
 538   };
 539
 540   double test_dy[31] = {
 541     6.275761975917336,
 542     6.039181349093287,
 543     5.382620652895025,
 544     4.306079887322550,
 545     3.957389777282752,
 546     3.591234754612763,
 547     3.207614819312586,
 548     2.473048186927024,
 549     1.702885029353488,
 550     0.897125346591982,
 551     0.513561009477969,
 552     0.125477545550710,
 553     -0.267125045189792,
 554     -1.773533414873785,
 555     -3.141450982859891,
 556     -4.370877749148106,
 557     -5.562104227060234,
 558     -6.171864888961167,
 559     -6.200159734850907,
 560     -5.781556055213110,
 561     -4.974725570010514,
 562     -3.779668279243120,
 563     -2.569220222282655,
 564     -1.254891157670244,
 565     0.163318914594122,
 566     2.074985916277011,
 567     3.711348850147548,
 568     5.072407716205733,
 569     5.754438923510391,
 570     6.155557010080926,
 571     6.275761975917336
 572   };
 573
 574   double test_iy[31] = {
 575     0.000000000000000,
 576     0.005864903144198,
 577     0.023030998930796,
 578     0.050262495763030,
 579     0.058902457844100,
 580     0.068062330564832,
 581     0.077693991819461,
 582     0.096340576055474,
 583     0.116070578226521,
 584     0.136546192283223,
 585     0.146181187290769,
 586     0.155864434185569,
 587     0.165559443629720,
 588     0.203636318965633,
 589     0.239075190181586,
 590     0.269828050745236,
 591     0.299428805999600,
 592     0.315560685785616,
 593     0.316734151903742,
 594     0.305773798930857,
 595     0.284537037103156,
 596     0.254466034797342,
 597     0.224146112780097,
 598     0.190643050847000,
 599     0.155590744566140,
 600     0.107448508851745,
 601     0.064263083957312,
 602     0.029987183298960,
 603     0.013618510529526,
 604     0.003506827320794,
 605     0.000090064010007,
 606   };
 607
 608   xy_table data_table = make_xy_table(data_x, data_y, 11);
 609   xy_table test_table = make_xy_table(test_x, test_y, 31);
 610   xy_table test_d_table = make_xy_table(test_x, test_dy, 31);
 611   xy_table test_i_table = make_xy_table(test_x, test_iy, 31);
 612
 613   s = test_interp (&data_table, gsl_interp_cspline_periodic, &test_table, &test_d_table, &test_i_table);
 614   gsl_test (s, "cspline periodic sine interpolation");
 615   return s;
 616 }
 617
 618 static int
 619 test_csplinep2 (void)
 620 {
 621   /* This data tests the periodic case n=3 */
 622
 623   int s;
 624
 625   double data_x[3] = {
 626     0.123,
 627     0.423,
 628     1.123
 629   };
 630
 631   double data_y[3] = {
 632     0.456000000000000,
 633     1.407056516295154,
 634     0.456000000000000
 635   } ;
 636
 637   double test_x[61] = {
 638     0.123000000000000,
 639     0.133000000000000,
 640     0.143000000000000,
 641     0.153000000000000,
 642     0.163000000000000,
 643     0.173000000000000,
 644     0.183000000000000,
 645     0.193000000000000,
 646     0.203000000000000,
 647     0.213000000000000,
 648     0.223000000000000,
 649     0.233000000000000,
 650     0.243000000000000,
 651     0.253000000000000,
 652     0.263000000000000,
 653     0.273000000000000,
 654     0.283000000000000,
 655     0.293000000000000,
 656     0.303000000000000,
 657     0.313000000000000,
 658     0.323000000000000,
 659     0.333000000000000,
 660     0.343000000000000,
 661     0.353000000000000,
 662     0.363000000000000,
 663     0.373000000000000,
 664     0.383000000000000,
 665     0.393000000000000,
 666     0.403000000000000,
 667     0.413000000000000,
 668     0.423000000000000,
 669     0.446333333333333,
 670     0.469666666666667,
 671     0.493000000000000,
 672     0.516333333333333,
 673     0.539666666666667,
 674     0.563000000000000,
 675     0.586333333333333,
 676     0.609666666666667,
 677     0.633000000000000,
 678     0.656333333333333,
 679     0.679666666666667,
 680     0.703000000000000,
 681     0.726333333333333,
 682     0.749666666666667,
 683     0.773000000000000,
 684     0.796333333333333,
 685     0.819666666666667,
 686     0.843000000000000,
 687     0.866333333333333,
 688     0.889666666666667,
 689     0.913000000000000,
 690     0.936333333333333,
 691     0.959666666666667,
 692     0.983000000000000,
 693     1.006333333333333,
 694     1.029666666666667,
 695     1.053000000000000,
 696     1.076333333333333,
 697     1.099666666666667,
 698     1.123000000000000
 699   };
 700
 701   double test_y[61] = {
 702     0.456000000000000,
 703     0.475443822110923,
 704     0.497423794931967,
 705     0.521758764840979,
 706     0.548267578215809,
 707     0.576769081434305,
 708     0.607082120874316,
 709     0.639025542913690,
 710     0.672418193930275,
 711     0.707078920301921,
 712     0.742826568406475,
 713     0.779479984621787,
 714     0.816858015325704,
 715     0.854779506896076,
 716     0.893063305710751,
 717     0.931528258147577,
 718     0.969993210584403,
 719     1.008277009399078,
 720     1.046198500969450,
 721     1.083576531673367,
 722     1.120229947888679,
 723     1.155977595993233,
 724     1.190638322364879,
 725     1.224030973381464,
 726     1.255974395420838,
 727     1.286287434860848,
 728     1.314788938079344,
 729     1.341297751454174,
 730     1.365632721363187,
 731     1.387612694184230,
 732     1.407056516295154,
 733     1.442092968697928,
 734     1.463321489456714,
 735     1.471728359403224,
 736     1.468299859369172,
 737     1.454022270186272,
 738     1.429881872686237,
 739     1.396864947700781,
 740     1.355957776061616,
 741     1.308146638600458,
 742     1.254417816149018,
 743     1.195757589539010,
 744     1.133152239602149,
 745     1.067588047170148,
 746     1.000051293074719,
 747     0.931528258147577,
 748     0.863005223220435,
 749     0.795468469125006,
 750     0.729904276693004,
 751     0.667298926756143,
 752     0.608638700146136,
 753     0.554909877694696,
 754     0.507098740233537,
 755     0.466191568594372,
 756     0.433174643608916,
 757     0.409034246108881,
 758     0.394756656925981,
 759     0.391328156891929,
 760     0.399735026838439,
 761     0.420963547597225,
 762     0.456000000000000
 763   };
 764
 765   double test_dy[61] = {
 766     1.8115362215145774,
 767     2.0742089736341924,
 768     2.3187663635386602,
 769     2.5452083912279826,
 770     2.7535350567021584,
 771     2.9437463599611897,
 772     3.1158423010050744,
 773     3.2698228798338147,
 774     3.4056880964474079,
 775     3.5234379508458549,
 776     3.6230724430291570,
 777     3.7045915729973125,
 778     3.7679953407503231,
 779     3.8132837462881874,
 780     3.8404567896109061,
 781     3.8495144707184790,
 782     3.8404567896109061,
 783     3.8132837462881874,
 784     3.7679953407503231,
 785     3.7045915729973125,
 786     3.6230724430291565,
 787     3.5234379508458549,
 788     3.4056880964474074,
 789     3.2698228798338147,
 790     3.1158423010050749,
 791     2.9437463599611897,
 792     2.7535350567021584,
 793     2.5452083912279830,
 794     2.3187663635386597,
 795     2.0742089736341924,
 796     1.8115362215145772,
 797     1.1986331332354823,
 798     0.6279992234583869,
 799     0.0996344921833026,
 800     -0.3864610605897765,
 801     -0.8302874348608467,
 802     -1.2318446306299125,
 803     -1.5911326478969707,
 804     -1.9081514866620208,
 805     -2.1829011469250670,
 806     -2.4153816286861041,
 807     -2.6055929319451341,
 808     -2.7535350567021584,
 809     -2.8592080029571765,
 810     -2.9226117707101862,
 811     -2.9437463599611893,
 812     -2.9226117707101862,
 813     -2.8592080029571760,
 814     -2.7535350567021593,
 815     -2.6055929319451354,
 816     -2.4153816286861045,
 817     -2.1829011469250656,
 818     -1.9081514866620242,
 819     -1.5911326478969716,
 820     -1.2318446306299160,
 821     -0.8302874348608498,
 822     -0.3864610605897769,
 823     0.0996344921832995,
 824     0.6279992234583824,
 825     1.1986331332354769,
 826     1.8115362215145772
 827   };
 828
 829   double test_iy[61] = {
 830     0.000000000000000,
 831     0.004655030170954,
 832     0.009517330277919,
 833     0.014611356059886,
 834     0.019959751719625,
 835     0.025583349923681,
 836     0.031501171802382,
 837     0.037730426949832,
 838     0.044286513423914,
 839     0.051183017746288,
 840     0.058431714902395,
 841     0.066042568341453,
 842     0.074023729976459,
 843     0.082381540184189,
 844     0.091120527805195,
 845     0.100243410143811,
 846     0.109751092968147,
 847     0.119642670510092,
 848     0.129915425465314,
 849     0.140564828993259,
 850     0.151584540717153,
 851     0.162966408723997,
 852     0.174700469564574,
 853     0.186774948253444,
 854     0.199176258268946,
 855     0.211889001553197,
 856     0.224895968512091,
 857     0.238178138015305,
 858     0.251714677396289,
 859     0.265482942452275,
 860     0.279458477444273,
 861     0.312726362409309,
 862     0.346648754292945,
 863     0.380914974633193,
 864     0.415237358187469,
 865     0.449351252932597,
 866     0.483015020064806,
 867     0.516010033999735,
 868     0.548140682372425,
 869     0.579234366037328,
 870     0.609141499068300,
 871     0.637735508758604,
 872     0.664912835620912,
 873     0.690592933387298,
 874     0.714718269009247,
 875     0.737254322657649,
 876     0.758189587722801,
 877     0.777535570814405,
 878     0.795326791761572,
 879     0.811620783612819,
 880     0.826498092636068,
 881     0.840062278318649,
 882     0.852439913367300,
 883     0.863780583708163,
 884     0.874256888486789,
 885     0.884064440068133,
 886     0.893421864036559,
 887     0.902570799195836,
 888     0.911775897569142,
 889     0.921324824399059,
 890     0.931528258147577
 891   };
 892
 893   xy_table data_table = make_xy_table(data_x, data_y, 3);
 894   xy_table test_table = make_xy_table(test_x, test_y, 61);
 895   xy_table test_d_table = make_xy_table(test_x, test_dy, 61);
 896   xy_table test_i_table = make_xy_table(test_x, test_iy, 61);
 897
 898   s = test_interp (&data_table, gsl_interp_cspline_periodic, &test_table, &test_d_table, &test_i_table);
 899   gsl_test (s, "cspline periodic 3pt interpolation");
 900   return s;
 901 }
 902
 903
 904
 905 static int
 906 test_akima (void)
 907 {
 908   int s;
 909
 910   double data_x[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
 911   double data_y[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
 912   double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
 913   double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
 914   double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
 915   double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };
 916
 917   xy_table data_table = make_xy_table(data_x, data_y, 5);
 918   xy_table test_table = make_xy_table(test_x, test_y, 4);
 919   xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
 920   xy_table test_i_table = make_xy_table(test_x, test_iy, 4);
 921
 922   s = test_interp (&data_table, gsl_interp_akima, &test_table, &test_d_table, &test_i_table);
 923   gsl_test (s, "akima interpolation");
 924   return s;
 925 }
 926
 927
 928 int
 929 main (int argc, char **argv)
 930 {
 931   int status = 0;
 932
 933   gsl_ieee_env_setup ();
 934
 935   argc = 0;    /* prevent warnings about unused parameters */
 936   argv = 0;
 937
 938   status += test_bsearch();
 939   status += test_linear();
 940   status += test_polynomial();
 941   status += test_cspline();
 942   status += test_cspline2();
 943   status += test_cspline3();
 944   status += test_csplinep();
 945   status += test_csplinep2();
 946   status += test_akima();
 947
 948   exit (gsl_test_summary());
 949 }