////-------------------------------------------------------------------//

//  Syntax:  besselj ( ALPHA, X )

//  Description:

//  Besselj is the bessel function of the first kind of order ALPHA.
//  Either ALPHA or X may a vector in which case a vector result of
//  the same dimension is returned.  If both ALPHA and X are vectors
//  a matrix result of dimension length(ALPHA) x length(X) is returned.
//  If either ALPHA or X is a scalar, the other argument may be a matrix
//  and a matrix result of the same dimension is returned.

//-------------------------------------------------------------------//

// See Numerical Recipes in C (second edition)
// Currently only integer order supported.

static (besseljn, besselj0, besselj1);

besselj = function(alph,x) 
{
  global (pi)

  // Checking the class will automatically result in an error if
  // the argument doesn't exist.
  if(class(alph) != "num" || class(x) != "num") {
    error("besselj: argument is non-numeric.");
  }

  if(min(size(alph)) > 1 && min(size(x)) > 1) {
    // Both arguments are matrices.
    error("besselj: only one argument may be a matrix.");
  }

  if(length(alph) == 1) {
    // x is a matrix (or scalar), alph a scalar
    if(alph == 0 || mod(alph,int(alph)) == 0) {
      return besseljn(alph,x);
    else
      error("besselj: fractional orders not yet supported.");
      // return besseljr(alph,x);
    }
  else if(length(x) == 1) {
    // alph is a matrix (or vector), x a scalar
    y = zeros(size(alph));

    // Use integer order routines for integer alph
    inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
    frac_ind = complement(inte_ind,1:alph.n);

    if(inte_ind.n) {
      for(index in inte_ind) {
        y[index] = besseljn(alph[index],x);
      }
    }
    if(frac_ind.n) {
      error("besselj: fractional orders not yet supported.");
      // for(index in frac_ind) {
        // y[index] = besseljr(alph[index],x);
      // }
    }
    return y;
  else
    // alph and x are both vectors
    y = zeros(length(alph),length(x));

    // Use integer order routines for integer alph
    inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
    frac_ind = complement(inte_ind,1:alph.n);

    if(inte_ind.n) {
      for(index in inte_ind) {
        y[index;] = besseljn(alph[index],x);
      }
    }
    if(frac_ind.n) {
      error("besselj: fractional orders not yet supported.");
      // for(index in frac_ind) {
        // y[index;] = besseljr(alph[index],x);
      // }
    }
    return y;
  }}
};

// In these functions x can be a vector (or matrix), but n
// must be a scalar.

// No argument checking is performed on static functions.

besseljn = function ( n, x ) {

  local(abs_x,y,index,index_1,index_2,index_3,arg,p_1,p_2,p_3, ...
        two_over_x,m,y_2,Norm,even,limit,lim_ind);

  if(n == 0) {
    return besselj0(x);
  else if (n == 1) {
    return besselj1(x);
  else
    // integer order greater than two
    abs_x = abs(x);
    y = zeros(size(x));
    index_1 = find(abs_x > n);
    index_2 = complement(index_1,1:x.n);
    index_3 = find(abs_x == 0);

    index_1 = complement(intersection(index_1,index_3),index_1);
    index_2 = complement(intersection(index_2,index_3),index_2);

    if(index_1.n) {
      arg = x[index_1];

      p_1 = besselj0(arg);
      p_2 = besselj1(arg);
      two_over_x = 2./arg;
      for(index in 1:(n-1)) {
        p_3 = index * (two_over_x .* p_2) - p_1;
        p_1 = p_2;
        p_2 = p_3;
      }
      if(mod(n,2) == 1) {
        index = find(arg < 0.);
        if(index.n) {
          p_3[index] = -p_3[index];
        }
      }
      y[index_1] = p_3;
    }

    if(index_2.n) {
      limit=1.e20;
      arg = x[index_2];
      two_over_x = 2./arg;
      // Downward recurrence from arbitrary starting values.
      // Find a starting upper even index (this give 16 digits
      // of accuracy).
      m = 2*int((n + 16*sqrt(n))/2);
      p_1 = zeros(size(arg));
      p_2 = ones(size(arg));
      Norm = zeros(size(arg));
      y_2 = zeros(size(arg));
      even = 0;
      for(index in m:1:-1) {
        p_3 = index * (two_over_x .* p_2) - p_1;
        p_1 = p_2;
        p_2 = p_3;
        if(even) {
          Norm = Norm + p_2;
        }
        lim_ind = find(abs(p_2) > limit);
        if(!isempty(lim_ind)) {
          p_2[lim_ind] = p_2[lim_ind] / limit;
          p_1[lim_ind] = p_1[lim_ind] / limit;
          Norm[lim_ind] = Norm[lim_ind] / limit;
          if(index < n) {
            y_2[lim_ind] = y_2[lim_ind] / limit;
          }
        }
        even = !even;

        if(index == n) {
          // Save the unnormalized answer, continue looping to get
          // normalization.
          y_2 = p_1;
        }
      }
      // Only need 1*j0(x) so subtract off p_2 after doubling.
      Norm = 2*Norm - p_2;
      y_2 = y_2 ./ Norm;
      
      if(mod(n,2) == 1) {
        index = find(arg < 0.);
        if(index.n) {
          y_2[index] = -y_2[index];
        }
      }
      y[index_2] = y_2;
    }

    if(index_3.n) {
      y[index_3] = zeros(size(index_3));
    }

    return y;
  }}
};

