theta_i = 30; % in degrees

lambda = .6328; % in microns



theta_i = pi/180*theta_i;



n1 = 1;

n2 = 1.5;

n21 = n2/n1;

theta_t = asin(n1/n2*sin(theta_i)); % Snell's law



k1 = 2*pi*n1/lambda;

k2 = 2*pi*n2/lambda;

k_i = k1*[sin(theta_i); -cos(theta_i)];

k_r = k1*[sin(theta_i);  cos(theta_i)];

k_t = k2*[sin(theta_t); -cos(theta_t)];



[x y1] = meshgrid(-2:0.04:2,0:0.04:2);

[x y2] = meshgrid(-2:0.04:2,-2:0.04:0);



N = 3;

T = linspace(0,2*pi*N,N*100);

figure;    



r = (cos(theta_i)-sqrt(n21^2-sin(theta_i)^2))/(cos(theta_i)+sqrt(n21^2-sin(theta_i)^2));

t = 1+r;



for k = 1:length(T)

    E_i = exp(i*(k_i(1)*x+k_i(2)*y1-T(k)));

    E_i = real(E_i);



    E_r = r*exp(i*(k_r(1)*x+k_r(2)*y1-T(k)));

    E_r = real(E_r);



    E_t = t*exp(i*(k_t(1)*x+k_t(2)*y2-T(k)));

    E_t = real(E_t);



    surf(x,y1,E_i+E_r);

    hold on;

    surf(x,y2,E_t);

    xlabel('x');

    ylabel('y');

    ylim([-2 2]);

    title('TE boundary conditions n_1 = 1, n_2 = 1.5');

    zlabel('Transverse electric field (continuous)'); view(2);

    drawnow;

    hold off;

end