import java.io.*;

import javax.swing.*;

class problem4bessel {

	public static double[] thomas(double a[][], double r[]) {
		//
		int n = a.length;
		double f[] = new double[n];
		double e[] = new double[n];
		double g[] = new double[n];
		double x[] = new double[n];
		for (int i = 0; i < n; i++) {
			f[i] = a[i][i];
		}
		for (int i = 0; i < (n - 1); i++) {
			g[i] = a[i][i + 1];
		}
		for (int i = 0; i < (n - 1); i++) {
			e[i + 1] = a[i + 1][i];
		}
		for (int k = 1; k < n; k++) {
			e[k] = e[k] / f[k - 1];
			f[k] = f[k] - e[k] * g[k - 1];
		}
		for (int k = 1; k < n; k++) {
			r[k] = r[k] - e[k] * r[k - 1];
		}
		x[n - 1] = r[n - 1] / f[n - 1];
		for (int k = (n - 2); k >= 0; k--) {
			x[k] = (r[k] - g[k] * x[k + 1]) / f[k];
		}
		return x;
	}

	public static double[] thomas(double f[], double e[], double g[],
			double r[]) {
		int n = f.length;
		double x[] = new double[n];
		for (int k = 1; k < n; k++) {
			e[k] = e[k] / f[k - 1];
			f[k] = f[k] - e[k] * g[k - 1];
		}
		for (int k = 1; k < n; k++) {
			r[k] = r[k] - e[k] * r[k - 1];
		}
		x[n - 1] = r[n - 1] / f[n - 1];
		for (int k = (n - 2); k >= 0; k--) {
			x[k] = (r[k] - g[k] * x[k + 1]) / f[k];
		}
		return x;
	}

	public static double[][] cubic_spline(double xi[], double yi[], double c0,
			double cn) {
		int n = xi.length;
		double h[] = new double[n];
		double w[] = new double[n];
		double f[] = new double[n];
		double e[] = new double[n];
		double g[] = new double[n];
		double d[] = new double[n];
		double x[] = new double[n];
		double S[][] = new double[4][n];
		int k;
		for (k = 0; k < (n - 1); k++) {
			h[k] = xi[k + 1] - xi[k];
			w[k] = (yi[k + 1] - yi[k]) / h[k];
		}
		d[0] = c0;
		d[n - 1] = cn;
		for (k = 1; k < (n - 1); k++) {
			d[k] = 6.0 * (w[k] - w[k - 1]);
		}
		f[0] = 1.0;
		f[n - 1] = 1.0;
		g[0] = 0.0;
		g[n - 1] = 0.0;
		e[0] = 0.0;
		e[n - 1] = 0.0;
		for (k = 1; k < (n - 1); k++) {
			f[k] = 2.0 * (h[k] + h[k - 1]);
			e[k] = h[k - 1];
			g[k] = h[k];
		}
		S[2] = thomas(f, e, g, d);
		S[3] = xi;
		for (k = 0; k < (n - 1); k++) {
			S[0][k] = (6. * yi[k + 1] - h[k] * h[k] * S[2][k + 1])
					/ (6.0 * h[k]);
			S[1][k] = (6. * yi[k] - h[k] * h[k] * S[2][k]) / (6.0 * h[k]);
		}
		return S;
	}

	public static double funcSpline(double S[][], double x) {
		int n = S[0].length;
		double xx1 = 0;
		double xx2 = 0;
		double y = 0;
		double hk = 0;
		for (int k = 0; k < (n - 1); k++) {
			if (S[3][k] <= x && x <= S[3][k + 1]) {
				hk = (S[3][k + 1] - S[3][k]);
				xx1 = (x - S[3][k]);
				xx2 = (S[3][k + 1] - x);
				y = S[0][k]
						* xx1
						+ S[1][k]
						* xx2
						+ (xx1 * xx1 * xx1 * S[2][k + 1] + xx2 * xx2 * xx2
								* S[2][k]) / (6.0 * hk);
				break;
			}
		}
		if (y == 0 && S[3][n - 2] <= x) {
			int k = n - 2;
			hk = (S[3][k + 1] - S[3][k]);
			xx1 = (x - S[3][k]);
			xx2 = (S[3][k + 1] - x);
			y = S[0][k]
					* xx1
					+ S[1][k]
					* xx2
					+ (xx1 * xx1 * xx1 * S[2][k + 1] + xx2 * xx2 * xx2
							* S[2][k]) / (6.0 * hk);
		}
		return y;
	}