besselj0 = function(x) {
  local(abs_x,index_1,index_2,y,c_1,c_2,arg_1,arg_2,p_1,p_2);

  abs_x = abs(x);

  index_1 = find(abs_x < 8);
  index_2 = complement(index_1,1:x.n);

  // First evaluate using rational approx. for small arguments.

  y = zeros(size(x));

  if(index_1.n) {
    c_1 = [ 57568490574, -13362590354, 651619640.7, ...
           -11214424.18, 77392.33017, -184.9052456];
    c_2 = [ 57568490411, 1029532985, 9494680.718, ...
            59272.64853, 267.8532712];

    arg_1 = x[index_1] .* x[index_1];

    p_1 = c_1[1] + arg_1 .* (c_1[2] + arg_1 .* (c_1[3] + arg_1 .* ...
          (c_1[4] + arg_1 .* (c_1[5] + arg_1 .* c_1[6]))));
    p_2 = c_2[1] + arg_1 .* (c_2[2] + arg_1 .* (c_2[3] + arg_1 .* ...
          (c_2[4] + arg_1 .* (c_2[5] + arg_1))));
    y[index_1] = p_1./p_2;
  }

  // Now evaluate for large arguments using approximating form.

  if(index_2.n) {
    c_1 = [ -0.1098628627e-2, 0.2734510407e-4, ...
            -0.2073370639e-5, 0.2093887211e-6 ];
    c_2 = [ -0.1562499995e-1, 0.1430488765e-3, ...
            -0.6911147651e-5, 0.7621095161e-6, 0.934935152e-7 ];

    arg_1 = 8./x[index_2];
    arg_2 = arg_1 .* arg_1;

    p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
          arg_2 .* c_1[4])));
    p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
          arg_2 .* c_2[5])));

    c_1 = abs_x[index_2];
    c_2 = c_1 - pi/4;
    y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*cos(c_2) - arg_1.*p_2.*sin(c_2));
  }

  return y;
};


besselj1 = function(x) {
  local(abs_x,index_1,index_2,y,c_1,c_2,arg_1,arg_2,p_1,p_2);

  abs_x = abs(x);

  index_1 = find(abs_x < 8);
  index_2 = complement(index_1,1:x.n);

  // First evaluate using rational approx. for small arguments.

  y = zeros(size(x));

  if(index_1.n) {
    c_1 = [ 72362614232, -7895059235, 242396853.1, ...
            -2972611.439, 15704.48260, -30.16036606 ]; 
           
    c_2 = [ 144725228442, 2300535178, 18583304.74, ...
            99447.43394, 376.9991397];

    
    arg_1 = x[index_1];
    arg_2 = arg_1 .* arg_1;

    p_1 = arg_1.*(c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
          (c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6])))));
    p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
          (c_2[4] + arg_2 .* (c_2[5] + arg_2))));
    y[index_1] = p_1./p_2;
  }

  // Now evaluate for large arguments using approximating form.

  if(index_2.n) {
    c_1 = [ 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, ...
            -0.240337019e-6 ];
    c_2 = [ 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, ...
            -0.88228987e-6, 0.105787412e-6 ];

    arg_1 = 8./x[index_2];
    arg_2 = arg_1 .* arg_1;

    p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
          arg_2 .* c_1[4])));
    p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
          arg_2 .* c_2[5])));

    c_1 = abs_x[index_2];
    c_2 = c_1 - 3*pi/4;
    arg_2 = sqrt(2./(pi*c_1)).*(p_1.*cos(c_2) - arg_1.*p_2.*sin(c_2));
    index_1 = find(arg_1 < 0);
    if(!isempty(index_1)) {
      arg_2[index_1] = -arg_2[index_1];
    }

    y[index_2] = arg_2;
  }

  return y;
};