/* specfunc/mathieu_charv.c

 * 

 * Copyright (C) 2002, 2009 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_eigen.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_mathieu.h>

/* prototypes */

static double solve_cubic(double c2, double c1, double c0);

static double ceer(int order, double qq, double aa, int nterms)

{

  double term, term1;

  int ii, n1;

  if (order == 0)

      term = 0.0;

  else

  {      

      term = 2.0*qq*qq/aa;

      if (order != 2)

      {

          n1 = order/2 - 1;

          for (ii=0; ii

              term = qq*qq/(aa - 4.0*(ii+1)*(ii+1) - term);

      }

  }

  term += order*order;

  term1 = 0.0;

  for (ii=0; ii

      term1 = qq*qq/

        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  if (order == 0)

      term1 *= 2.0;

  return (term + term1 - aa);

}

static double ceor(int order, double qq, double aa, int nterms)

{

  double term, term1;

  int ii, n1;

  term = qq;

  n1 = (int)((float)order/2.0 - 0.5);

  for (ii=0; ii

      term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);

  term += order*order;

  term1 = 0.0;

  for (ii=0; ii

      term1 = qq*qq/

        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);

}

static double seer(int order, double qq, double aa, int nterms)

{

  double term, term1;

  int ii, n1;

  term = 0.0;

  n1 = order/2 - 1;

  for (ii=0; ii

      term = qq*qq/(aa - 4.0*(ii + 1)*(ii + 1) - term);

  term += order*order;

  term1 = 0.0;

  for (ii=0; ii

      term1 = qq*qq/

        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);

}

static double seor(int order, double qq, double aa, int nterms)

{

  double term, term1;

  int ii, n1;

  term = -1.0*qq;

  n1 = (int)((float)order/2.0 - 0.5);

  for (ii=0; ii

      term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);

  term += order*order;

  term1 = 0.0;

  for (ii=0; ii

      term1 = qq*qq/

        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);

}

/*----------------------------------------------------------------------------

 * Asymptotic and approximation routines for the characteristic value.

 *

 * Adapted from F.A. Alhargan's paper,

 * "Algorithms for the Computation of All Mathieu Functions of Integer

 * Orders," ACM Transactions on Mathematical Software, Vol. 26, No. 3,

 * September 2000, pp. 390-407.

 *--------------------------------------------------------------------------*/

static double asymptotic(int order, double qq)

{

  double asymp;

  double nn, n2, n4, n6;

  double hh, ah, ah2, ah3, ah4, ah5;

  /* Set up temporary variables to simplify the readability. */

  nn = 2*order + 1;

  n2 = nn*nn;

  n4 = n2*n2;

  n6 = n4*n2;

  hh = 2*sqrt(qq);

  ah = 16*hh;

  ah2 = ah*ah;

  ah3 = ah2*ah;

  ah4 = ah3*ah;

  ah5 = ah4*ah;

  /* Equation 38, p. 397 of Alhargan's paper. */

  asymp = -2*qq + nn*hh - 0.125*(n2 + 1);

  asymp -= 0.25*nn*(                          n2 +     3)/ah;

  asymp -= 0.25*   (             5*n4 +    34*n2 +     9)/ah2;

  asymp -= 0.25*nn*(            33*n4 +   410*n2 +   405)/ah3;

  asymp -=         ( 63*n6 +  1260*n4 +  2943*n2 +   486)/ah4;

  asymp -=      nn*(527*n6 + 15617*n4 + 69001*n2 + 41607)/ah5;

  return asymp;

}

/* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */

static double solve_cubic(double c2, double c1, double c0)

{

  double qq, rr, ww, ss, tt;

  qq = (3*c1 - c2*c2)/9;

  rr = (9*c2*c1 - 27*c0 - 2*c2*c2*c2)/54;

  ww = qq*qq*qq + rr*rr;

  if (ww >= 0)

  {

      double t1 = rr + sqrt(ww);

      ss = fabs(t1)/t1*pow(fabs(t1), 1/3.);

      t1 = rr - sqrt(ww);

      tt = fabs(t1)/t1*pow(fabs(t1), 1/3.);

  }

  else

  {

      double theta = acos(rr/sqrt(-qq*qq*qq));

      ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.);

      tt = 0.0;

  }

  return (ss + tt - c2/3);

}

/* Compute an initial approximation for the characteristic value. */

static double approx_c(int order, double qq)

{

  double approx;

  double c0, c1, c2;

  if (order < 0)

  {

    GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);

  }

  switch (order)

  {

      case 0:

          if (qq <= 4)

              return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */

          else

              return asymptotic(order, qq);

          break;

      case 1:

          if (qq <= 4)

              return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */

          else

              return asymptotic(order, qq);

          break;

      case 2:

          if (qq <= 3)

          {

              c2 = -8.0;  /* Eqn. 33 */

              c1 = -48 - 3*qq*qq;

              c0 = 20*qq*qq;

          }

          else

              return asymptotic(order, qq);

          break;

      case 3:

          if (qq <= 6.25)

          {

              c2 = -qq - 8;  /* Eqn. 34 */

              c1 = 16*qq - 128 - 2*qq*qq;

              c0 = qq*qq*(qq + 8);

          }

          else

              return asymptotic(order, qq);

          break;

      default:

          if (order < 70)

          {

              if (1.7*order > 2*sqrt(qq))

              {

                  /* Eqn. 30 */

                  double n2 = (double)(order*order);

                  double n22 = (double)((n2 - 1)*(n2 - 1));

                  double q2 = qq*qq;

                  double q4 = q2*q2;

                  approx = n2 + 0.5*q2/(n2 - 1);

                  approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));

                  approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/

                      (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));

                  if (1.4*order < 2*sqrt(qq))

                  {

                      approx += asymptotic(order, qq);

                      approx *= 0.5;

                  }

              }

              else

                  approx = asymptotic(order, qq);

              return approx;

          }

          else

              return order*order;

  }

  /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */

  approx = solve_cubic(c2, c1, c0);

  if ( approx < 0 && sqrt(qq) > 0.1*order )

      return asymptotic(order-1, qq);

  else

      return (order*order + fabs(approx));

}

