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

   1 /* specfunc/mathieu_angfunc.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 <stdio.h>
  25 #include <math.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_sf_mathieu.h>
  28
  29
  30 int gsl_sf_mathieu_ce(int order, double qq, double zz, gsl_sf_result *result)
  31 {
  32   int even_odd, ii, status;
  33   double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
  34   gsl_sf_result aa;
  35
  36
  37   norm = 0.0;
  38   even_odd = 0;
  39   if (order % 2 != 0)
  40       even_odd = 1;
  41
  42   /* Handle the trivial case where q = 0. */
  43   if (qq == 0.0)
  44   {
  45       norm = 1.0;
  46       if (order == 0)
  47           norm = sqrt(2.0);
  48
  49       fn = cos(order*zz)/norm;
  50
  51       result->val = fn;
  52       result->err = 2.0*GSL_DBL_EPSILON;
  53       factor = fabs(fn);
  54       if (factor > 1.0)
  55           result->err *= factor;
  56
  57       return GSL_SUCCESS;
  58   }
  59
  60   /* Use symmetry characteristics of the functions to handle cases with
  61      negative order. */
  62   if (order < 0)
  63       order *= -1;
  64
  65   /* Compute the characteristic value. */
  66   status = gsl_sf_mathieu_a(order, qq, &aa);
  67   if (status != GSL_SUCCESS)
  68   {
  69       return status;
  70   }
  71
  72   /* Compute the series coefficients. */
  73   status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
  74   if (status != GSL_SUCCESS)
  75   {
  76       return status;
  77   }
  78
  79   if (even_odd == 0)
  80   {
  81       fn = 0.0;
  82       norm = coeff[0]*coeff[0];
  83       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
  84       {
  85           fn += coeff[ii]*cos(2.0*ii*zz);
  86           norm += coeff[ii]*coeff[ii];
  87       }
  88   }
  89   else
  90   {
  91       fn = 0.0;
  92       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
  93       {
  94           fn += coeff[ii]*cos((2.0*ii + 1.0)*zz);
  95           norm += coeff[ii]*coeff[ii];
  96       }
  97   }
  98
  99   norm = sqrt(norm);
 100   fn /= norm;
 101
 102   result->val = fn;
 103   result->err = 2.0*GSL_DBL_EPSILON;
 104   factor = fabs(fn);
 105   if (factor > 1.0)
 106       result->err *= factor;
 107
 108   return GSL_SUCCESS;
 109 }
 110
 111
 112 int gsl_sf_mathieu_se(int order, double qq, double zz, gsl_sf_result *result)
 113 {
 114   int even_odd, ii, status;
 115   double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
 116   gsl_sf_result aa;
 117
 118
 119   norm = 0.0;
 120   even_odd = 0;
 121   if (order % 2 != 0)
 122       even_odd = 1;
 123
 124   /* Handle the trivial cases where order = 0 and/or q = 0. */
 125   if (order == 0)
 126   {
 127       result->val = 0.0;
 128       result->err = 0.0;
 129       return GSL_SUCCESS;
 130   }
 131
 132   if (qq == 0.0)
 133   {
 134       norm = 1.0;
 135       fn = sin(order*zz);
 136
 137       result->val = fn;
 138       result->err = 2.0*GSL_DBL_EPSILON;
 139       factor = fabs(fn);
 140       if (factor > 1.0)
 141           result->err *= factor;
 142
 143       return GSL_SUCCESS;
 144   }
 145
 146   /* Use symmetry characteristics of the functions to handle cases with
 147      negative order. */
 148   if (order < 0)
 149       order *= -1;
 150
 151   /* Compute the characteristic value. */
 152   status = gsl_sf_mathieu_b(order, qq, &aa);
 153   if (status != GSL_SUCCESS)
 154   {
 155       return status;
 156   }
 157
 158   /* Compute the series coefficients. */
 159   status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
 160   if (status != GSL_SUCCESS)
 161   {
 162       return status;
 163   }
 164
 165   if (even_odd == 0)
 166   {
 167       fn = 0.0;
 168       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 169       {
 170           norm += coeff[ii]*coeff[ii];
 171           fn += coeff[ii]*sin(2.0*(ii + 1)*zz);
 172       }
 173   }
 174   else
 175   {
 176       fn = 0.0;
 177       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 178       {
 179           norm += coeff[ii]*coeff[ii];
 180           fn += coeff[ii]*sin((2.0*ii + 1)*zz);
 181       }
 182   }
 183   norm = sqrt(norm);
 184   fn /= norm;
 185
 186   result->val = fn;
 187   result->err = 2.0*GSL_DBL_EPSILON;
 188   factor = fabs(fn);
 189   if (factor > 1.0)
 190       result->err *= factor;
 191
 192   return GSL_SUCCESS;
 193 }
 194
 195
 196 int gsl_sf_mathieu_ce_array(int nmin, int nmax, double qq, double zz,
 197                             gsl_sf_mathieu_workspace *work,
 198                             double result_array[])
 199 {
 200   int even_odd, order, ii, jj, status;
 201   double coeff[GSL_SF_MATHIEU_COEFF], *aa = work->aa, norm;
 202
 203
 204   /* Initialize the result array to zeroes. */
 205   for (ii=0; ii<nmax-nmin+1; ii++)
 206       result_array[ii] = 0.0;
 207
 208   /* Ensure that the workspace is large enough to accomodate. */
 209   if (work->size < (unsigned int)nmax)
 210   {
 211       GSL_ERROR("Work space not large enough", GSL_EINVAL);
 212   }
 213
 214   if (nmin < 0 || nmax < nmin)
 215   {
 216       GSL_ERROR("domain error", GSL_EDOM);
 217   }
 218
 219   /* Compute all of the eigenvalues up to nmax. */
 220   gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
 221
 222   for (ii=0, order=nmin; order<=nmax; ii++, order++)
 223   {
 224       norm = 0.0;
 225       even_odd = 0;
 226       if (order % 2 != 0)
 227           even_odd = 1;
 228
 229       /* Handle the trivial case where q = 0. */
 230       if (qq == 0.0)
 231       {
 232           norm = 1.0;
 233           if (order == 0)
 234               norm = sqrt(2.0);
 235
 236           result_array[ii] = cos(order*zz)/norm;
 237
 238           continue;
 239       }
 240
 241       /* Compute the series coefficients. */
 242       status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
 243       if (status != GSL_SUCCESS)
 244           return status;
 245
 246       if (even_odd == 0)
 247       {
 248           norm = coeff[0]*coeff[0];
 249           for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
 250           {
 251               result_array[ii] += coeff[jj]*cos(2.0*jj*zz);
 252               norm += coeff[jj]*coeff[jj];
 253           }
 254       }
 255       else
 256       {
 257           for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
 258           {
 259               result_array[ii] += coeff[jj]*cos((2.0*jj + 1.0)*zz);
 260               norm += coeff[jj]*coeff[jj];
 261           }
 262       }
 263
 264       norm = sqrt(norm);
 265       result_array[ii] /= norm;
 266   }
 267
 268   return GSL_SUCCESS;
 269 }
 270
 271
 272 int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,
 273                             gsl_sf_mathieu_workspace *work,
 274                             double result_array[])
 275 {
 276   int even_odd, order, ii, jj, status;
 277   double coeff[GSL_SF_MATHIEU_COEFF], *bb = work->bb, norm;
 278
 279
 280   /* Initialize the result array to zeroes. */
 281   for (ii=0; ii<nmax-nmin+1; ii++)
 282       result_array[ii] = 0.0;
 283
 284   /* Ensure that the workspace is large enough to accomodate. */
 285   if (work->size < (unsigned int)nmax)
 286   {
 287       GSL_ERROR("Work space not large enough", GSL_EINVAL);
 288   }
 289
 290   if (nmin < 0 || nmax < nmin)
 291   {
 292       GSL_ERROR("domain error", GSL_EDOM);
 293   }
 294
 295   /* Compute all of the eigenvalues up to nmax. */
 296   gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
 297
 298   for (ii=0, order=nmin; order<=nmax; ii++, order++)
 299   {
 300       norm = 0.0;
 301       even_odd = 0;
 302       if (order % 2 != 0)
 303           even_odd = 1;
 304
 305       /* Handle the trivial case where q = 0. */
 306       if (qq == 0.0)
 307       {
 308           norm = 1.0;
 309           result_array[ii] = sin(order*zz);
 310
 311           continue;
 312       }
 313
 314       /* Compute the series coefficients. */
 315       status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
 316       if (status != GSL_SUCCESS)
 317       {
 318           return status;
 319       }
 320
 321       if (even_odd == 0)
 322       {
 323           for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
 324           {
 325               result_array[ii] += coeff[jj]*sin(2.0*(jj + 1)*zz);
 326               norm += coeff[jj]*coeff[jj];
 327           }
 328       }
 329       else
 330       {
 331           for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
 332           {
 333               result_array[ii] += coeff[jj]*sin((2.0*jj + 1.0)*zz);
 334               norm += coeff[jj]*coeff[jj];
 335           }
 336       }
 337
 338       norm = sqrt(norm);
 339       result_array[ii] /= norm;
 340   }
 341
 342   return GSL_SUCCESS;
 343 }