#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