static double approx_s(int order, double qq)

{

  double approx;

  double c0, c1, c2;

  if (order < 1)

  {

    GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);

  }

  switch (order)

  {

      case 1:

          if (qq <= 4)

              return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */

          else

              return asymptotic(order-1, qq);

          break;

      case 2:

          if (qq <= 5)

              return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */

          else

              return asymptotic(order-1, qq);

          break;

      case 3:

          if (qq <= 6.25)

          {

              c2 = qq - 8; /* Eqn. 37 */

              c1 = -128 - 16*qq - 2*qq*qq;

              c0 = qq*qq*(8 - qq);

          }

          else

              return asymptotic(order-1, qq);

          break;

      default:

          if (order < 70)

          {

              if (1.7*order > 2*sqrt(qq))

              {

                  /* Eqn. 30 */

                  double n2 = (double)(order*order);

                  double n22 = (double)((n2 - 1)*(n2 - 1));

                  double q2 = qq*qq;

                  double q4 = q2*q2;

                  approx = n2 + 0.5*q2/(n2 - 1);

                  approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));

                  approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/

                      (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));

                  if (1.4*order < 2*sqrt(qq))

                  {

                      approx += asymptotic(order-1, qq);

                      approx *= 0.5;

                  }

              }

              else

                  approx = asymptotic(order-1, qq);

              return approx;

          }

          else

              return order*order;

  }

  /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */

  approx = solve_cubic(c2, c1, c0);

  if ( approx < 0 && sqrt(qq) > 0.1*order )

      return asymptotic(order-1, qq);

  else

      return (order*order + fabs(approx));

}

int gsl_sf_mathieu_a_e(int order, double qq, gsl_sf_result *result)

