#include "sph_bessel.h"
#include <limits>

using namespace std;
using namespace sph_bessel;

#define DG 12

Complex sph_bessel::besj(Complex z, int n) {
	if (n == 0) return besj0(z);
	else if (n == 1) return besj1(z);
	//	else if (n < 0) return (n%2) ? besy(z,-n-1) : -besy(z,-n-1);

	if (abs(z) < 0.6) {
		int i, ip;
		Real tvv;
		Complex tc, tcc;

		tc = 1.; for (i=1; i<n+1; i++) tc *= z/Real(2*i+1); i = 1; tcc = tc;
		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
		do {tcc += ( tc *= -Real(0.5)*z*z/Real(i*(2*(n+i)+1)) ); i++;} while (abs(tc) > tvv);

		return tcc;
	}
	else if (abs(z)/Real(n) < 1.) {
		int nn, i, pw;
		Real tv1, tv2, tv3, tx, tx1, tx2, ty1, ty2;
		Complex tc, tc1, tc2, tcc;

		if (n > int(abs(z))) pw = abs(int( (n*log(0.5*_e*abs(z)/n) - 0.5*log(n*abs(z)) - _ln2)/_ln10 )) + DG;
		else pw = DG;

		tv1 = abs(z); tv2 = 2./_e/tv1; tv3 = pw*_ln10 - _ln2;
		tx2 = (tx1 = n) + 2;
		ty1 = tx1*log(tv2*tx1) + 0.5*log(tv1*tx1) - tv3;
		ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
		do {
			tx = tx2 - ty2*(tx2-tx1)/(ty2-ty1); tx1 = tx2; tx2 = tx;
			ty1 = ty2; ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
		} while (fabs(tx2-tx1) > 0.5);
		nn = int(tx2);

		tc1 = 0.; tc2 = exp(-pw*_ln10);
		for (i=nn; i>0; i--) {
			tc = Real(2*i+1)/z*tc2 - tc1; tc1 = tc2; tc2 = tc;
			if (i == n+1) tcc = tc2;
		}
		return (abs(besj1(z)) > abs(besj0(z))) ? tcc*besj1(z)/tc1 : tcc*besj0(z)/tc2;
	}
	else {
		int i, tm;
		Real tvv;
		Complex tc, tp, tq, pp, qq;
		i = -int(0.4343*log(abs(z)/n)) - DG; tvv = (i < -100) ? 1.e-100 : exp(2.3*i);
		tp = pp = 1.; tm = (2*n+1)*(2*n+1); tq = qq = Real(tm-1)*(tc = Real(0.125)/z); tc *= tc; i = 1;
		do {
			pp += ( tp *= -tc*Real((tm - (4*i-1)*(4*i-1))*(tm - (4*i-3)*(4*i-3)))/Real(2*i*(2*i-1)) );
			qq += ( tq *= -tc*Real((tm - (4*i+1)*(4*i+1))*(tm - (4*i-1)*(4*i-1)))/Real(2*i*(2*i+1)) );
			i++;
		} while ((abs(tp) > tvv) && (abs(tq) > tvv));
		return (pp*cos(z-_pi_2*(n+1)) - qq*sin(z-_pi_2*(n+1)))/z;
	}
}

