/* specfunc/mathieu_radfunc.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 <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_mathieu.h>


int gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz,
                      gsl_sf_result *result)
{
  int even_odd, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1pc, z2pc;
  double u1, u2;
  gsl_sf_result aa;


  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);
  
  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 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)
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
          }
              
          fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
          fn += fc*j1c*z2c;
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }
  else
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+1, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+1, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+1, u2);
          }
          fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
          fn += fc*(j1c*z2pc + j1pc*z2c);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }

  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_Ms(int kind, int order, double qq, double zz,
                      gsl_sf_result *result)
{
  int even_odd, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
  double u1, u2;
  gsl_sf_result aa;


  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);
  
  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 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)
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1mc = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+2, u1);
          if (kind == 1)
          {
              z2mc = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+2, u2);
          }
          else /* kind = 2 */
          {
              z2mc = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+2, u2);
          }
          
          fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
          fn += fc*(j1mc*z2pc - j1pc*z2mc);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }
  else
  {
      for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
      {
          amax = GSL_MAX(amax, fabs(coeff[kk]));
          if (fabs(coeff[kk])/amax < maxerr)
              break;

          j1c = gsl_sf_bessel_Jn(kk, u1);
          j1pc = gsl_sf_bessel_Jn(kk+1, u1);
          if (kind == 1)
          {
              z2c = gsl_sf_bessel_Jn(kk, u2);
              z2pc = gsl_sf_bessel_Jn(kk+1, u2);
          }
          else /* kind = 2 */
          {
              z2c = gsl_sf_bessel_Yn(kk, u2);
              z2pc = gsl_sf_bessel_Yn(kk+1, u2);
          }
          
          fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
          fn += fc*(j1c*z2pc - j1pc*z2c);
      }

      fn *= sqrt(pi/2.0)/coeff[0];
  }

  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_Mc_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1pc, z2pc;
  double u1, u2;
  double *aa = work->aa;


  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;
  
  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);
  
  /* Compute all eigenvalues up to nmax. */
  gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
  
  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;

      /* 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)
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
              }
              
              fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
              fn += fc*j1c*z2c;
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }
      else
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+1, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+1, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+1, u2);
              }
              fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
              fn += fc*(j1c*z2pc + j1pc*z2c);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }

      result_array[ii] = fn;
  } /* order loop */
  
  return GSL_SUCCESS;
}


int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[])
{
  int even_odd, order, ii, kk, mm, status;
  double maxerr = 1e-14, amax, pi = M_PI, fn;
  double coeff[GSL_SF_MATHIEU_COEFF], fc;
  double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
  double u1, u2;
  double *bb = work->bb;


  /* Initialize the result array to zeroes. */
  for (ii=0; ii<nmax-nmin+1; ii++)
      result_array[ii] = 0.0;
  
  /* Check for out of bounds parameters. */
  if (qq <= 0.0)
  {
      GSL_ERROR("q must be greater than zero", GSL_EINVAL);
  }
  if (kind < 1 || kind > 2)
  {
      GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
  }

  mm = 0;
  amax = 0.0;
  fn = 0.0;
  u1 = sqrt(qq)*exp(-1.0*zz);
  u2 = sqrt(qq)*exp(zz);
  
  /* Compute all eigenvalues up to nmax. */
  gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
  
  for (ii=0, order=nmin; order<=nmax; ii++, order++)
  {
      even_odd = 0;
      if (order % 2 != 0)
          even_odd = 1;
  
      /* 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 (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1mc = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+2, u1);
              if (kind == 1)
              {
                  z2mc = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+2, u2);
              }
              else /* kind = 2 */
              {
                  z2mc = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+2, u2);
              }
          
              fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
              fn += fc*(j1mc*z2pc - j1pc*z2mc);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }
      else
      {
          for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
          {
              amax = GSL_MAX(amax, fabs(coeff[kk]));
              if (fabs(coeff[kk])/amax < maxerr)
                  break;

              j1c = gsl_sf_bessel_Jn(kk, u1);
              j1pc = gsl_sf_bessel_Jn(kk+1, u1);
              if (kind == 1)
              {
                  z2c = gsl_sf_bessel_Jn(kk, u2);
                  z2pc = gsl_sf_bessel_Jn(kk+1, u2);
              }
              else /* kind = 2 */
              {
                  z2c = gsl_sf_bessel_Yn(kk, u2);
                  z2pc = gsl_sf_bessel_Yn(kk+1, u2);
              }
          
              fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
              fn += fc*(j1c*z2pc - j1pc*z2c);
          }

          fn *= sqrt(pi/2.0)/coeff[0];
      }

      result_array[ii] = fn;
  } /* order loop */
  
  return GSL_SUCCESS;
}