optics / imie

  __device__ double g(double v0, double v1)

  {
  	return (v0 + v1)*log(abs(v0+v1)) + (v0-v1)*log(abs(v0-v1));
  }
  

  __device__ double hfin(double v0, double v1, double dv)

  {

  	double e = 0.001;
  	double t0 = g(v0+e, v1-dv)/dv;
  	double t1 = 2*g(v0+e, v1)/dv;
  	double t2 = g(v0+e, v1+dv)/dv;

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

  __global__ void devKramersKronig(double* gpuN, double* gpuK, int numVals, double nuStart, double nuEnd, double nOffset)

  {
  	int i = blockIdx.x * blockDim.x + threadIdx.x;
  
  	if(i >= numVals) return;

  	double nuDelta = (nuEnd - nuStart)/(numVals - 1);

  	double nu = nuStart + i*nuDelta;
  	double n = 0.0;
  	double jNu;
  	double 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(double* cpuN, double* cpuK, int nVals, double nuStart, double nuEnd, double nOffset)

  {

  	double* gpuK;
  	HANDLE_ERROR(cudaMalloc(&gpuK, sizeof(double)*nVals));
  	HANDLE_ERROR(cudaMemcpy(gpuK, cpuK, sizeof(double)*nVals, cudaMemcpyHostToDevice));
  	double* gpuN;
  	HANDLE_ERROR(cudaMalloc(&gpuN, sizeof(double)*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(double)*nVals, cudaMemcpyDeviceToHost));

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

  __global__ void devComputeSpectrum(double* I, double2* B, double* alpha, int Nl, 
  								   int nSamples, double oThetaI, double oThetaO, double cThetaI, double cThetaO)

  {
  	int i = blockIdx.x * blockDim.x + threadIdx.x;
  
  	//compute the delta-theta value

  	double dTheta = (oThetaO - oThetaI)/nSamples;

  
  	//allocate space for the Legendre polynomials

  	double Ptheta[2];	
  
  	double cosTheta, theta;
  	cuDoubleComplex Us;
  	cuDoubleComplex UsSample;
  	cuDoubleComplex U;
  	//cuComplex Ui;
  	//Ui.x = 2*PI;
  	//Ui.y = 0.0;
  	cuDoubleComplex numer;

  	numer.x = 0.0;

  	cuDoubleComplex exp_numer;
  	cuDoubleComplex iL;
  	cuDoubleComplex imag;

  	imag.x = 0.0; imag.y = 1.0;

  	double realFac;
  	cuDoubleComplex complexFac;
  	double PlTheta;
  	double 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);
  

  		}
  
  		//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(double* cpuI, double* cpuB, double* cpuAlpha,
  						 int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples)

  {
  	//copy everything to the GPU

  	double2* gpuB;
  	HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(double2) * nLambda * Nl));
  	HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(double2) * nLambda * Nl, cudaMemcpyHostToDevice));

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

  	double* gpuI;
  	HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(double) * 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, (double2*)gpuB, gpuAlpha, Nl,

  										nSamples, oThetaI, oThetaO, cThetaI, cThetaO);
  

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

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