Complex sph_bessel::besjd(Complex z, int n) {
	if (n == 0) return besj0d(z);
	else if (n == 1) return besj1d(z);
	//	else if (n < 0) return (n%2) ? besy(z,-n-1) : -besy(z,-n-1);

	if (abs(z) < 0.6) {
		int i, ip;
		Real tv, tvv;
		Complex tc, tcc;

		tc = 1/3.; for (i=2; i<n+1; i++) tc *= z/Real(2*i+1); i = 1; tcc = Real(n)*tc;
		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
		do {tv = abs( tc *= -Real(0.5)*z*z/Real(i*(2*(n+i)+1)) ); tcc += Real(n+2*i)*tc; i++;} while (tv > tvv);

		return tcc;
	}
	else if (abs(z)/n < 1.) {
		int nn, i, pw;
		Real tv1, tv2, tv3, tx, tx1, tx2, ty1, ty2;
		Complex tc, tc1, tc2, tcc;

		if (n > int(abs(z))) pw = abs(int( (n*log(0.5*_e*abs(z)/n) - 0.5*log(n*abs(z)) - _ln2)/_ln10 )) + DG;
		else pw = DG;

		tv1 = abs(z); tv2 = 2./_e/tv1; tv3 = pw*_ln10 - _ln2;
		tx2 = (tx1 = n) + 2;
		ty1 = tx1*log(tv2*tx1) + 0.5*log(tv1*tx1) - tv3;
		ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
		do {
			tx = tx2 - ty2*(tx2-tx1)/(ty2-ty1); tx1 = tx2; tx2 = tx;
			ty1 = ty2; ty2 = tx2*log(tv2*tx2) + 0.5*log(tv1*tx2) - tv3;
		} while (fabs(tx2-tx1) > 0.5);
		nn = int(tx2);

		tc1 = 0.; tc2 = exp(-pw*_ln10);
		for (i=nn; i>0; i--) {
			tc = Real(2*i+1)/z*tc2 - tc1; tc1 = tc2; tc2 = tc;
			if (i == n+1) tcc = Real(n)/z*tc2 - tc1;
		}
		return (abs(besj1(z)) > abs(besj0(z))) ? tcc*besj1(z)/tc1 : tcc*besj0(z)/tc2;
	}
	else {
		return (Real(n)*besj(z,n-1) - Real(n+1)*besj(z,n+1))/Real(2*n+1);
	}
}

Complex sph_bessel::besy(Complex z, int n) {
	if (n == 0) return besy0(z);
	else if (n == 1) return besy1(z);
	//	else if (n < 0) return (abs(n)%2) ? besj(z,n) : -besj(z,n);

	if (abs(z) < 0.6) {
		int i, ip;
		Real tv, tvv;
		Complex tc, tcc;

		tc = -Real(1)/z; for (i=1; i<n+1; i++) tc *= Real(2*i-1)/z; i = 1; tcc = tc;
		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
		do {tv = abs( tc *= Real(0.5)*z*z/Real(2*(n-i+1)-1)/Real(i) ); tcc += tc; i++;} while (tv > tvv);

		return tcc;
	}
	else if (abs(z)/n < 1.) {
		Complex tc, tc1 = besy0(z), tc2 = besy1(z);
		for (int i=2; i<n+1; i++) {
			tc = Real(2*i-1)*tc2/z - tc1; tc1 = tc2; tc2 = tc;
		}
		return tc2;
	}
	else {
		int i, tm;
		Real tvv;
		Complex tc, tp, tq, pp, qq;
		i = -int(0.4343*log(abs(z))) - DG; tvv = (i < -100) ? 1.e-100 : exp(2.3*i);
		tp = pp = 1.; tm = (2*n+1)*(2*n+1); tq = qq = Real(tm-1)*(tc = Real(0.125)/z); tc *= tc; i = 1;
		do {
			pp += ( tp *= -tc*Real((tm - (4*i-1)*(4*i-1))*(tm - (4*i-3)*(4*i-3)))/Real(2*i*(2*i-1)) );
			qq += ( tq *= -tc*Real((tm - (4*i+1)*(4*i+1))*(tm - (4*i-1)*(4*i-1)))/Real(2*i*(2*i+1)) );
			i++;
		} while ((abs(tp) > tvv) && (abs(tq) > tvv));
		return (pp*sin(z-_pi_2*(n+1)) + qq*cos(z-_pi_2*(n+1)))/z;
	}
}

Complex sph_bessel::besyd(Complex z, int n) {
	if (n == 0) return besy0d(z);
	else if (n == 1) return besy1d(z);
	//	else if (n < 0) return (abs(n)%2) ? besj(z,n) : -besj(z,n);

	if (abs(z) < 0.6) {
		int i, ip;
		Real tvv;
		Complex tc, tcc;
		tc = Real(1)/z/z; for (i=1; i<n+1; i++) tc *= Real(2*i-1)/z; i = 1; tcc = Real(n+1)*tc;
		ip = int(0.4343*log(abs(tc))) - DG; tvv = (ip < -100) ? 1.e-100 : exp(2.3*ip);
		do {tcc -= Real(2*i-n-1)*( tc *= Real(0.5)*z*z/Real(2*(n-i+1)-1)/Real(i) ); i++;} while (abs(tc) > tvv);
		return tcc;
	}
	else if (abs(z)/n < 1.) {
		Complex tc, tc1 = besy0(z), tc2 = besy1(z);
		for (int i=2; i<n+1; i++) {
			tc = Real(2*i-1)*tc2/z - tc1; tc1 = tc2; tc2 = tc;
		}
		return (tc1 - Real(n+1)/z*tc2);
	}
	else return (Real(n)*besy(z,n-1) - Real(n+1)*besy(z,n+1))/Real(2*n+1);
}

