Funciones de Bessel

public abstract class Raiz {
	private final double CERO=1e-10;
	private final double ERROR=0.0001;
	private final int MAXITER=100;
	protected final int MAXRAICES=20;
	protected double raices[]=new double[MAXRAICES];
	protected int iRaiz=0;

protected double puntoMedio(double a, double b)throws RaizExcepcion{
	double m, ym;
	int iter=0;
	do{
		m=(a+b)/2;
		ym=f(m);
		if(Math.abs(ym)<CERO) break;
		if(Math.abs((a-b)/m)<ERROR) break;
		if((f(a)*ym)<0) b=m;
		else a=m;
		iter++;
	}while(iter<MAXITER);
	if(iter==MAXITER){
		throw new RaizExcepcion("No se ha encontrado una raíz");
	}
	return m;
}

private void explorar(double xIni, double xFin, double dx){
	double y1, y2;
	iRaiz=0;
	y1=f(xIni);
	for(double x=xIni; x<xFin; x+=dx){
		y2=f(x+dx);
//Uno de los extremos del intervalo es raíz
		if(Math.abs(y1)<CERO && iRaiz<MAXRAICES){
			raices[iRaiz++]=x;
			y1=y2;
			continue;
		}
//no hay raíz en este intervalo
		if(y1*y2>=0.0){
			y1=y2;
			continue;
		}
//hay una raíz en este intervalo
		if(iRaiz<MAXRAICES){
			try{
				raices[iRaiz]=puntoMedio(x, x+dx);
				iRaiz++;
			}catch(RaizExcepcion ex){
				System.out.println(ex.getMessage());
			}
		}
		y1=y2;
	}
}

public double[] hallarRaices(double ini, double fin, double paso){
	explorar(ini, fin, paso);
	double solucion[]=new double[iRaiz];
	for(int i=0; i<iRaiz; i++){
		solucion[i]=(double)Math.round(raices[i]*1000)/1000;
	}
	return solucion;
}

abstract protected double f(double x);
}

class RaizExcepcion extends Exception {
	public RaizExcepcion(String s) {
	super(s);
}
}
public class Maximos extends Raiz {
protected double f(double x){
	return Bessel.bessj(2, x);
}

}
public class Minimos extends Raiz {
protected double f(double x){
	return Bessel.bessj1(x);
}


public class Difraccion {

public static void main(String[] args) {
	double[] raices=new Maximos().hallarRaices(0.0, 30.0, 1.0);
	System.out.print("Máximos");
	for(int i=0; i<raices.length; i++){
		System.out.print(" "+raices[i]);
	}
	System.out.println("");
	raices=new Minimos().hallarRaices(0.0, 30.0, 1.0);
	System.out.print("Mínimos");
	for(int i=0; i<raices.length; i++){
		System.out.print(" "+raices[i]);
	}
}
}
public class Bessel {
public static double bessj0(double x){
//Returns the Bessel function J0(x) for any real x.
	double ax,z;
	double xx,y,ans,ans1,ans2; //Accumulate polynomials in double precision.
	if ((ax=Math.abs(x)) < 8.0) { //Direct rational function t.
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7+y*(
-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+
y*(9494680.718+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans=ans1/ans2;
	} else { //Fitting function (6.5.9).
		z=8.0/ax;
		y=z*z;
		xx=ax-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4 +
y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3 +y*(-0.6911147651e-5+
		y*(0.7621095161e-6-y*0.934935152e-7)));
		ans=Math.sqrt(0.636619772/ax)*(Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
	}
	return ans;
}
}
public static double bessj1(double x){
//Returns the Bessel function J1(x) for any real x.
	double ax,z;
	double xx,y,ans,ans1,ans2; //Accumulate polynomials in double precision.
	if ((ax=Math.abs(x)) < 8.0) { //Direct rational approximation.
		y=x*x;
		ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+
			y*(15704.48260+y*(-30.16036606))))));
		ans2=144725228442.0+y*(2300535178.0
+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
		ans=ans1/ans2;
	} else { //Fitting function (6.5.9).
		z=8.0/ax;
		y=z*z;
		xx=ax-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+
		y*(-0.88228987e-6+y*0.105787412e-6)));
		ans=Math.sqrt(0.636619772/ax)*(Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
		if (x < 0.0) ans = -ans;
	}
	return ans;
}
  public static double bessj(int n, double x){
//Returns the Bessel function Jn(x) for any real x and n > 2.
    final double ACC=40.0; // Make larger to increase accuracy.
    final double BIGNO=1.0e10;
    final double BIGNI=1.0e-10;
    //int j,jsum,m;
    boolean jsum;
    int j, m;
    double ax,bj,bjm,bjp,sum,tox,ans;
    if(n==0) return bessj0(x);
    if(n==1) return bessj1(x);
    if(n==-1) return -bessj1(x);
//    if (n < 2) nrerror("Index n less than 2 in bessj");
    ax=Math.abs(x);
    if (ax == 0.0)
        return 0.0;
    else if (ax > (double) n) { //Upwards recurrence from J0 and J1.
        tox=2.0/ax;
        bjm=bessj0(ax);
        bj=bessj1(ax);
        for (j=1;j<n;j++) {
            bjp=j*tox*bj-bjm;
            bjm=bj;
            bj=bjp;
        }
        ans=bj;
    } else { //Downwards recurrence from an even m here computed.
        tox=2.0/ax;
        m=2*((n+(int) Math.sqrt(ACC*n))/2);
        jsum=false; //jsum will alternate between 0 and 1; 
//when it is 1, we accumulate in sum the even terms in (5.5.16).
        bjp=ans=sum=0.0;
        bj=1.0;
        for (j=m;j>0;j--) {// The downward recurrence.
            bjm=j*tox*bj-bjp;
            bjp=bj;
            bj=bjm;
            if (Math.abs(bj) > BIGNO) { //Renormalize to prevent overflows.
                bj *= BIGNI;
                bjp *= BIGNI;
                ans *= BIGNI;
                sum *= BIGNI;
            }
            if (jsum) sum += bj;  // Accumulate the sum.
            jsum=!jsum;         //Change 0 to 1 or vice versa.
            if (j == n) ans=bjp; //Save the unnormalized answer.
        }
        sum=2.0*sum-bj;         //Compute (5.5.16)
        ans /= sum;             //and use it to normalize the answer.
    }
    boolean bRes=((n&1)==0)? false:true;
    return (x < 0.0) && bRes ? -ans : ans;
 }
}