#include <iostream>
#include <iomanip>
using namespace std;

#include "MpIeee.hh"
#include "CMpIeee.hh"
#include "ArithmosIO.hh"

/* specfunc/bessel_Jn.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * 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:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_pow_int.h>
#include "bessel.h"
#include "bessel_amp_phase.h"
#include "bessel_olver.h"
#include <gsl/gsl_sf_bessel.h>



/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/


int  gsl_sf_bessel_Jn_e(int  n, MpIeee x, gsl_sf_result * result)
{
  int  sign=  1;

  if(n < 0) {
    /* reduce to case n >= 0 */
    n = -n;
    if(GSL_IS_ODD(n)) sign = -sign;
  }  

  if(x < MpIeee( "0.0" )) {
    /* reduce to case x >= 0. */
    x = -x;
    if(GSL_IS_ODD(n)) sign = -sign;
  }

  /* CHECK_POINTER(result) */

  if(n == 0) {
    gsl_sf_result b0;
    int  stat_J0=  gsl_sf_bessel_J0_e(x, &b0);
    result->val = sign * b0.val;
    result->err = b0.err;
    return stat_J0;
  }
  else if(n == 1) {
    gsl_sf_result b1;
    int  stat_J1=  gsl_sf_bessel_J1_e(x, &b1);
    result->val = sign * b1.val;
    result->err = b1.err;
    return stat_J1;
  }
  else {
    if(x == MpIeee( "0.0" )) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else if(x*x < MpIeee( "10.0" )*(n+MpIeee( "1.0" ))*GSL_ROOT5_DBL_EPSILON) {
      gsl_sf_result b;
      int  status=  gsl_sf_bessel_IJ_taylor_e((double)n, x, -1, 50, GSL_DBL_EPSILON, &b);
      result->val  = sign * b.val;
      result->err  = b.err;
      result->err += GSL_DBL_EPSILON * fabs(result->val);
      return status;
    }
    else if(GSL_ROOT4_DBL_EPSILON * x > (n*n+1.0)) {
      int  status=  gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else if(n > 50) {
      int  status=  gsl_sf_bessel_Jnu_asymp_Olver_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else if(x > 1000.0)
    {
      /* We need this to avoid feeding large x to CF1; note that
       * due to the above check, we know that n <= 50.
       */
      int  status=  gsl_sf_bessel_Jnu_asympx_e((double)n, x, result);
      result->val *= sign;
      return status;      
    }
    else {
      MpIeee ans;
      MpIeee err;
      MpIeee ratio;
      MpIeee sgn;
      int  stat_b;
      int  stat_CF1=  gsl_sf_bessel_J_CF1((double)n, x, &ratio, &sgn);

      /* backward recurrence */
      MpIeee Jkp1=  GSL_SQRT_DBL_MIN * ratio;
      MpIeee Jk=  GSL_SQRT_DBL_MIN;
      MpIeee Jkm1;
      int  k;

      for(k=n; k>0; k--) {
        Jkm1 = MpIeee( "2.0" )*k/x * Jk - Jkp1;
        Jkp1 = Jk;
        Jk   = Jkm1;
      }

      if(fabs(Jkp1) > fabs(Jk)) {
        gsl_sf_result b1;
        stat_b = gsl_sf_bessel_J1_e(x, &b1);
        ans = b1.val/Jkp1 * GSL_SQRT_DBL_MIN;
        err = b1.err/Jkp1 * GSL_SQRT_DBL_MIN;
      }
      else {
        gsl_sf_result b0;
        stat_b = gsl_sf_bessel_J0_e(x, &b0);
        ans = b0.val/Jk * GSL_SQRT_DBL_MIN;
        err = b0.err/Jk * GSL_SQRT_DBL_MIN;
      }

      result->val = sign * ans;
      result->err = fabs(err);
      return GSL_ERROR_SELECT_2(stat_CF1, stat_b);
    }
  }
}


int
 gsl_sf_bessel_Jn_array(int  nmin, int  nmax, MpIeee x, MpIeee * result_array)
{
  /* CHECK_POINTER(result_array) */

  if(nmin < 0 || nmax < nmin) {
    int  n;
    for(n=nmax; n>=nmin; n--) {
      result_array[n-nmin] = MpIeee( "0.0" );
    }
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(x == MpIeee( "0.0" )) {
    int  n;
    for(n=nmax; n>=nmin; n--) {
      result_array[n-nmin] = MpIeee( "0.0" );
    }
    if(nmin == 0) result_array[0] = MpIeee( "1.0" );
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result r_Jnp1;
    gsl_sf_result r_Jn;
    int  stat_np1=  gsl_sf_bessel_Jn_e(nmax+1, x, &r_Jnp1);
    int  stat_n=  gsl_sf_bessel_Jn_e(nmax,   x, &r_Jn);
    int  stat=  GSL_ERROR_SELECT_2(stat_np1, stat_n);

    MpIeee Jnp1=  r_Jnp1.val;
    MpIeee Jn=  r_Jn.val;
    MpIeee Jnm1;
    int  n;

    if(stat == GSL_SUCCESS) {
      for(n=nmax; n>=nmin; n--) {
        result_array[n-nmin] = Jn;
        Jnm1 = -Jnp1 + MpIeee( "2.0" )*n/x * Jn;
        Jnp1 = Jn;
        Jn   = Jnm1;
      }
    }
    else {
      for(n=nmax; n>=nmin; n--) {
        result_array[n-nmin] = MpIeee( "0.0" );
      }
    }

    return stat;
  }
}

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

#include "eval.h"

MpIeee gsl_sf_bessel_Jn(const int n, const MpIeee x)
{
  EVAL_RESULT(gsl_sf_bessel_Jn_e(n, x, &result));
}