Complex sph_bessel::besh1(Complex z, int n) {
	if (n == 0) return besh10(z);
	else if (n == 1) return besh11(z);
	//	return (besj(z,n) + j_*besy(z,n));

	if (z.imag() > -1.e-16*abs(z)) {
		Complex tc, tc1 = besh10(z), tc2 = besh11(z); int i = 1;
		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
		return tc2;
	}
	else return (Real(2)*besj(z,n) - besh2(z,n));
}

Complex sph_bessel::besh2(Complex z, int n) {
	if (n == 0) return besh20(z);
	else if (n == 1) return besh21(z);
	//	return (besj(z,n) - j_*besy(z,n));

	if (z.imag() < 1.e-16*abs(z)) {
		Complex tc, tc1 = besh20(z), tc2 = besh21(z); int i = 1;
		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
		return tc2;
	}
	else return (Real(2)*besj(z,n) + besh1(z,n));
}

Complex sph_bessel::besh1d(Complex z, int n) {
	if (n == 0) return besh10d(z);
	else if (n == 1) return besh11d(z);
	//	return (besjd(z,n) + j_*besyd(z,n));

	if (z.imag() > -1.e-16*abs(z)) {
		Complex tc, tc1 = besh10(z), tc2 = besh11(z); int i = 1;
		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
		return (tc1 - Real(n+1)*tc2/z);
	}
	else return (Real(2)*besjd(z,n) - besh2d(z,n));
}

Complex sph_bessel::besh2d(Complex z, int n) {
	if (n == 0) return besh20d(z);
	else if (n == 1) return besh21d(z);
	//	return (besjd(z,n) - j_*besyd(z,n));

	if (z.imag() < 1.e-16*abs(z)) {
		Complex tc, tc1 = besh20(z), tc2 = besh21(z); int i = 1;
		do {tc = tc2*Real(2*i+1)/z-tc1; tc1 = tc2; tc2 = tc;} while (++i < n);
		return (tc1 - Real(n+1)*tc2/z);
	}
	else return (Real(2)*besjd(z,n) + besh1d(z,n));
}

//###################################################

void get_pitau_m01(Real x, Real *pi1, Real *tau0, Real *tau1, int deg_max) {
	if (fabs(x-1) < numeric_limits<Real>::epsilon()) {
		memset(tau0,0,(deg_max-1)*sizeof(Real));
		for (int n=0; n<deg_max; ++n) {
			pi1[n] = tau1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3));
		}
	}
	else if (fabs(x+1) < numeric_limits<Real>::epsilon()) {
		memset(tau0,0,(deg_max-1)*sizeof(Real));
		for (int n=0; n<deg_max; n+=2) {
			tau1[n] = -(pi1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3)));
		}
		for (int n=1; n<deg_max; n+=2) {
			pi1[n] = -(tau1[n] = 0.5*sqrt(0.5*(n+1)*(n+2)*(2*n+3)));
		}
	}
	else {
		Real pl1, pl2, tmp;
		Real tsin2 = Real(1)-x*x, tsin = sqrt(tsin2);

		pl1 = _1sq2;
		pl2 = sqrt(1.5)*x;
		pi1[0] = 0.5*sqrt(3);
		tau0[0] = -sqrt(1.5)*tsin;
		tau1[0] = 0.5*sqrt(3)*x;

		for (int n=1; n<deg_max; ++n) {
			tmp = sqrt(Real(2*n+3))/Real(n+1)*(x*pl2*sqrt(Real(2*n+1)) - Real(n)/sqrt(Real(2*n-1))*pl1);
			pl1 = pl2; pl2 = tmp;
			tmp = sqrt(Real(n+1)/Real(n+2))/tsin2;
			tau1[n] = sqrt(Real(2*n+3)/Real(2*n+1))*pl1;
			pi1[n] = -tmp*( tau0[n] = x*pl2 - tau1[n] );
			tau0[n] *= Real(n+1)/tsin;
			tau1[n] = tmp*( (Real(1)+(n+1)*tsin2)*pl2 - x*tau1[n] );
		}
	}
}