/* specfunc/mathieu_angfunc.c

 * 

 * Copyright (C) 2002 Lowell Johnson

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

 */

/* Author:  L. Johnson */

#include <config.h>

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_sf_mathieu.h>

int gsl_sf_mathieu_ce_e(int order, double qq, double zz, gsl_sf_result *result)

{

  int even_odd, ii, status;

  double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;

  gsl_sf_result aa;

  norm = 0.0;

  even_odd = 0;

  if (order % 2 != 0)

      even_odd = 1;

  /* Handle the trivial case where q = 0. */

  if (qq == 0.0)

  {

      norm = 1.0;

      if (order == 0)

          norm = sqrt(2.0);

      fn = cos(order*zz)/norm;

      result->val = fn;

      result->err = 2.0*GSL_DBL_EPSILON;

      factor = fabs(fn);

      if (factor > 1.0)

          result->err *= factor;

      return GSL_SUCCESS;

  }

  /* Use symmetry characteristics of the functions to handle cases with

     negative order. */

  if (order < 0)

      order *= -1;

  /* Compute the characteristic value. */

  status = gsl_sf_mathieu_a_e(order, qq, &aa);

  if (status != GSL_SUCCESS)

  {

      return status;

  }

  /* Compute the series coefficients. */

  status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);

  if (status != GSL_SUCCESS)

  {

      return status;

  }

  if (even_odd == 0)

  {

      fn = 0.0;

      norm = coeff[0]*coeff[0];

      for (ii=0; ii

      {

          fn += coeff[ii]*cos(2.0*ii*zz);

          norm += coeff[ii]*coeff[ii];

      }

  }

  else

  {

      fn = 0.0;

      for (ii=0; ii

      {

          fn += coeff[ii]*cos((2.0*ii + 1.0)*zz);

          norm += coeff[ii]*coeff[ii];

      }

  }

  norm = sqrt(norm);

  fn /= norm;

  result->val = fn;

  result->err = 2.0*GSL_DBL_EPSILON;

  factor = fabs(fn);

  if (factor > 1.0)

      result->err *= factor;

  return GSL_SUCCESS;

}

int gsl_sf_mathieu_se_e(int order, double qq, double zz, gsl_sf_result *result)

{

  int even_odd, ii, status;

  double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;

  gsl_sf_result aa;

  norm = 0.0;

  even_odd = 0;

  if (order % 2 != 0)

      even_odd = 1;

  /* Handle the trivial cases where order = 0 and/or q = 0. */

  if (order == 0)

  {

      result->val = 0.0;

      result->err = 0.0;

      return GSL_SUCCESS;

  }

  if (qq == 0.0)

  {

      norm = 1.0;

      fn = sin(order*zz);

      result->val = fn;

      result->err = 2.0*GSL_DBL_EPSILON;

      factor = fabs(fn);

      if (factor > 1.0)

          result->err *= factor;

      return GSL_SUCCESS;

  }

  /* Use symmetry characteristics of the functions to handle cases with

     negative order. */

  if (order < 0)

      order *= -1;

  /* Compute the characteristic value. */

  status = gsl_sf_mathieu_b_e(order, qq, &aa);

  if (status != GSL_SUCCESS)

  {

      return status;

  }

  /* Compute the series coefficients. */

  status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);

  if (status != GSL_SUCCESS)

  {

      return status;

  }

  if (even_odd == 0)

  {

      fn = 0.0;

      for (ii=0; ii

      {

          norm += coeff[ii]*coeff[ii];

          fn += coeff[ii]*sin(2.0*(ii + 1)*zz);

      }

  }

  else

  {

      fn = 0.0;

      for (ii=0; ii

      {

          norm += coeff[ii]*coeff[ii];

          fn += coeff[ii]*sin((2.0*ii + 1)*zz);

      }

  }

  norm = sqrt(norm);

  fn /= norm;

  result->val = fn;

  result->err = 2.0*GSL_DBL_EPSILON;

  factor = fabs(fn);

  if (factor > 1.0)

      result->err *= factor;

  return GSL_SUCCESS;

}

int gsl_sf_mathieu_ce_array(int nmin, int nmax, double qq, double zz,

                            gsl_sf_mathieu_workspace *work,

                            double result_array[])

{

  int even_odd, order, ii, jj, status;

  double coeff[GSL_SF_MATHIEU_COEFF], *aa = work->aa, norm;

  /* Initialize the result array to zeroes. */

  for (ii=0; ii

      result_array[ii] = 0.0;

  /* Ensure that the workspace is large enough to accomodate. */

  if (work->size < (unsigned int)nmax)

  {

      GSL_ERROR("Work space not large enough", GSL_EINVAL);

  }

  if (nmin < 0 || nmax < nmin)

  {

      GSL_ERROR("domain error", GSL_EDOM);

  }

  /* Compute all of the eigenvalues up to nmax. */

  gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)

  {

      norm = 0.0;

      even_odd = 0;

      if (order % 2 != 0)

          even_odd = 1;

      /* Handle the trivial case where q = 0. */

      if (qq == 0.0)

      {

          norm = 1.0;

          if (order == 0)

              norm = sqrt(2.0);

          result_array[ii] = cos(order*zz)/norm;

          continue;

      }

      /* Compute the series coefficients. */

      status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);

      if (status != GSL_SUCCESS)

          return status;

      if (even_odd == 0)

      {

          norm = coeff[0]*coeff[0];

          for (jj=0; jj

          {

              result_array[ii] += coeff[jj]*cos(2.0*jj*zz);

              norm += coeff[jj]*coeff[jj];

          }

      }

      else

      {

          for (jj=0; jj

          {

              result_array[ii] += coeff[jj]*cos((2.0*jj + 1.0)*zz);

              norm += coeff[jj]*coeff[jj];

          }

      }

      norm = sqrt(norm);

      result_array[ii] /= norm;

  }

  return GSL_SUCCESS;

}

int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,

                            gsl_sf_mathieu_workspace *work,

                            double result_array[])

{

  int even_odd, order, ii, jj, status;

  double coeff[GSL_SF_MATHIEU_COEFF], *bb = work->bb, norm;

  /* Initialize the result array to zeroes. */

  for (ii=0; ii

      result_array[ii] = 0.0;

  /* Ensure that the workspace is large enough to accomodate. */

  if (work->size < (unsigned int)nmax)

  {

      GSL_ERROR("Work space not large enough", GSL_EINVAL);

  }

  if (nmin < 0 || nmax < nmin)

  {

      GSL_ERROR("domain error", GSL_EDOM);

  }

  /* Compute all of the eigenvalues up to nmax. */

  gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);

  for (ii=0, order=nmin; order<=nmax; ii++, order++)

  {

      norm = 0.0;

      even_odd = 0;

      if (order % 2 != 0)

          even_odd = 1;

      /* Handle the trivial cases where order = 0 and/or q = 0. */

      if (order == 0)

      {

          norm = 1.0;

          result_array[ii] = 0.0;

          continue;

      }

      if (qq == 0.0)

      {

          norm = 1.0;

          result_array[ii] = sin(order*zz);

          continue;

      }

      /* Compute the series coefficients. */

      status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);

      if (status != GSL_SUCCESS)

      {

          return status;

      }

      if (even_odd == 0)

      {

          for (jj=0; jj

          {

              result_array[ii] += coeff[jj]*sin(2.0*(jj + 1)*zz);

              norm += coeff[jj]*coeff[jj];

          }

      }

      else

      {

          for (jj=0; jj

          {

              result_array[ii] += coeff[jj]*sin((2.0*jj + 1.0)*zz);

              norm += coeff[jj]*coeff[jj];

          }

      }

      norm = sqrt(norm);

      result_array[ii] /= norm;

  }

  return GSL_SUCCESS;

}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"                                                           

double gsl_sf_mathieu_ce(int order, double qq, double zz)

{

	EVAL_RESULT(gsl_sf_mathieu_ce_e(order, qq, zz, &result));

}

double gsl_sf_mathieu_se(int order, double qq, double zz)

{

	EVAL_RESULT(gsl_sf_mathieu_se_e(order, qq, zz, &result));

}