epsilon0=8.85e-12;
Ne=1e12;
e=1.6e-19;
me=9.1e-31;
c=3e8;
fc=sqrt(Ne*e^2/(epsilon0*me))/(2*%pi)
            

delta=c/(2*%pi)/sqrt(fc^2*(1-0.999^2))
            

function R=reflexionTE(n1,n2,theta1)
    cosTheta2 = sqrt(1-(n1*sin(theta1)/n2)^2);
    R = abs((n1*cos(theta1)-n2*cosTheta2)/(n1*cos(theta1)+n2*cosTheta2))^2
endfunction
function R=reflexionTM(n1,n2,theta1)
    cosTheta2 = sqrt(1-(n1*sin(theta1)/n2)^2);
    R = abs((n2*cos(theta1)-n1*cosTheta2)/(n2*cos(theta1)+n1*cosTheta2))^2
endfunction
function n=indice(f)
    n = sqrt(1-1/f^2)
endfunction
            

f = linspace(0.1,3,200);
np = length(f);
R = zeros(np);
for i=1:np,
    R(i)=reflexionTE(1,indice(f(i)),0);
end,
            
plotA=scf();
plot2d(f,R,style=5);
xtitle("Incidence theta = 0","f/fc","R");
            

RTE=zeros(np);
RTM=zeros(np);
for i=1:np,
    RTE(i)=reflexionTE(1,indice(f(i)),45/180*%pi);
    RTM(i)=reflexionTM(1,indice(f(i)),45/180*%pi);
end,
            
plotB=scf();
plot2d(f,RTE,style=4);
plot2d(f,RTM,style=5);
legend("TE","TM");
xtitle("Incidence theta = 45","f/fc","R");
            

function r=reflexion(f,theta)
    k=2*%pi*f;
    h=100;
    alpha2=sin(theta)^2;
    function n2=indice2(z)
        n2=1-1/f^2*(1-4/h^2*(z-h/2)^2)
    endfunction
    function [deriv]=systeme(z,y)
        deriv(1)=k*y(4),
        deriv(2)=-k*y(3),
        deriv(3)=k*(indice2(z)-alpha2)*y(2),
        deriv(4)=-k*(indice2(z)-alpha2)*y(1),
    endfunction 
    z=[0,h];
    tolA=1d-4;
    tolR=1d-10;
    %ODEOPTIONS=[1,0,0,%inf,0,0,1e9,12,5,0,-1,-1];
    V1=ode("stiff",[1;0;0;0],0,z,tolR,tolA,systeme);
    V2=ode("stiff",[0;0;1;0],0,z,tolR,tolA,systeme);
    p=length(z);
    M=[V2(3,p)+%i*V2(4,p), -V2(1,p)-%i*V2(2,p); -V1(3,p)-%i*V1(4,p), V1(1,p)+%i*V1(2,p)];
    r=((M(1,1)+M(1,2))-(M(2,1)+M(2,2)))/((M(1,1)+M(1,2))+(M(2,1)+M(2,2)))
endfunction
            
f=linspace(0.5,1.5,300);
np=length(f);
R=zeros(np);
angle=0;
for i=1:np,
    R(i)=abs(reflexion(f(i),angle))^2;
end,
            
plotC=scf();
plot2d(f,R,style=4);
plot2d(f,R,style=5);
xtitle("Incidence theta = 0","f/fc","R");
            
f=linspace(1,2,1000);
np=length(f);
R=zeros(np);
angle=45*%pi/180;
for i=1:np,
    R(i)=abs(reflexion(f(i),angle))^2;
end,
            
plotD=scf();
plot2d(f,R,style=4);
plot2d(f,R,style=5);
xtitle("Incidence theta = 45 (onde TE)","f/fc","R");
            
