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, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO)

  {

  	int i = blockIdx.x * blockDim.x + threadIdx.x;
  	if(i >= nLambda-1)
          return;

  

  	//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));
  	HANDLE_ERROR(cudaMemset(gpuI, 0, 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, nLambda, oThetaI, oThetaO, cThetaI, cThetaO);

  

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

  

  	//printf("Final array value: %f\n", cpuI[nLambda-1]);

  

  	HANDLE_ERROR(cudaFree(gpuB));

  	HANDLE_ERROR(cudaFree(gpuAlpha));

  	HANDLE_ERROR(cudaFree(gpuI));

  

  

  }