///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Bessel.java                                        //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of using the Bessel functions generated. 

import java.lang.*;
import java.util.Random;
public class Bessel {
  public static void main(String argv[]) {
    int n = 10, nb = 5;
    Random r = new Random();
    double x = r.nextDouble();

    double[][] f = b(x, n, nb);
    for (int i=0; i<=n; ++i) {
      System.out.println("J_" + i + " is: " + f[0][i]);
      System.out.println("Y_" + i + " is: " + f[1][i]);
    }
  }

// Method to create Bessel functions for a given x, index n,
// and index buffer nb.

  public static double[][] b(double x, int n, int nb) {
    int nmax = n+nb;
    double g = 0.5772156649;
    double y[][] = new double[2][n+1];
    double z[] = new double[nmax+1];

 // Generate the Bessel function of 1st kind J_n(x)
    z[nmax-1] = 1;
    double s = 0; 
    for (int i=nmax-1; i>0; --i) {
      z[i-1] = 2*i*z[i]/x-z[i+1];
      if (i%2 == 0) s += 2*z[i];
    }
    s += z[0];
    for (int i=0; i<=n; ++i) y[0][i] = z[i]/s;

 // Generate the Bessel function of 2nd kind Y_n(x)
    double t = 0;
    int sign = -1;
    for (int i=1; i<nmax/2; ++i) {
      t += sign*z[2*i]/i;
      sign *= -1;
    }
    t *= -4/(Math.PI*s);
    y[1][0] = 2*(Math.log(x/2)+g)*y[0][0]/Math.PI+t;
    y[1][1] = (y[0][1]*y[1][0]-2/(Math.PI*x))/y[0][0];
    for (int i=1; i<n; ++i)
      y[1][i+1] = 2*i*y[1][i]/x-y[1][i-1];
  
    return y;
  }
}