halfspace.cuh 8.83 KB
``````#ifndef	RTS_HALFSPACE_H
#define	RTS_HALFSPACE_H

#include "../math/plane.h"

namespace stim{

//GPU kernel to compute the electric field
template<typename T>
__global__ void gpu_halfspace_pw2ef(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
plane<T> P, planewave<T> w, rect<T> q, bool at_surface = false)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;

//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;

//compute the index into the field
int i = iv*r0 + iu;

//get the current position
vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );

if(at_surface){
if(P.side(p) > 0)
return;
}
else{
if(P.side(p) >= 0)
return;
}

//if the current position is on the wrong side of the plane

//vec<T> r(p[0], p[1], p[2]);

complex<T> x( 0.0f, w.kvec().dot(p) );
//complex<T> phase( 0.0f, w.phase());

if(Y == NULL)                       //if this is a scalar simulation
X[i] += w.E().len() * exp(x);    //use the vector magnitude as the plane wave amplitude
else
{
//X[i] = Y[i] = Z[i] = 1;
X[i] += w.E()[0] * exp(x);// * exp(phase);
Y[i] += w.E()[1] * exp(x);// * exp(phase);
Z[i] += w.E()[2] * exp(x);// * exp(phase);
}
}

//GPU kernel to compute the electric displacement field
template<typename T>
__global__ void gpu_halfspace_pw2df(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
plane<T> P, planewave<T> w, rect<T> q, T n, bool at_surface = false)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;

//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;

//compute the index into the field
int i = iv*r0 + iu;

//get the current position
vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );

if(at_surface){
if(P.side(p) > 0)
return;
}
else{
if(P.side(p) >= 0)
return;
}

//if the current position is on the wrong side of the plane

//vec<T> r(p[0], p[1], p[2]);

complex<T> x( 0.0f, w.kvec().dot(p) );
//complex<T> phase( 0.0f, w.phase());

//vec< complex<T> > testE(1, 0, 0);

if(Y == NULL)                       //if this is a scalar simulation
X[i] += w.E().len() * exp(x);    //use the vector magnitude as the plane wave amplitude
else
{
plane< complex<T> > cplane = plane< complex<T>, 3>(P);
vec< complex<T> > E_para;// = cplane.parallel(w.E());
vec< complex<T> > E_perp;// = cplane.perpendicular(w.E()) * (n*n);
cplane.decompose(w.E(), E_para, E_perp);
T epsilon = n*n;

X[i] += (E_para[0] + E_perp[0] * epsilon) * exp(x);
Y[i] += (E_para[1] + E_perp[1] * epsilon) * exp(x);
Z[i] += (E_para[2] + E_perp[2] * epsilon) * exp(x);
}
}

//computes a scalar field containing the refractive index of the half-space at each point
template<typename T>
__global__ void gpu_halfspace_n(T* n, unsigned int r0, unsigned int r1, rect<T> q, plane<T> P, T n0, T n1){

int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;

//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;

//compute the index into the field
int i = iv*r0 + iu;

//get the current position
vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );

//set the appropriate refractive index
if(P.side(p) < 0) n[i] = n0;
else n[i] = n1;
}

