% Propagation d'un paquet d'ondes (ondelette) dans un milieu peu dense
%
% Licence de Physique, UL5A
% Universite Paul Sabatier
% Licence de Physique UL5, TD numeriques, 2003-2004
% Version 08/02/2004 LK
% -----------------------------------------------

clear all

% Votre nom
% ---------
initiales='??';

% constantes
% ----------
% pi = 3.14159265359;
e= 1.60217733e-19;                    % moins la charge de l'electron
m= 0.91093897e-30;                    % masse de l'electron
c= 299792458;                         % Celerite de la lumiere
epsilon0= 1e7/(4*pi*c^2);

% parametres du probleme (a completer)
% ----------------------------
% Pulsation de rappel de l'electron
w1 = ??  % w1= 0 (ionosph�re), ou w1= 1.6e15 radians/s (milieu dielectrique)
w1_2 = w1^2;

% nombre d'electrons par m3
n= ?? % n=1.22e12  (ionosph�re), ou  n=1e24 (milieu dielectrique)

% pulsation de plasma
wp_2 = e^2*n/(epsilon0*m);
wp  = sqrt (wp_2)

W_2 = w1_2 + wp_2;                     % Grand Omega W^2 = w1^2 + wp^2
W   = sqrt (W_2)

% choix de la pulsation de l'onde a etudier
w0= W * ??                                % Pulsation centrale
w0_2= w0^2;

Dw= 0.05 * w0                              % Delta w
Cte= 1/Dw/sqrt(2*pi);

% Nombre d'onde, vitesse de phase et vitesse de groupe (a completer)
% ------------------------------------------------------------------
vphase = c * sqrt((w0_2 - w1_2) / (w0_2 - W_2))    % Vitesse de phase
k0 = w0 / vphase                                   % Nb d'onde
lambda0= 2*pi / k0				  % Longueur d'onde

vgroupe= vphase/(1 + w0_2*(1/(w1_2 - w0_2) - 1/(W_2 - w0_2)))
                                       
Dt= 4*pi/Dw                   % Temps que prend le paquet pour passer en 1 pt
Dz= vgroupe/Dw                % largeur du paquet en z

ti= 0                     % Temps initial
tf= Dt                    % Temps final

Nt= 2					  % Nb de fois de calcul de l'integrale
dt= (tf-ti)/(Nt-1)        % Pas de temps


% Calcul  de l'integration (a completer)
%----------------------------------------

% increment en temps
% --------------------
for lt= 1:Nt
  t(lt)= ti + (lt-1)*dt;
end

% Definition de la grille en z
% ----------------------------
zi= -4*Dz                    % grille en z a gauche
zf= 20*Dz                    % grille en z a droite
dz= lambda0/10               % pas en z
Nz= fix((zf-zi)/dz)+1

for lz= 1:Nz
  z(lz)= zi+(lz-1)*dz;
end

% Definition des pulsations
% -------------------------
wi= w0-4*Dw              % limite d'integration inferieure
if (wi<W & wi>w1)
  wi=W+1e-5
end
if (wi < 0)
  wi=0;
end

wf=w0+4*Dw              % Limite d'integration superieure
if (wf >w1 & wf<W)
  wf=w1-1e-5
end

% Integration (a completer) 
% -------------------------
Psi= zeros(Nz,Nt);
dw= Dw/20              % pour l'integrale

for w=wi:dw:wf
  w_2= w^2;
  C= Cte*exp(-((w-w0)/Dw)^2/2);
  k = (w / c) * sqrt ((W_2 - w_2) / (w1_2 - w_2));    % Nb d'onde

  for lt=1:Nt
  for lz=1:Nz
      Psi(lz, lt) = ??;
  end
  end
end

% Partie graphique
% ----------------
plot(z, Psi(:,1), z, Psi(:,2))
xlabel('Distance')
ylabel('Psi')
title(['w_1=',num2str(w1), ', w_0=', num2str(w0), ...
       ', Dw=', num2str(Dw), '/', initiales])

