/* 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 2 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(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(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<GSL_SF_MATHIEU_COEFF; ii++)
      {
          fn += coeff[ii]*cos(2.0*ii*zz);
          norm += coeff[ii]*coeff[ii];
      }
  }
  else
  {
      fn = 0.0;
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; 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(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(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<GSL_SF_MATHIEU_COEFF; ii++)
      {
          norm += coeff[ii]*coeff[ii];
          fn += coeff[ii]*sin(2.0*(ii + 1)*zz);
      }
  }
  else
  {
      fn = 0.0;
      for (ii=0; ii<GSL_SF_MATHIEU_COEFF; 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<nmax-nmin+1; 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<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*cos(2.0*jj*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }
      else
      {
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; 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<nmax-nmin+1; 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 case where q = 0. */
      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<GSL_SF_MATHIEU_COEFF; jj++)
          {
              result_array[ii] += coeff[jj]*sin(2.0*(jj + 1)*zz);
              norm += coeff[jj]*coeff[jj];
          }
      }
      else
      {
          for (jj=0; jj<GSL_SF_MATHIEU_COEFF; 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;
}