__device__ float g(float v0, float v1)
{
	return (v0 + v1)*log(abs(v0+v1)) + (v0-v1)*log(abs(v0-v1));
}

__device__ float hfin(float v0, float v1, float dv)
{
	float e = 0.001;
	float t0 = g(v0+e, v1-dv)/dv;
	float t1 = 2*g(v0+e, v1)/dv;
	float t2 = g(v0+e, v1+dv)/dv;

	return -1.0/PI * (t0 - t1 + t2);
}

__global__ void devKramersKronig(float* gpuN, float* gpuK, int numVals, float nuStart, float nuEnd, float nOffset)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;

	if(i >= numVals) return;
	float nuDelta = (nuEnd - nuStart)/(numVals - 1);

	float nu = nuStart + i*nuDelta;
	float n = 0.0;
	float jNu;
	float jK;
	for(int j=1; j<numVals-1; j++)
	{
		jNu = nuStart + j*nuDelta;
		jK = gpuK[j];
		n += hfin(nu, jNu, nuDelta) * jK;
	}
	gpuN[i] = n + nOffset;


}

void cudaKramersKronig(float* cpuN, float* cpuK, int nVals, float nuStart, float nuEnd, float nOffset)
{
	float* gpuK;
	HANDLE_ERROR(cudaMalloc(&gpuK, sizeof(float)*nVals));
	HANDLE_ERROR(cudaMemcpy(gpuK, cpuK, sizeof(float)*nVals, cudaMemcpyHostToDevice));
	float* gpuN;
	HANDLE_ERROR(cudaMalloc(&gpuN, sizeof(float)*nVals));

	dim3 block(BLOCK_SIZE*BLOCK_SIZE);
	dim3 grid(nVals/block.x + 1);
	devKramersKronig<<<grid, block>>>(gpuN, gpuK, nVals, nuStart, nuEnd, nOffset);

	HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(float)*nVals, cudaMemcpyDeviceToHost));

	//free resources
	HANDLE_ERROR(cudaFree(gpuK));
	HANDLE_ERROR(cudaFree(gpuN));
}

__global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl, 
								   int nSamples, float oThetaI, float oThetaO, float cThetaI, float cThetaO)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;

	//compute the delta-theta value
	float dTheta = (oThetaO - oThetaI)/nSamples;

	//allocate space for the Legendre polynomials
	float Ptheta[2];	

	float cosTheta, theta;
	cuComplex Us;
	cuComplex UsSample;
	cuComplex U;
	cuComplex Ui;
	Ui.x = 2*PI;
	Ui.y = 0.0;
	cuComplex numer;
	numer.x = 0.0;
	cuComplex exp_numer;
	cuComplex iL;
	cuComplex imag;
	imag.x = 0.0; imag.y = 1.0;
	float realFac;
	cuComplex complexFac;
	float PlTheta;
	float Isum = 0.0;
	float maxVal = 0;
	float val;
	for(int iTheta = 0; iTheta < nSamples; iTheta++)
	{
		//calculate theta
		theta = iTheta * dTheta + oThetaI;
		cosTheta = cos(theta);

		//initialize the theta Legendre polynomial
		Ptheta[0] = 1.0;
		Ptheta[1] = cosTheta;

		//initialize the scattered field
		Us.x = Us.y = 0.0;
		iL.x = 1.0;
		iL.y = 0.0;
		for(int l = 0; l<Nl; l++)
		{
			//compute the theta legendre polynomial
			if(l == 0)
				PlTheta = Ptheta[0];
			else if(l == 1)
				PlTheta = Ptheta[1];
			else
			{
				PlTheta = ((2*l - 1)*cosTheta*Ptheta[1] - (l - 1)*Ptheta[0])/l;
				Ptheta[0] = Ptheta[1];
				Ptheta[1] = PlTheta;
			}

			//compute the real components of the scattered field
			realFac = alpha[l] * PlTheta;

			//compute the complex components of the scattered field
			numer.x = 0.0;
			numer.y = -(l*PI)/2.0;
			exp_numer = cExp(numer);

			complexFac = cMult(B[Nl * i + l], exp_numer);
			complexFac = cMult(complexFac, iL);
			

			//combine the real and complex components
			UsSample = cMult(complexFac, realFac);
			Us = cAdd(Us, UsSample);

			//increment the imaginary exponent i^l
			iL = cMult(iL, imag);

			//val = cMag(Us);
			//if(val > maxVal)
			//	maxVal = val;

			//Us += B[l] * exp(numer) * Ptheta[l] * alpha * M * pow(complex<float>(0.0, 1.0), l);
			
		}

		//sum the scattered and incident fields
		if(theta >= cThetaI && theta <= cThetaO)
			U = cAdd(Us, 2*PI);
		else
			U = Us;
		Isum += (U.x*U.x + U.y*U.y) * sin(theta) * 2 * PI * dTheta;
	}

	I[i] = Isum;
}

void cudaComputeSpectrum(float* cpuI, float* cpuB, float* cpuAlpha,
						 int Nl, int nLambda, float oThetaI, float oThetaO, float cThetaI, float cThetaO, int nSamples)
{
	//copy everything to the GPU
	float2* gpuB;
	HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(float2) * nLambda * Nl));
	HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(float2) * nLambda * Nl, cudaMemcpyHostToDevice));

	float* gpuAlpha;
	HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(float) * Nl));
	HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(float) * Nl, cudaMemcpyHostToDevice));

	float* gpuI;
	HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(float) * nLambda));
	

	//call the kernel to compute the spectrum
	dim3 block(BLOCK_SIZE*BLOCK_SIZE);
	dim3 grid(nLambda/block.x + 1);

	//devComputeSpectrum
	devComputeSpectrum<<<grid, block>>>(gpuI, (float2*)gpuB, gpuAlpha, Nl,
										nSamples, oThetaI, oThetaO, cThetaI, cThetaO);

	HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(float) * nLambda, cudaMemcpyDeviceToHost));

	HANDLE_ERROR(cudaFree(gpuB));
	HANDLE_ERROR(cudaFree(gpuAlpha));
	HANDLE_ERROR(cudaFree(gpuI));

	


}