propagation d'ondes planes dans une cavité parallélépipédique

[jackboer63]

Salut à tous,

Je suis nouveau sur le forum. Je sollicite votre aide pour programmer la propagation d'une onde plane dans une cavité électromagnétique par la méthode FDTD.
Je ne sais pas comment programmer une onde plane avec des paramètres en entrée réglables. Je souhaiterais avoir une onde plane de type Gaussienne 3D et temporelle qui se propage à travers les mailles de calcul.
Dans un second temps, j'aimerais appliquer le formalisme champ total/diffracté (TFSF en anglais) en programmant une surface fictive :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%VECTORISE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

% Constantes physiques
eps0 = 8.8541878e-12;        % Permittivite du vide
mu0  = 4e-7 * pi;            % Permeabilite du vide
c0  = 299792458;            % Vitesse de la lumiere dans le vide

% Initialisation parametrique
Lx = 6.7; Ly = 8.4; Lz = 3.5; % Dimensions de la cavite en metres
Nx =  60; Ny =  60; Nz =  60; % Nombre de cellules dans chaque direction
Cx = Nx / Lx;                % Dimensions inverses des cellules
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 1000;                    % Nombre de pas temporels
Dt = 1/(c0*norm([Cx Cy Cz])); % Definition du pas temporel
Dt = Dt*0.9;

% Allocation des champs de matrices
Ex = zeros(Nx  , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny  , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz  );
Hx = zeros(Nx+1, Ny  , Nz  );
Hy = zeros(Nx  , Ny+1, Nz  );
Hz = zeros(Nx  , Ny  , Nz+1);

% Déclaration de la variable à la même taille que t
Ets = zeros(Nt+1,3);
Ets1 = zeros(Nt+1,3);
Ets2 = zeros(Nt+1,3);

% Initialisation des champs (pas sur les bords)
Ex(1,2,2) = 0;
Ey(2,1,2) = 0;
Ez(2,2,1) = 0;

tic

% Boucle pour le remplissage des matrices
for n = 1:Nt;
    for X = 1:Nx;
        for Y = 1:Ny;
            for Z = 1:Nz;
               
        %source d'excitation Gaussienne
        t0 = 2;
        tau = 5; 
        x0 = -0.5;
        y0 = -0.5;
        z0 = -0.5;
        tc = n +(X+x0)/c0 + (Y+y0)/c0 + (Z+z0)/c0;
        gama = (tc+t0)/tau;
        Ex(30,30,30) = 377.*exp(-gama.*gama);
   
            end;
        end;
    end;
       
  % Mise a jour de H partout
  Hx = Hx + (Dt/mu0)*((Ey,:,2:Nz+1)-Ey,:,1:Nz))*Cz ...
                    - (Ez,2:Ny+1,-Ez,1:Ny,)*Cy);
 
  Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,-Ez(1:Nx,:,)*Cx ...
                    - (Ex,:,2:Nz+1)-Ex,:,1:Nz))*Cz);
               
  Hz = Hz + (Dt/mu0)*((Ex,2:Ny+1,-Ex,1:Ny,)*Cy ...
                    - (Ey(2:Nx+1,:,-Ey(1:Nx,:,)*Cx);
 
 
  % Mise a jour de E partout sauf sur les bords
  Ex,2:Ny,2:Nz) = Ex,2:Ny,2:Nz) + (Dt /eps0) * ...
      ((Hz,2:Ny,2:Nz)-Hz,1:Ny-1,2:Nz))*Cy ...
    - (Hy,2:Ny,2:Nz)-Hy,2:Ny,1:Nz-1))*Cz);

  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
      ((Hx(2:Nx,:,2:Nz)-Hx(2:Nx,:,1:Nz-1))*Cz ...
    - (Hz(2:Nx,:,2:Nz)-Hz(1:Nx-1,:,2:Nz))*Cx);

  Ez(2:Nx,2:Ny, = Ez(2:Nx,2:Ny, + (Dt /eps0) * ...
      ((Hy(2:Nx,2:Ny,-Hy(1:Nx-1,2:Ny,)*Cx ...
    - (Hx(2:Nx,2:Ny,-Hx(2:Nx,1:Ny-1,)*Cy);
 
  % Echantillon de champ electrique dans un point precis du DC
  Ets(n, = Ex(10,10,10);% Remplissage de la matrice/vecteur avec les valeurs de Ex, Ey, ...
  Ets1(n, = Ey(10,10,10);
  Ets2(n, = Ez(10,10,10);
  t = [0t:Nt*Dt]; 
end

  %Affichage des résultats
  figure(1)
  plot(t,Ets,'b');
  legend('Ex(10,10,10)')
  title('Allure temporelle de la source Gaussienne')
  xlabel('temps(s)')
  ylabel('Amplitude du champ electrique (V/m)')
 
  figure(2)
  plot(t,Ets1,'r',t,Ets2,'g');
  legend('Ey(10,10,10)','Ez(10,10,10)')
  title('Propagation de la Gaussienne à travers le domaine de calcul')
  xlabel('temps(s)')
  ylabel('Amplitude du champ electrique (V/m)')
 
toc



S'il y a des connaisseurs en électromagnétisme cela serait super

Je vous donne mon adresse mail : jackboer63@gmail.com

Merci infiniment