htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/mathieu_radfunc.c

   1 /* specfunc/mathieu_radfunc.c
   2  *
   3  * Copyright (C) 2002 Lowell Johnson
   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., 675 Mass Ave, Cambridge, MA 02139, USA.
  18  */
  19
  20 /* Author:  L. Johnson */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_sf.h>
  27 #include <gsl/gsl_sf_mathieu.h>
  28
  29
  30 int gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz,
  31                       gsl_sf_result *result)
  32 {
  33   int even_odd, kk, mm, status;
  34   double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
  35   double coeff[GSL_SF_MATHIEU_COEFF], fc;
  36   double j1c, z2c, j1pc, z2pc;
  37   double u1, u2;
  38   gsl_sf_result aa;
  39
  40
  41   /* Check for out of bounds parameters. */
  42   if (qq <= 0.0)
  43   {
  44       GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  45   }
  46   if (kind < 1 || kind > 2)
  47   {
  48       GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  49   }
  50
  51   mm = 0;
  52   amax = 0.0;
  53   fn = 0.0;
  54   u1 = sqrt(qq)*exp(-1.0*zz);
  55   u2 = sqrt(qq)*exp(zz);
  56
  57   even_odd = 0;
  58   if (order % 2 != 0)
  59       even_odd = 1;
  60
  61   /* Compute the characteristic value. */
  62   status = gsl_sf_mathieu_a(order, qq, &aa);
  63   if (status != GSL_SUCCESS)
  64   {
  65       return status;
  66   }
  67
  68   /* Compute the series coefficients. */
  69   status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
  70   if (status != GSL_SUCCESS)
  71   {
  72       return status;
  73   }
  74
  75   if (even_odd == 0)
  76   {
  77       for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
  78       {
  79           amax = GSL_MAX(amax, fabs(coeff[kk]));
  80           if (fabs(coeff[kk])/amax < maxerr)
  81               break;
  82
  83           j1c = gsl_sf_bessel_Jn(kk, u1);
  84           if (kind == 1)
  85           {
  86               z2c = gsl_sf_bessel_Jn(kk, u2);
  87           }
  88           else /* kind = 2 */
  89           {
  90               z2c = gsl_sf_bessel_Yn(kk, u2);
  91           }
  92
  93           fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
  94           fn += fc*j1c*z2c;
  95       }
  96
  97       fn *= sqrt(pi/2.0)/coeff[0];
  98   }
  99   else
 100   {
 101       for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 102       {
 103           amax = GSL_MAX(amax, fabs(coeff[kk]));
 104           if (fabs(coeff[kk])/amax < maxerr)
 105               break;
 106
 107           j1c = gsl_sf_bessel_Jn(kk, u1);
 108           j1pc = gsl_sf_bessel_Jn(kk+1, u1);
 109           if (kind == 1)
 110           {
 111               z2c = gsl_sf_bessel_Jn(kk, u2);
 112               z2pc = gsl_sf_bessel_Jn(kk+1, u2);
 113           }
 114           else /* kind = 2 */
 115           {
 116               z2c = gsl_sf_bessel_Yn(kk, u2);
 117               z2pc = gsl_sf_bessel_Yn(kk+1, u2);
 118           }
 119           fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
 120           fn += fc*(j1c*z2pc + j1pc*z2c);
 121       }
 122
 123       fn *= sqrt(pi/2.0)/coeff[0];
 124   }
 125
 126   result->val = fn;
 127   result->err = 2.0*GSL_DBL_EPSILON;
 128   factor = fabs(fn);
 129   if (factor > 1.0)
 130       result->err *= factor;
 131
 132   return GSL_SUCCESS;
 133 }
 134
 135
 136 int gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz,
 137                       gsl_sf_result *result)
 138 {
 139   int even_odd, kk, mm, status;
 140   double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
 141   double coeff[GSL_SF_MATHIEU_COEFF], fc;
 142   double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
 143   double u1, u2;
 144   gsl_sf_result aa;
 145
 146
 147   /* Check for out of bounds parameters. */
 148   if (qq <= 0.0)
 149   {
 150       GSL_ERROR("q must be greater than zero", GSL_EINVAL);
 151   }
 152   if (kind < 1 || kind > 2)
 153   {
 154       GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
 155   }
 156
 157   mm = 0;
 158   amax = 0.0;
 159   fn = 0.0;
 160   u1 = sqrt(qq)*exp(-1.0*zz);
 161   u2 = sqrt(qq)*exp(zz);
 162
 163   even_odd = 0;
 164   if (order % 2 != 0)
 165       even_odd = 1;
 166
 167   /* Compute the characteristic value. */
 168   status = gsl_sf_mathieu_b(order, qq, &aa);
 169   if (status != GSL_SUCCESS)
 170   {
 171       return status;
 172   }
 173
 174   /* Compute the series coefficients. */
 175   status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
 176   if (status != GSL_SUCCESS)
 177   {
 178       return status;
 179   }
 180
 181   if (even_odd == 0)
 182   {
 183       for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 184       {
 185           amax = GSL_MAX(amax, fabs(coeff[kk]));
 186           if (fabs(coeff[kk])/amax < maxerr)
 187               break;
 188
 189           j1mc = gsl_sf_bessel_Jn(kk, u1);
 190           j1pc = gsl_sf_bessel_Jn(kk+2, u1);
 191           if (kind == 1)
 192           {
 193               z2mc = gsl_sf_bessel_Jn(kk, u2);
 194               z2pc = gsl_sf_bessel_Jn(kk+2, u2);
 195           }
 196           else /* kind = 2 */
 197           {
 198               z2mc = gsl_sf_bessel_Yn(kk, u2);
 199               z2pc = gsl_sf_bessel_Yn(kk+2, u2);
 200           }
 201
 202           fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
 203           fn += fc*(j1mc*z2pc - j1pc*z2mc);
 204       }
 205
 206       fn *= sqrt(pi/2.0)/coeff[0];
 207   }
 208   else
 209   {
 210       for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 211       {
 212           amax = GSL_MAX(amax, fabs(coeff[kk]));
 213           if (fabs(coeff[kk])/amax < maxerr)
 214               break;
 215
 216           j1c = gsl_sf_bessel_Jn(kk, u1);
 217           j1pc = gsl_sf_bessel_Jn(kk+1, u1);
 218           if (kind == 1)
 219           {
 220               z2c = gsl_sf_bessel_Jn(kk, u2);
 221               z2pc = gsl_sf_bessel_Jn(kk+1, u2);
 222           }
 223           else /* kind = 2 */
 224           {
 225               z2c = gsl_sf_bessel_Yn(kk, u2);
 226               z2pc = gsl_sf_bessel_Yn(kk+1, u2);
 227           }
 228
 229           fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
 230           fn += fc*(j1c*z2pc - j1pc*z2c);
 231       }
 232
 233       fn *= sqrt(pi/2.0)/coeff[0];
 234   }
 235
 236   result->val = fn;
 237   result->err = 2.0*GSL_DBL_EPSILON;
 238   factor = fabs(fn);
 239   if (factor > 1.0)
 240       result->err *= factor;
 241
 242   return GSL_SUCCESS;
 243 }
 244
 245
 246 int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
 247                             double zz, gsl_sf_mathieu_workspace *work,
 248                             double result_array[])
 249 {
 250   int even_odd, order, ii, kk, mm, status;
 251   double maxerr = 1e-14, amax, pi = M_PI, fn;
 252   double coeff[GSL_SF_MATHIEU_COEFF], fc;
 253   double j1c, z2c, j1pc, z2pc;
 254   double u1, u2;
 255   double *aa = work->aa;
 256
 257
 258   /* Initialize the result array to zeroes. */
 259   for (ii=0; ii<nmax-nmin+1; ii++)
 260       result_array[ii] = 0.0;
 261
 262   /* Check for out of bounds parameters. */
 263   if (qq <= 0.0)
 264   {
 265       GSL_ERROR("q must be greater than zero", GSL_EINVAL);
 266   }
 267   if (kind < 1 || kind > 2)
 268   {
 269       GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
 270   }
 271
 272   mm = 0;
 273   amax = 0.0;
 274   fn = 0.0;
 275   u1 = sqrt(qq)*exp(-1.0*zz);
 276   u2 = sqrt(qq)*exp(zz);
 277
 278   /* Compute all eigenvalues up to nmax. */
 279   gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
 280
 281   for (ii=0, order=nmin; order<=nmax; ii++, order++)
 282   {
 283       even_odd = 0;
 284       if (order % 2 != 0)
 285           even_odd = 1;
 286
 287       /* Compute the series coefficients. */
 288       status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
 289       if (status != GSL_SUCCESS)
 290       {
 291           return status;
 292       }
 293
 294       if (even_odd == 0)
 295       {
 296           for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 297           {
 298               amax = GSL_MAX(amax, fabs(coeff[kk]));
 299               if (fabs(coeff[kk])/amax < maxerr)
 300                   break;
 301
 302               j1c = gsl_sf_bessel_Jn(kk, u1);
 303               if (kind == 1)
 304               {
 305                   z2c = gsl_sf_bessel_Jn(kk, u2);
 306               }
 307               else /* kind = 2 */
 308               {
 309                   z2c = gsl_sf_bessel_Yn(kk, u2);
 310               }
 311
 312               fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
 313               fn += fc*j1c*z2c;
 314           }
 315
 316           fn *= sqrt(pi/2.0)/coeff[0];
 317       }
 318       else
 319       {
 320           for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 321           {
 322               amax = GSL_MAX(amax, fabs(coeff[kk]));
 323               if (fabs(coeff[kk])/amax < maxerr)
 324                   break;
 325
 326               j1c = gsl_sf_bessel_Jn(kk, u1);
 327               j1pc = gsl_sf_bessel_Jn(kk+1, u1);
 328               if (kind == 1)
 329               {
 330                   z2c = gsl_sf_bessel_Jn(kk, u2);
 331                   z2pc = gsl_sf_bessel_Jn(kk+1, u2);
 332               }
 333               else /* kind = 2 */
 334               {
 335                   z2c = gsl_sf_bessel_Yn(kk, u2);
 336                   z2pc = gsl_sf_bessel_Yn(kk+1, u2);
 337               }
 338               fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
 339               fn += fc*(j1c*z2pc + j1pc*z2c);
 340           }
 341
 342           fn *= sqrt(pi/2.0)/coeff[0];
 343       }
 344
 345       result_array[ii] = fn;
 346   } /* order loop */
 347
 348   return GSL_SUCCESS;
 349 }
 350
 351
 352 int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
 353                             double zz, gsl_sf_mathieu_workspace *work,
 354                             double result_array[])
 355 {
 356   int even_odd, order, ii, kk, mm, status;
 357   double maxerr = 1e-14, amax, pi = M_PI, fn;
 358   double coeff[GSL_SF_MATHIEU_COEFF], fc;
 359   double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
 360   double u1, u2;
 361   double *bb = work->bb;
 362
 363
 364   /* Initialize the result array to zeroes. */
 365   for (ii=0; ii<nmax-nmin+1; ii++)
 366       result_array[ii] = 0.0;
 367
 368   /* Check for out of bounds parameters. */
 369   if (qq <= 0.0)
 370   {
 371       GSL_ERROR("q must be greater than zero", GSL_EINVAL);
 372   }
 373   if (kind < 1 || kind > 2)
 374   {
 375       GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
 376   }
 377
 378   mm = 0;
 379   amax = 0.0;
 380   fn = 0.0;
 381   u1 = sqrt(qq)*exp(-1.0*zz);
 382   u2 = sqrt(qq)*exp(zz);
 383
 384   /* Compute all eigenvalues up to nmax. */
 385   gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
 386
 387   for (ii=0, order=nmin; order<=nmax; ii++, order++)
 388   {
 389       even_odd = 0;
 390       if (order % 2 != 0)
 391           even_odd = 1;
 392
 393       /* Compute the series coefficients. */
 394       status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
 395       if (status != GSL_SUCCESS)
 396       {
 397           return status;
 398       }
 399
 400       if (even_odd == 0)
 401       {
 402           for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 403           {
 404               amax = GSL_MAX(amax, fabs(coeff[kk]));
 405               if (fabs(coeff[kk])/amax < maxerr)
 406                   break;
 407
 408               j1mc = gsl_sf_bessel_Jn(kk, u1);
 409               j1pc = gsl_sf_bessel_Jn(kk+2, u1);
 410               if (kind == 1)
 411               {
 412                   z2mc = gsl_sf_bessel_Jn(kk, u2);
 413                   z2pc = gsl_sf_bessel_Jn(kk+2, u2);
 414               }
 415               else /* kind = 2 */
 416               {
 417                   z2mc = gsl_sf_bessel_Yn(kk, u2);
 418                   z2pc = gsl_sf_bessel_Yn(kk+2, u2);
 419               }
 420
 421               fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
 422               fn += fc*(j1mc*z2pc - j1pc*z2mc);
 423           }
 424
 425           fn *= sqrt(pi/2.0)/coeff[0];
 426       }
 427       else
 428       {
 429           for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
 430           {
 431               amax = GSL_MAX(amax, fabs(coeff[kk]));
 432               if (fabs(coeff[kk])/amax < maxerr)
 433                   break;
 434
 435               j1c = gsl_sf_bessel_Jn(kk, u1);
 436               j1pc = gsl_sf_bessel_Jn(kk+1, u1);
 437               if (kind == 1)
 438               {
 439                   z2c = gsl_sf_bessel_Jn(kk, u2);
 440                   z2pc = gsl_sf_bessel_Jn(kk+1, u2);
 441               }
 442               else /* kind = 2 */
 443               {
 444                   z2c = gsl_sf_bessel_Yn(kk, u2);
 445                   z2pc = gsl_sf_bessel_Yn(kk+1, u2);
 446               }
 447
 448               fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
 449               fn += fc*(j1c*z2pc - j1pc*z2c);
 450           }
 451
 452           fn *= sqrt(pi/2.0)/coeff[0];
 453       }
 454
 455       result_array[ii] = fn;
 456   } /* order loop */
 457
 458   return GSL_SUCCESS;
 459 }