template<class T>
class halfspace
{
private:
rts::plane<T> S;		//surface plane splitting the half space
rts::complex<T> n0;		//refractive index at the front of the plane
rts::complex<T> n1;		//refractive index at the back of the plane

//lists of waves in front (pw0) and behind (pw1) the plane
std::vector< rts::planewave<T> > w0;
std::vector< rts::planewave<T> > w1;

//rts::planewave<T> pi;	//incident plane wave
//rts::planewave<T> pr;	//plane wave reflected from the surface
//rts::planewave<T> pt;	//plane wave transmitted through the surface

void init(){
n0 = 1.0;
n1 = 1.5;
}

//compute which side of the interface is hit by the incoming plane wave (0 = front, 1 = back)
bool facing(planewave<T> p){
if(p.kvec().dot(S.norm()) > 0)
return 1;
else
return 0;
}

T calc_theta_i(vec<T> v){

vec<T> v_hat = v.norm();

//compute the cosine of the angle between k_hat and the plane normal
T cos_theta_i = v_hat.dot(S.norm());

return acos(abs(cos_theta_i));
}

T calc_theta_t(T ni_nt, T theta_i){

T sin_theta_t = ni_nt * sin(theta_i);
return asin(sin_theta_t);
}

public:

//constructors
halfspace(){
init();
}

halfspace(T na, T nb){
init();
n0 = na;
n1 = nb;
}

halfspace(T na, T nb, plane<T> p){
n0 = na;
n1 = nb;
S = p;
}

//compute the transmitted and reflective waves given the incident (vacuum) plane wave p
void incident(rts::planewave<T> p){

planewave<T> r, t;
p.scatter(S, n1.real()/n0.real(), r, t);

//std::cout<<"i+r: "<<p.r()[0] + r.r()[0]<<std::endl;
//std::cout<<"t:   "<<t.r()[0]<<std::endl;

if(facing(p)){
w1.push_back(p);

if(r.E().len() != 0)
w1.push_back(r);
if(t.E().len() != 0)
w0.push_back(t);
}
else{
w0.push_back(p);

if(r.E().len() != 0)
w0.push_back(r);
if(t.E().len() != 0)
w1.push_back(t);
}
}

void incident(rts::beam<T> b, unsigned int N = 10000){

//generate a plane wave decomposition for the beam
std::vector< planewave<T> > pw_list = b.mc(N);

//calculate the reflected and refracted waves for each incident wave
for(unsigned int w = 0; w < pw_list.size(); w++){
incident(pw_list[w]);
}
}

//return the electric field at the specified resolution and position
rts::efield<T> E(unsigned int r0, unsigned int r1, rect<T> R){

int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);

//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((r0 + SQRT_BLOCK -1)/SQRT_BLOCK, (r1 + SQRT_BLOCK - 1)/SQRT_BLOCK);

//create an electric field
rts::efield<T> ef(r0, r1);
ef.position(R);

//render each plane wave

//plane waves at the surface front
for(unsigned int w = 0; w < w0.size(); w++){
//std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.x(), ef.y(), ef.z(), r0, r1, S, w0[w], ef.p());
}

//plane waves at the surface back
for(unsigned int w = 0; w < w1.size(); w++){
//std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
rts::gpu_halfspace_pw2ef<T> <<<dimGrid, dimBlock>>> (ef.x(), ef.y(), ef.z(), r0, r1, -S, w1[w], ef.p(), true);
}

return ef;
}

//return the electric displacement at the specified resolution and position
rts::efield<T> D(unsigned int r0, unsigned int r1, rect<T> R){

int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);

//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((r0 + SQRT_BLOCK -1)/SQRT_BLOCK, (r1 + SQRT_BLOCK - 1)/SQRT_BLOCK);

//create a complex vector field
rts::efield<T> df(r0, r1);
df.position(R);

//render each plane wave

//plane waves at the surface front
for(unsigned int w = 0; w < w0.size(); w++){
//std::cout<<"w0 "<<w<<": "<<hs.w0[w].str()<<std::endl;
rts::gpu_halfspace_pw2df<T> <<<dimGrid, dimBlock>>> (df.x(), df.y(), df.z(), r0, r1, S, w0[w], df.p(), n0.real());
}

//plane waves at the surface back
for(unsigned int w = 0; w < w1.size(); w++){
//std::cout<<"w1 "<<w<<": "<<hs.w1[w].str()<<std::endl;
rts::gpu_halfspace_pw2df<T> <<<dimGrid, dimBlock>>> (df.x(), df.y(), df.z(), r0, r1, -S, w1[w], df.p(), n1.real(), true);
}

return df;
}

realfield<T, 1, true> nfield(unsigned int Ru, unsigned int Rv, rect<T> p){

int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);

//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((Ru + SQRT_BLOCK -1)/SQRT_BLOCK, (Rv + SQRT_BLOCK - 1)/SQRT_BLOCK);

realfield<T, 1, true> n(Ru, Rv);	//create a scalar field to store the result of n

rts::gpu_halfspace_n<T> <<<dimGrid, dimBlock>>> (n.ptr(), Ru, Rv, p, S, n0.real(), n1.real());

return n;
}

std::string str(){
std::stringstream ss;
ss<<"Half Space-------------"<<std::endl;
ss<<"n0: "<<n0<<std::endl;
ss<<"n1: "<<n1<<std::endl;
ss<<S.str();

if(w0.size() != 0 || w1.size() != 0){
ss<<std::endl<<"Plane Waves:"<<std::endl;
for(unsigned int w = 0; w < w0.size(); w++){
ss<<"w0 = "<<w<<": "<<w0[w]<<std::endl;
}
for(unsigned int w = 0; w < w1.size(); w++){
ss<<"w1 = "<<w<<": "<<w1[w]<<std::endl;
}
}

return ss.str();
}
////////FRIENDSHIP
//friend void operator<< <> (rts::efield<T> &ef, rts::halfspace<T> hs);

};

}

#endif``````