esphere.cuh 5.33 KB
``````#ifndef RTS_ESPHERE
#define RTS_ESPHERE

#include "../math/complex.h"
#include "../math/bessel.h"
#include "../visualization/colormap.h"
#include "../optics/planewave.h"
#include "../cuda/devices.h"
#include "../optics/efield.cuh"

namespace stim{

/*  This class implements a discrete representation of an electromagnetic field
in 2D scattered by a sphere. This class implements Mie scattering.
*/
template<typename P>
class esphere : public efield<P>
{
private:
stim::complex<P> n;		//sphere refractive index
P a;					//sphere radius

//parameters dependent on wavelength
unsigned int Nl;		//number of orders for the calculation
rts::complex<P>* A;		//internal scattering coefficients
rts::complex<P>* B;		//external scattering coefficients

void calcNl(P kmag)
{
//return ceil( ((P)6.282 * a) / lambda + 4 * pow( ((P)6.282 * a) / lambda, ((P)1/(P)3)) + 2);
Nl = ceil( kmag*a + 4 * pow(kmag * a, (P)1/(P)3) + 2);
}

void calcAB(P k, unsigned int Nl, rts::complex<P>* A, rts::complex<P>* B)
{
/*  These calculations require double precision, so they are computed
using doubles and converted to P at the end.

Input:

k	=	magnitude of the k vector (tau/lambda)
Nl	=	highest order coefficient ([0 Nl] are computed)
*/

//clear the previous coefficients
rts::complex<P>* cpuA = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
rts::complex<P>* cpuB = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));

//convert to an std complex value
complex<double> nc = (rts::complex<double>)n;

//compute the magnitude of the k vector
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)
unsigned 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);

//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
rts::complex<double> An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
cpuA[l] = (rts::complex<P>)An;

//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
rts::complex<double> Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
cpuB[l] = (rts::complex<P>)Bn;

std::cout<<"A: "<<cpuA[l]<<"     B: "<<cpuB[l]<<std::endl;
}

if(A != NULL)	cudaFree(A);	//free any previous coefficients
if(B != NULL)	cudaFree(B);
cudaMalloc(&A, sizeof(rts::complex<P>) * (Nl+1));	//allocate memory for new coefficients
cudaMalloc(&B, sizeof(rts::complex<P>) * (Nl+1));

//copy the calculations from the CPU to the GPU
cudaMemcpy(A, cpuA, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
cudaMemcpy(B, cpuB, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
}

public:

esphere(unsigned int res0, unsigned int res1, P _a, rts::complex<P> _n, bool _scalar = false) :
efield<P>(res0, res1, _scalar)
{
std::cout<<"Sphere scattered field created."<<std::endl;
n = _n;		//save the refractive index
a = _a;		//save the radius
}

//assignment operator: build an electric field from a plane wave
efield<P> & operator= (const planewave<P> & rhs)
{
calcNl(rhs.kmag());		//compute the number of scattering coefficients
std::cout<<"Nl: "<<Nl<<std::endl;
calcAB(rhs.kmag(), Nl, A, B);	//compute the scattering coefficients

//determine important parameters for the scattering domain
unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
//unsigned int thetaR = 256;

/////////////////////continue scattering code here/////////////////////////

esphere::clear();				//clear any previous field data
from_planewave(rhs);	//create a field from the planewave
return *this;			//return the current object
}

string str()
{
stringstream ss;
ss<<"Mie Scattered Field"<<std::endl;
ss<<(*this).efield<P>::str()<<std::endl;
ss<<"a = "<<a<<std::endl;
ss<<"n = "<<n;

return ss.str();
}

};

}	//end namespace rts

#endif``````