#include <complex>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include "rts/complex.h"
#include "rts/sbessel.h"

#include "compare.h"

typedef float precision;

//typedef std::complex<double> ptype;
typedef rts::complex<precision> ptype;

void cpuValidateBessel()
{

	//order
	precision v = 5;
	precision vm;

	//parameter
	precision max_z = 5;
	int nz = 20;
	precision dz = max_z / nz;

	//bessel function results (first and second kind)
	int S = sizeof(ptype) * (v + 1);

	ptype* jv = (ptype*)malloc(S);
	ptype* yv = (ptype*)malloc(S);
	ptype* hv = (ptype*)malloc(S);


    std::cout<<"---------j_v(x)-------------"<<std::endl;
	rts::complex<precision> z;
	for(int iz = 0; iz < nz; iz++)
	{
		z = iz * dz;

		rts::sbesselj<precision>(v, z, jv);

		std::cout<<z.toStr()<<", "<<jv[0].toStr()<<", "<<jv[1].toStr()<<", "<<jv[2].toStr()<<std::endl;

	}

	std::cout<<"---------y_v(x)-------------"<<std::endl;
	for(int iz = 0; iz < nz; iz++)
	{
		z = iz * dz;

		rts::sbessely<precision>(v, z, yv);

		std::cout<<z.toStr()<<", "<<yv[0].toStr()<<", "<<yv[1].toStr()<<", "<<yv[2].toStr()<<std::endl;

	}

	std::cout<<"---------h(1)_v(x)-------------"<<std::endl;
	for(int iz = 0; iz < nz; iz++)
	{
		z = iz * dz;

		rts::sbesselh1<precision>(v, z, hv);

		std::cout<<z.toStr()<<", "<<hv[0].toStr()<<", "<<hv[1].toStr()<<", "<<hv[2].toStr()<<std::endl;

	}








}