{

  int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;

  int dir = 0;  /* step direction for new search */

  double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;

  double aa_approx;  /* current approximation for solution */

  even_odd = 0;

  if (order % 2 != 0)

      even_odd = 1;

  /* If the argument is 0, then the coefficient is simply the square of

     the order. */

  if (qq == 0)

  {

      result->val = order*order;

      result->err = 0.0;

      return GSL_SUCCESS;

  }

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

     negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */

  if (order < 0)

      order *= -1;

  if (qq < 0.0)

  {

      if (even_odd == 0)

          return gsl_sf_mathieu_a_e(order, -qq, result);

      else

          return gsl_sf_mathieu_b_e(order, -qq, result);

  }

  /* Compute an initial approximation for the characteristic value. */

  aa_approx = approx_c(order, qq);

  /* Save the original approximation for later comparison. */

  aa_orig = aa = aa_approx;

  /* Loop as long as the final value is not near the approximate value

     (with a max limit to avoid potential infinite loop). */

  while (counter < maxcount)

  {

      a1 = aa + 0.001;

      ii = 0;

      if (even_odd == 0)

          fa1 = ceer(order, qq, a1, nterms);

      else

          fa1 = ceor(order, qq, a1, nterms);

      for (;;)

      {

          if (even_odd == 0)

              fa = ceer(order, qq, aa, nterms);

          else

              fa = ceor(order, qq, aa, nterms);

          a2 = a1;

          a1 = aa;

          if (fa == fa1)

          {

              result->err = GSL_DBL_EPSILON;

              break;

          }

          aa -= (aa - a2)/(fa - fa1)*fa;

          dela = fabs(aa - a2);

          if (dela < GSL_DBL_EPSILON)

          {

              result->err = GSL_DBL_EPSILON;

              break;

          }

          if (ii > 40)

          {

              result->err = dela;

              break;

          }

          fa1 = fa;

          ii++;

      }

      /* If the solution found is not near the original approximation,

         tweak the approximate value, and try again. */

      if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||

          (order > 10 && fabs(aa - aa_orig) > 1.5*order))

      {

          counter++;

          if (counter == maxcount)

          {

              result->err = fabs(aa - aa_orig);

              break;

          }

          if (aa > aa_orig)

          {

              if (dir == 1)

                  da /= 2;

              dir = -1;

          }

          else

          {

              if (dir == -1)

                  da /= 2;

              dir = 1;

          }

          aa_approx += dir*da*counter;

          aa = aa_approx;

          continue;

      }

      else

          break;

  }

  result->val = aa;

  /* If we went through the maximum number of retries and still didn't

     find the solution, let us know. */

  if (counter == maxcount)

  {

      GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);

  }

  return GSL_SUCCESS;

}

int gsl_sf_mathieu_b_e(int order, double qq, gsl_sf_result *result)

{

  int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;

  int dir = 0;  /* step direction for new search */

  double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;

  double aa_approx;  /* current approximation for solution */

  even_odd = 0;

  if (order % 2 != 0)

      even_odd = 1;

  /* The order cannot be 0. */

  if (order == 0)

  {

      GSL_ERROR("Characteristic value undefined for order 0", GSL_EFAILED);

  }

  /* If the argument is 0, then the coefficient is simply the square of

     the order. */

  if (qq == 0)

  {

      result->val = order*order;

      result->err = 0.0;

      return GSL_SUCCESS;

  }

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

     negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */

  if (order < 0)

      order *= -1;

  if (qq < 0.0)

  {

      if (even_odd == 0)

          return gsl_sf_mathieu_b_e(order, -qq, result);

      else

          return gsl_sf_mathieu_a_e(order, -qq, result);

  }

  /* Compute an initial approximation for the characteristic value. */

  aa_approx = approx_s(order, qq);

  /* Save the original approximation for later comparison. */

  aa_orig = aa = aa_approx;

  /* Loop as long as the final value is not near the approximate value

     (with a max limit to avoid potential infinite loop). */

  while (counter < maxcount)

  {

      a1 = aa + 0.001;

      ii = 0;

      if (even_odd == 0)

          fa1 = seer(order, qq, a1, nterms);

      else

          fa1 = seor(order, qq, a1, nterms);

      for (;;)

      {

          if (even_odd == 0)

              fa = seer(order, qq, aa, nterms);

          else

              fa = seor(order, qq, aa, nterms);

          a2 = a1;

          a1 = aa;

          if (fa == fa1)

          {

              result->err = GSL_DBL_EPSILON;

              break;

          }

          aa -= (aa - a2)/(fa - fa1)*fa;

          dela = fabs(aa - a2);

          if (dela < 1e-18)

          {

              result->err = GSL_DBL_EPSILON;

              break;

          }

          if (ii > 40)

          {

              result->err = dela;

              break;

          }

          fa1 = fa;

          ii++;

      }

      /* If the solution found is not near the original approximation,

         tweak the approximate value, and try again. */

      if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||

          (order > 10 && fabs(aa - aa_orig) > 1.5*order))

      {

          counter++;

          if (counter == maxcount)

          {

              result->err = fabs(aa - aa_orig);

              break;

          }

          if (aa > aa_orig)

          {

              if (dir == 1)

                  da /= 2;

              dir = -1;

          }

          else

          {

              if (dir == -1)

                  da /= 2;

              dir = 1;

          }

          aa_approx += dir*da*counter;

          aa = aa_approx;

          continue;

      }

      else

          break;

  }

  result->val = aa;

  /* If we went through the maximum number of retries and still didn't

     find the solution, let us know. */

  if (counter == maxcount)

  {

      GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);

  }

  return GSL_SUCCESS;

}

