MIRAGES : SIMULATION DANS UN MILIEU NON PERTURBE

function [ xf, zf, af ] = snell_descartes( xi, zi, ai )
%SNELL_DESCARTES float * float * float -> float*float*float
%   Retourne une position et un angle refracte a partir d une position
% initiale et un angle  incident

dx = 0.0001; %pas pour le calcul
ni = indice_exp(zi); %indice en zi

dz = dx/tan(ai); %calcul de l'incrémentage de l'altitude
zf = zi - dz; %altitude final
nf = indice_exp(zf); %indice optique à l'altitude zf

af=asin(ni*sin(ai)/nf);%calcul de l'ange refracte grâce a la loi de Snell-Descartes

xf = xi + dx; %abscisse final

end

%-------------------------------------------------------------------------%
%        Simulation de propagation du laser dans le mileu non perturbé
                     

z = 0.066; %z initial
x=0; %abscisse initiale
a= pi/2 - 0.018; %angle incident initial (interface horizontale)
xmax = 0.6; %limite de la cuve en abscisse
zmax = 0.1; %limite de la cuve en altitude
res(1,1)=x; %
res(1,2)=z; %
res(1,3)=a; %
i=2;

while and(x<xmax, and(z>0,z<zmax))
    [x,z,a]=snell_descartes(x,z,a); %calcul des nouvelles valeurs
    res(i,1)=x;     %remplissage
    res(i,2)=z;
    res(i,3)=a;
    i = i +1; 
end
plot(res(:,1), res(:,2)) %affichage de l'altitude en fonction de x
title('Height of the theoritical beam as a function of x')
xlabel('x (m)')
ylabel('Height (m)')

function [ n ] = indice_exp( z )
%INDICE_EXP float -> float
%   Retourne l'indice optique n correspondant à l'altitude z


load('ft_ind.mat'); %chargement de l'ajustement de l'indice pour le calcul

n=ft_ind(z);
end

function [ dz ] = devfite( x )
%DEVFITE float -> float
%   Retourne la valeur de la derivee de l'approximation de l'altitude en x

       p1 =  -1.041e+04;  %(-2.659e+04, 5763)
       p2 =   1.972e+04;  %(-6919, 4.636e+04)
       p3 =  -1.496e+04 ; %(-3.337e+04, 3460)
       p4 =        5926 ; %(-1012, 1.286e+04)
       p5 =       -1345 ; %(-2893, 201.9)
       p6 =         180 ; %(-27.59, 387.7)
       p7 =      -14.18  ;%(-30.47, 2.11)
       p8 =      0.5147 ; %(-0.1748, 1.204)
       p9 =    -0.02478  ;%(-0.03808, -0.01147)
       p10 =     0.06588 ;%(0.0658, 0.06597)

       
dz = 9*p1*x^8 + 8*p2*x^7 + 7*p3*x^6 + 6*p4*x^5 + 5*p5*x^4 + 4*p6*x^3 + 3*p7*x^2 + 2*p8*x + p9;
end