% Schneider book pg. 196 first 2D sim
% using 'yee' algorithm
clear
sizeX = 161;  
sizeY = 161;
maxTime = 50;
si = 79; sj=79;

%sizeX = 20;  
%sizeY = 16;
%maxTime =3;
%si = 11; sj=9;

ez=zeros(sizeX, sizeY);
hx=zeros(sizeX, sizeY);
hy=zeros(sizeX, sizeY);

ezn=zeros(sizeX, sizeY);ezn1=zeros(sizeX, sizeY);
hxn=zeros(sizeX, sizeY);hxn1=zeros(sizeX, sizeY);
hyn=zeros(sizeX, sizeY);hyn1=zeros(sizeX, sizeY);



imp0 = 377;
Cdtds = 1/sqrt(2);


for n=1:maxTime  
    
    if mod(n,2) == 0    
       for i=3:2:sizeX-1
           for j=2:2:sizeY-1
               hx(i,j) = hx(i,j) - Cdtds*(ez(i,j+1) - ez(i,j-1)) / imp0;
           end
       end
    
        for i=2:2:sizeX-1
           for j=3:2:sizeY-1
               hy(i,j) = hy(i,j) +Cdtds*(ez(i+1,j) - ez(i-1,j)) / imp0;
           end
        end
        % apply hx and hy boundary conditions here
    end
    
    if mod(n,2) == 1
       for i=3:2:sizeX-1
           for j=3:2:sizeY-1
               ez(i,j) = ez(i,j)+Cdtds*(hy(i+1,j)-hy(i-1,j)- ... 
               (hx(i,j+1)-hx(i,j-1)))*imp0; 
           end
       end
       
       % apply ez boundary conditions here
       arg = pi * Cdtds*n / 40 - 1;
       arg = arg * arg;
       ez(si,sj) = (1 - 2*arg)*exp(-arg);    
    end
    
   
    
    
end
     
for i=1:2:sizeX
    for j=1:2:sizeY
        ezyee(fix(i/2+1),fix(j/2+1)) = ez(i,j);
    end
end

for i=1:2:sizeX
    for j=2:2:sizeY
        hxyee(fix(i/2+1),fix(j/2)) = hx(i,j);
    end
end

for i=2:2:sizeX
    for j=1:2:sizeY
        hyyee(fix(i/2),fix(j/2+1)) = hy(i,j);
    end
end
hold off
figure(1)
surf(ezyee)

xlabel('x')
ylabel('y')
zlabel('z')