
function [M]=matriceCouche(lambda,ni,thetai,pola,n,e)
    cosinus = sqrt(1-(ni/n*sin(thetai))^2);
    if pola==1 then
        p = n*cosinus;
    else
        p = cosinus /n;
    end,
    phi = 2*%pi/lambda*n*e*cosinus;
    M=[cos(phi),%i*sin(phi)/p;%i*p*sin(phi),cos(phi)]
endfunction
            
function [M]=matriceMultiCouches(lambda,ni,thetai,pola,couches)
    M = [1,0;0,1];
    s = size(couches);
    nc = s(1);
    for k=1:nc,
        M = M*matriceCouche(lambda,ni,thetai,pola,couches(k,1),couches(k,2));
    end,
endfunction
            
function [r,t]=coefficientsRT(lambda,ni,thetai,nt,pola,couches)
    if pola==1 then 
        pi = ni*cos(thetai);
        pt = nt*sqrt(1-(ni/nt*sin(thetai))^2);
    else 
        pi = cos(thetai)/ni;
        pt = sqrt(1-(ni/nt*sin(thetai))^2)/nt;
    end,
    m = matriceMultiCouches(lambda,ni,thetai,pola,couches);
    a = (m(1,1)+m(1,2)*pt)*pi;
    b = m(2,1)+m(2,2)*pt;
    r = (a-b)/(a+b);
    if pola==1 then
        t = 2*pi/(a+b);
    else
        t = 2*pi/(a+b)*ni/nt;
    end,
endfunction

            
function [R,T]=coefficientsRTpuissances(lambda,ni,thetai,nt,pola,couches)
    [r,t] = coefficientsRT(lambda,ni,thetai,nt,pola,couches);
    R=abs(r)^2;
    if pola==1 then 
            pi = ni*cos(thetai);
            pt = nt*sqrt(1-(ni/nt*sin(thetai))^2);
    else 
            pi = cos(thetai)/ni;
            pt = sqrt(1-(ni/nt*sin(thetai))^2)/nt;
    end,
    T = pt/pi*abs(t)^2;
endfunction 
            
function [r,t]=coefRT(lambda,ni,listThetai,nt,pola,couches)
    np = length(listThetai);
    r = zeros(np);
    t = zeros(np);
    for i=1:np,
        [r(i),t(i)] = coefficientsRT(lambda,ni,listThetai(i),nt,pola,couches);
    end,
endfunction
            
function [r,t]=coefRTpuissances(lambda,ni,listThetai,nt,pola,couches)
    np = length(listThetai);
    r = zeros(np);
    t = zeros(np);
    for i=1:np,
        [r(i),t(i)] = coefficientsRTpuissances(lambda,ni,listThetai(i),nt,pola,couches);
    end,
endfunction
            