#ifndef _SPH_BESSEL_H
#define _SPH_BESSEL_H

#include <iostream>
#include <complex>

typedef double Real;
typedef std::complex<Real> Complex;

const Complex im1(Real(0),Real(1));
const Real _pi = 3.14159265358979323846;
const Real _pi_2 = 1.57079632679489661923;
const Real _e = 2.71828182845904523536;
const Real _ln10 = 2.30258509299404568402;
const Real _ln2 = 0.693147180559945309417;
const Real _1sq2 = 0.707106781186547524401;

inline Complex get_pow_i1(int p) {
	int i = p%4;
	if (i == 0) return Real(1);
	else if (i == 1) return im1;
	else if (i == 2) return Real(-1);
	else return -im1;
}

namespace sph_bessel {
	inline Complex besj0(Complex z) {if (abs(z) < 1.e-8) return 1.; else return sin(z)/z;};
	inline Complex besj1(Complex z) {if (abs(z) < 1.e-15) return z/Real(3); else return (sin(z)-z*cos(z))/z/z;}
	inline Complex besj0d(Complex z) {if (abs(z) < 1.e-15) return -z/Real(3); else return (z*cos(z)-sin(z))/z/z;}
	inline Complex besj1d(Complex z) {
		if (abs(z) < 1.e-8) return Real(1)/Real(3);
		else return ((z*z - Real(2))*sin(z) + Real(2)*z*cos(z))/(z*z*z);
	}
	inline Complex besy0(Complex z) {return -cos(z)/z;}
	inline Complex besy1(Complex z) {return (-cos(z) - z*sin(z))/z/z;}
	inline Complex besy0d(Complex z) {return (cos(z) + z*sin(z))/z/z;;}
	inline Complex besy1d(Complex z) {return -( (z*z - Real(2))*cos(z) - Real(2)*z*sin(z) )/z/z/z;}
	inline Complex besh10(Complex z) {return -im1*exp(im1*z)/z;}
	inline Complex besh11(Complex z) {return (-z-im1)*exp(im1*z)/z/z;}
	inline Complex besh10d(Complex z) {return (z+im1)*exp(im1*z)/z/z;}
	inline Complex besh11d(Complex z) {return (-im1 + Real(2)*(z+im1)/z/z )*exp(im1*z)/z;}
	inline Complex besh20(Complex z) {return im1*exp(-im1*z)/z;}
	inline Complex besh21(Complex z) {return (-z+im1)*exp(-im1*z)/z/z;}
	inline Complex besh20d(Complex z) {return (z-im1)*exp(-im1*z)/z/z;}
	inline Complex besh21d(Complex z) {return (im1 + Real(2)*(z-im1)/z/z)*exp(-im1*z)/z;}

	Complex besj(Complex, int);
	Complex besjd(Complex, int);
	Complex besy(Complex, int);
	Complex besyd(Complex, int);
	Complex besh1(Complex, int);
	Complex besh2(Complex, int);
	Complex besh1d(Complex, int);
	Complex besh2d(Complex, int);

	inline Complex bes_dzj(Complex z, int n) {return besj(z,n)+z*besjd(z,n);};
	inline Complex bes_dzy(Complex z, int n) {return besy(z,n)+z*besyd(z,n);};
	inline Complex bes_dzh1(Complex z, int n) {return besh1(z,n)+z*besh1d(z,n);};
	inline Complex bes_dzh2(Complex z, int n) {return besh2(z,n)+z*besh2d(z,n);};
}

void get_pitau_m01(Real x, Real *pi1, Real *tau0, Real *tau1, int deg_max);

#endif // _SPH_BESSEL_H