FD-TD Leap Frog function in matlab

  1. function [Ezo, Hxo, Hyo] = leapFrog2D(Ezi, Hxi, Hyi, parameters)
  2.  
  3.     Ezo = Ezi;
  4.     Hxo = Hxi;
  5.     Hyo = Hyi;
  6.  
  7.  
  8.     c1 = parameters.c1;
  9.     c2 = parameters.c2;
  10.     c3 = parameters.c3;
  11.     c4 = parameters.c4;
  12.     dx = parameters.dx;
  13.     dy = parameters.dy;
  14.     boundary = parameters.boundary; 
  15.     wire = parameters.wire;
  16.     
  17.    
  18.     Nx = size(Ezi,1);
  19.     Ny = size(Ezi,2);
  20.     
  21.     j_p_0_E = 1:Nx-1;
  22.     j_p_0_E = 1:Ny-1;
  23.  
  24.     i_p_1_E = 2:Nx;
  25.     j_p_1_E = 2:Ny;
  26.  
  27.     i_p_0_H = 1:Nx-1;
  28.     j_p_0_H = 1:Ny-1;
  29.  
  30.     i_p_1_H = 2:Nx;
  31.     j_p_1_H = 2:Ny;
  32.  
  33.    
  34.     % transverse electric
  35.   
  36.     Ezo(i_p_1_E, j_p_1_E) = c1(i_p_1_E, j_p_1_E).*Ezi(i_p_1_E, j_p_1_E) + c2(i_p_1_E, j_p_1_E).*(  ( wire(i_p_1_E, j_p_1_E).*Hyi(i_p_1_E, j_p_1_E) - wire(i_p_0_E, j_p_1_E).*Hyi(i_p_0_E, j_p_1_E) )/dx - ( wire(i_p_1_E, j_p_1_E).*Hxi(i_p_1_E, j_p_1_E) - wire(i_p_1_E, j_p_0_E).*Hxi(i_p_1_E, j_p_0_E) )/dy   );
  37.  
  38.     Hxo(i_p_0_H, j_p_0_H) = Hxi(i_p_0_H, j_p_0_H) + c3(i_p_0_H, j_p_0_H).*(   ( wire(i_p_0_H, j_p_0_H).*Ezo(i_p_0_H, j_p_0_H) - wire(i_p_0_H, j_p_1_H).*Ezo(i_p_0_H, j_p_1_H) )/dy   );
  39.     Hyo(i_p_0_H, j_p_0_H) = Hyi(i_p_0_H, j_p_0_H) + c4(i_p_0_H, j_p_0_H).*(   ( wire(i_p_1_H, j_p_0_H).*Ezo(i_p_1_H, j_p_0_H) - wire(i_p_0_H, j_p_0_H).*Ezo(i_p_0_H, j_p_0_H) )/dx   );
  40. end