/* Eigenvalue solutions for characteristic values below. */

/*  figi.c converted from EISPACK Fortran FIGI.F.

 *

 *   given a nonsymmetric tridiagonal matrix such that the products

 *    of corresponding pairs of off-diagonal elements are all

 *    non-negative, this subroutine reduces it to a symmetric

 *    tridiagonal matrix with the same eigenvalues.  if, further,

 *    a zero product only occurs when both factors are zero,

 *    the reduced matrix is similar to the original matrix.

 *

 *    on input

 *

 *       n is the order of the matrix.

 *

 *       t contains the input matrix.  its subdiagonal is

 *         stored in the last n-1 positions of the first column,

 *         its diagonal in the n positions of the second column,

 *         and its superdiagonal in the first n-1 positions of

 *         the third column.  t(1,1) and t(n,3) are arbitrary.

 *

 *    on output

 *

 *       t is unaltered.

 *

 *       d contains the diagonal elements of the symmetric matrix.

 *

 *       e contains the subdiagonal elements of the symmetric

 *         matrix in its last n-1 positions.  e(1) is not set.

 *

 *       e2 contains the squares of the corresponding elements of e.

 *         e2 may coincide with e if the squares are not needed.

 *

 *       ierr is set to

 *         zero       for normal return,

 *         n+i        if t(i,1)*t(i-1,3) is negative,

 *         -(3*n+i)   if t(i,1)*t(i-1,3) is zero with one factor

 *                    non-zero.  in this case, the eigenvectors of

 *                    the symmetric matrix are not simply related

 *                    to those of  t  and should not be sought.

 *

 *    questions and comments should be directed to burton s. garbow,

 *    mathematics and computer science div, argonne national laboratory

 *

 *    this version dated august 1983.

 */

static int figi(int nn, double *tt, double *dd, double *ee,

                double *e2)

{

  int ii;

  for (ii=0; ii

  {

      if (ii != 0)

      {

          e2[ii] = tt[3*ii]*tt[3*(ii-1)+2];

          if (e2[ii] < 0.0)

          {

              /* set error -- product of some pair of off-diagonal

                 elements is negative */

              return (nn + ii);

          }

          if (e2[ii] == 0.0 && (tt[3*ii] != 0.0 || tt[3*(ii-1)+2] != 0.0))

          {

              /* set error -- product of some pair of off-diagonal

                 elements is zero with one member non-zero */

              return (-1*(3*nn + ii));

          }

          ee[ii] = sqrt(e2[ii]);

      }

      dd[ii] = tt[3*ii+1];

  }

  return 0;

}

int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])

