#ifndef	RTS_HALFSPACE_H
#define	RTS_HALFSPACE_H

namespace rts{
template<typename T> class halfspace;

template<typename T> class efield;
}

template<typename T>
void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);

namespace rts{

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]);
		}
	}

	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