#include "sphere.h"

#include "rts/math/complex.h"
#include <complex>
#include <stdlib.h>

using namespace rts;
using namespace std;

int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
    complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);

int cbessjyva_sph(int v,complex<double> z,double &vm,complex<double>*cjv,
    complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);

void sphere::calcCoeff(ptype lambda, rtsComplex<ptype> ri)
{
    /*  These calculations are done at high-precision on the CPU
        since they are only required once for each wavelength.

        Input:

        lambda  =   wavelength of the incident field
        n       =   complex refractive index of the sphere
    */

    //clear the previous coefficients
    A.clear();
    B.clear();

	//convert to an std complex value
	complex<double> nc(ri.real(), ri.imag());
	n = ri;

    //compute the magnitude of the k vector
    double k = 2 * PI / lambda;
    complex<double> kna = nc * k * (double)a;

    //compute the arguments k*a and k*n*a
    complex<double> ka(k * a, 0.0);

    //allocate space for the Bessel functions of the first and second kind (and derivatives)
    int bytes = sizeof(complex<double>) * (Nl + 1);
    complex<double>* cjv_ka = (complex<double>*)malloc(bytes);
    complex<double>* cyv_ka = (complex<double>*)malloc(bytes);
    complex<double>* cjvp_ka = (complex<double>*)malloc(bytes);
    complex<double>* cyvp_ka = (complex<double>*)malloc(bytes);
    complex<double>* cjv_kna = (complex<double>*)malloc(bytes);
    complex<double>* cyv_kna = (complex<double>*)malloc(bytes);
    complex<double>* cjvp_kna = (complex<double>*)malloc(bytes);
    complex<double>* cyvp_kna = (complex<double>*)malloc(bytes);

    //allocate space for the spherical Hankel functions and derivative
    complex<double>* chv_ka = (complex<double>*)malloc(bytes);
    complex<double>* chvp_ka = (complex<double>*)malloc(bytes);

    //compute the bessel functions using the CPU-based algorithm
    double vm;
    cbessjyva_sph(Nl, ka, vm, cjv_ka, cyv_ka, cjvp_ka, cyvp_ka);
    cbessjyva_sph(Nl, kna, vm, cjv_kna, cyv_kna, cjvp_kna, cyvp_kna);


    //cout<<"Begin Sphere---------"<<endl;
    //cout<<"Nl =  "<<Nl<<endl;
    //cout<<"ka =  "<<ka<<endl;
    //cout<<"kna = "<<kna<<endl;

    //compute A for each order
    complex<double> i(0, 1);
    complex<double> a, b, c, d;
    complex<double> An, Bn;
    for(int l=0; l<=Nl; l++)
    {
        //compute the Hankel functions from j and y
        chv_ka[l] = cjv_ka[l] + i * cyv_ka[l];
        chvp_ka[l] = cjvp_ka[l] + i * cyvp_ka[l];

        //Compute A (internal scattering coefficient)
        //compute the numerator and denominator for A
        a = cjv_ka[l] * chvp_ka[l] - cjvp_ka[l] * chv_ka[l];
        b = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

        //calculate A and add it to the list
        An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
        A.push_back(bsComplex(An.real(), An.imag()));
        //cout<<"A:  "<<An<<endl;

        //Compute B (external scattering coefficient)
        c = cjv_ka[l] * cjvp_kna[l] * nc - cjv_kna[l] * cjvp_ka[l];
        d = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

        //calculate B and add it to the list
        Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
        B.push_back(bsComplex(Bn.real(), Bn.imag()));
        //cout<<"B:  "<<Bn<<endl;

    }
}