{

  unsigned int even_order = work->even_order, odd_order = work->odd_order,

      extra_values = work->extra_values, ii, jj;

  int status;

  double *tt = work->tt, *dd = work->dd, *ee = work->ee, *e2 = work->e2,

         *zz = work->zz, *aa = work->aa;

  gsl_matrix_view mat, evec;

  gsl_vector_view eval;

  gsl_eigen_symmv_workspace *wmat = work->wmat;

  if (order_max > work->size || order_max <= order_min || order_min < 0)

    {

      GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);

    }

  /* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal

     form. */

  tt[0] = 0.0;

  tt[1] = 0.0;

  tt[2] = qq;

  for (ii=1; ii

  {

      tt[3*ii] = qq;

      tt[3*ii+1] = 4*ii*ii;

      tt[3*ii+2] = qq;

  }

  tt[3*even_order-3] = qq;

  tt[3*even_order-2] = 4*(even_order - 1)*(even_order - 1);

  tt[3*even_order-1] = 0.0;

  tt[3] *= 2;

  status = figi((signed int)even_order, tt, dd, ee, e2);

  if (status) 

    {

      GSL_ERROR("Internal error in tridiagonal Mathieu matrix", GSL_EFAILED);

    }

  /* Fill the period \pi matrix. */

  for (ii=0; ii

      zz[ii] = 0.0;

  zz[0] = dd[0];

  zz[1] = ee[1];

  for (ii=1; ii

  {

      zz[ii*even_order+ii-1] = ee[ii];

      zz[ii*even_order+ii] = dd[ii];

      zz[ii*even_order+ii+1] = ee[ii+1];

  }

  zz[even_order*(even_order-1)+even_order-2] = ee[even_order-1];

  zz[even_order*even_order-1] = dd[even_order-1];

  /* Compute (and sort) the eigenvalues of the matrix. */

  mat = gsl_matrix_view_array(zz, even_order, even_order);

  eval = gsl_vector_subvector(work->eval, 0, even_order);

  evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);

  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);

  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  for (ii=0; ii

      aa[2*ii] = gsl_vector_get(&eval.vector, ii);

  /* Fill the period 2\pi matrix. */

  for (ii=0; ii

      zz[ii] = 0.0;

  for (ii=0; ii

      for (jj=0; jj

      {

          if (ii == jj)

              zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);

          else if (ii == jj + 1 || ii + 1 == jj)

              zz[ii*odd_order+jj] = qq;

      }

  zz[0] += qq;

  /* Compute (and sort) the eigenvalues of the matrix. */

  mat = gsl_matrix_view_array(zz, odd_order, odd_order);

  eval = gsl_vector_subvector(work->eval, 0, odd_order);

  evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);

  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);

  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  for (ii=0; ii

      aa[2*ii+1] = gsl_vector_get(&eval.vector, ii);

  for (ii = order_min ; ii <= order_max ; ii++)

    {

      result_array[ii - order_min] = aa[ii];

    }

  return GSL_SUCCESS;

}

int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])

{

  unsigned int even_order = work->even_order-1, odd_order = work->odd_order,

      extra_values = work->extra_values, ii, jj;

  double *zz = work->zz, *bb = work->bb;

  gsl_matrix_view mat, evec;

  gsl_vector_view eval;

  gsl_eigen_symmv_workspace *wmat = work->wmat;

  if (order_max > work->size || order_max <= order_min || order_min < 0)

    {

      GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);

    }

  /* Fill the period \pi matrix. */

  for (ii=0; ii

      zz[ii] = 0.0;

  for (ii=0; ii

      for (jj=0; jj

      {

          if (ii == jj)

              zz[ii*even_order+jj] = 4*(ii + 1)*(ii + 1);

          else if (ii == jj + 1 || ii + 1 == jj)

              zz[ii*even_order+jj] = qq;

      }

  /* Compute (and sort) the eigenvalues of the matrix. */

  mat = gsl_matrix_view_array(zz, even_order, even_order);

  eval = gsl_vector_subvector(work->eval, 0, even_order);

  evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);

  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);

  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  bb[0] = 0.0;

  for (ii=0; ii

      bb[2*(ii+1)] = gsl_vector_get(&eval.vector, ii);

  /* Fill the period 2\pi matrix. */

  for (ii=0; ii

      zz[ii] = 0.0;

  for (ii=0; ii

      for (jj=0; jj

      {

          if (ii == jj)

              zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);

          else if (ii == jj + 1 || ii + 1 == jj)

              zz[ii*odd_order+jj] = qq;

      }

  zz[0] -= qq;

  /* Compute (and sort) the eigenvalues of the matrix. */

  mat = gsl_matrix_view_array(zz, odd_order, odd_order);

  eval = gsl_vector_subvector(work->eval, 0, odd_order);

  evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);

  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);

  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  for (ii=0; ii

      bb[2*ii+1] = gsl_vector_get(&eval.vector, ii);  

  for (ii = order_min ; ii <= order_max ; ii++)

    {

      result_array[ii - order_min] = bb[ii];

    }

  return GSL_SUCCESS;

}

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

#include "eval.h"