	public static double[][] funcSpline(double xi[], double yi[],
			int aradegersayisi) {
		// aradegersayisi: x--o--o--x--o--o--x zincirinde x deneysel noktalar
		// ise
		// ara de?er say?s? 2 dir
		int n = xi.length;
		int nn = (n - 1) * (aradegersayisi + 1) + 1;
		double z[][] = new double[2][nn];
		double S[][] = cubic_spline(xi, yi, 0, 0);
		double dx = 0;
		int k = 0;
		int i;
		for (i = 0; i < (n - 1); i++) {
			z[0][k] = xi[i];
			z[1][k] = funcSpline(S, z[0][k]);
			k++;
			for (int j = 0; j < aradegersayisi; j++) {
				dx = (xi[i + 1] - xi[i]) / ((double) aradegersayisi + 1.0);
				z[0][k] = z[0][k - 1] + dx;
				z[1][k] = funcSpline(S, z[0][k]);
				k++;
			}
		}
		z[0][k] = xi[i];
		z[1][k] = funcSpline(S, z[0][k]);
		return z;
	}

	public static void main(String args[]) {
		try {
			besselCiz();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	private static void besselCiz() throws IOException {
		double x[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
		double y[];
		double d1[][] = new double[5][];
		double d2[][] = new double[5][];
		
		// v = 0
		y = new double[11];
		for (int i = 0; i <=10; i++) {
			y[i] = besselJ(0, i);
			System.out.println(y[i]);
		}
		int n = x.length;
		double t[] = new double[n];
		for (int i = 0; i < n; i++)
			t[i] = i;
		double z1[][] = funcSpline(t, x, 20);
		double z2[][] = funcSpline(t, y, 20);
		d1[0] = z1[1];
		d2[0] = z2[1];		
		
		// v = 0.5
		y = new double[11];
		for (int i = 0; i <=10; i++) {
			y[i] = besselJ(0.5, i);
			System.out.println(y[i]);
		}
		n = x.length;
		double t1[] = new double[n];
		for (int i = 0; i < n; i++)
			t1[i] = i;
		double z11[][] = funcSpline(t1, x, 20);
		double z21[][] = funcSpline(t1, y, 20);
		d1[1] = z11[1];
		d2[1] = z21[1];
		//
		
		// v = 1
		y = new double[11];
		for (int i = 0; i <=10; i++) {
			y[i] = besselJ(1, i);
			System.out.println(y[i]);
		}
		n = x.length;
		double t2[] = new double[n];
		for (int i = 0; i < n; i++)
			t2[i] = i;
		double z12[][] = funcSpline(t2, x, 20);
		double z22[][] = funcSpline(t2, y, 20);
		d1[2] = z12[1];
		d2[2] = z22[1];
		//
		// v = 1.5
		y = new double[11];
		for (int i = 0; i <=10; i++) {
			y[i] = besselJ(1.5, i);
			System.out.println(y[i]);
		}
		n = x.length;
		double t3[] = new double[n];
		for (int i = 0; i < n; i++)
			t3[i] = i;
		double z13[][] = funcSpline(t3, x, 20);
		double z23[][] = funcSpline(t3, y, 20);
		d1[3] = z13[1];
		d2[3] = z23[1];
		//
		
		// v = 2.0
		y = new double[11];
		for (int i = 0; i <=10; i++) {
			y[i] = besselJ(2.0, i);
			System.out.println(y[i]);
		}
		n = x.length;
		double t4[] = new double[n];
		for (int i = 0; i < n; i++)
			t4[i] = i;
		double z14[][] = funcSpline(t4, x, 20);
		double z24[][] = funcSpline(t4, y, 20);
		d1[4] = z14[1];
		d2[4] = z24[1];
		//
		// System.out.println(" z1\n"+Matrix.toString(z1));
		// System.out.println(" z2\n"+Matrix.toString(z1));
		Plot pp = new Plot(d1, d2);
		pp.setPlabel("kübik serit yöntemi ile parametrik egri uydurma");
		pp.setXlabel("x");
		pp.setYlabel("y");
		pp.setPlotType(0, 0);
		double aa[][] = new double[2][z1[0].length];
		for (int i = 0; i < z1[1].length; i++) {
			aa[0][i] = z1[1][i];
			aa[1][i] = z2[1][i];
		}
		Text.print(Text.T(aa));
		pp.plot();
	}
	
	public static double besselJ(double v,double x) {
		double sonuc = 0;
		double toplam = 0;
		for (int k=0; k<100; k++) {
			double pay = Math.pow((-(x*x) / 4 ), k);
			double b = sayisal.gammln(v + k + 1);
			double a = Math.pow(Math.E, b);
			double payda = fact(k) * a;		
			double bolum = pay / payda;			
			toplam += bolum;
		}
		sonuc = Math.pow((x/2), v) * toplam;
		return sonuc;
	}
	
	public static double fact(int n) {
		double a[]=new double[33];
		int ntop=4;
		a[0]=1.0;
		a[1]=1.0;
		a[2]=2.0;
		a[3]=6.0;
		int j;
		if (n < 0) System.out.println("Negative factorial in routine FACTLN");
		if (n <= 1) return 1.0;
		if (n > 32) return sayisal.gammln(n+1.0);
		while (ntop<=n) {
			j=ntop-1;
			a[ntop]=a[j]*ntop;
			ntop++;
		}
		return a[n];
	}
}
