function reflexp(teta1,n1,n2)

% Plots the electric field amplitude at an interface between
% incident medium n1 and transmitted medium n2 for p-polarization

wavel=1;
teta1=teta1*pi/180;

k0=2*pi/wavel;
k1=k0*n1;
k2=k0*n2;

k1x=k1*sin(teta1);
k1z=k1*cos(teta1);

teta2=asin(n1*sin(teta1))/n2;

k2x=k1x;
k2z=k2*cos(teta2);

%p-polarization
t=2*n1*cos(teta1)/(n2*cos(teta1)+n1*cos(teta2))
r=(n2*cos(teta1)-n1*cos(teta2))/(n2*cos(teta1)+n1*cos(teta2))

x=0:0.02:5;
z=-3:0.02:3;

xx=[0 5];
yy=[1 1];
zz=[0 0];

[dum,nx]=size(x);
[dum,nz]=size(z);
for iz=1:nz
    if z(iz) < 0
        for ix=1:nx
            E(ix,iz)=exp(i*(k1x*x(ix)+k1z*z(iz)))+r*exp(i*(k1x*x(ix)-k1z*z(iz)));
% Uncommment the following lines to see the incident
% or reflected fields only:
            %E(ix,iz)=exp(i*(k1x*x(ix)+k1z*z(iz)));
            %E(ix,iz)=r*exp(i*(k1x*x(ix)-k1z*z(iz)));
            P(ix,iz)=real(E(ix,iz));
        end
    else
        for ix=1:nx
            E(ix,iz)=t*exp(i*(k2x*x(ix)+k2z*z(iz)));
            P(ix,iz)=real(E(ix,iz));
        end
    end
end
[X,Z]=meshgrid(x,z);
surf(X,Z,P');
hold on;
shading interp;
view([0,90]);
axis equal;
axis off;
h=line(xx,zz,yy);
set(h,'Color',[1 1 1],'LineWidth',2);
set(gca,'position',[0 0 1 1],'units','normalized');
hold off;
