% Schneider book pg. 196 first 2D sim
clear
%sizeX = 201;  
%sizeY = 161;
%maxTime = 250;
lc=1;

sizeX = 120;  
sizeY = 60;
maxTime = 150;
%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);

eznm1=zeros(sizeX, sizeY);ezn1=zeros(sizeX, sizeY);
hxnm1=zeros(sizeX, sizeY);hxn1=zeros(sizeX, sizeY);
hynm1=zeros(sizeX, sizeY);hyn1=zeros(sizeX, sizeY);

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

CE = imp0/sqrt(2);
CH = 1/(sqrt(2)*imp0);

figure(1)
        

for n=1:maxTime  
    

    eznm1=ezn; ezn=ez; hynm1=hyn; hyn=hy; hxnm1=hxn; hxn=hx;
    
    for i=2:sizeX-1
        for j=2:sizeY-1
            hy(i,j) = hynm1(i,j) + CH*(ezn(i+1,j) - ezn(i-1,j));
            hx(i,j) = hxnm1(i,j) - CH*(ezn(i,j+1) - ezn(i,j-1));
        end
    end
    for i=2:sizeX-1
        for j=2:sizeY-1         
            ez(i,j) = eznm1(i,j)+CE*(hyn(i+1,j)-hyn(i-1,j)- ... 
            (hxn(i,j+1)-hxn(i,j-1))); 
        end
    end
       
    % apply boundary conditions here along top and bottom
       
       for i=2:sizeX-1
        ez(i,1) = ez(i,2);
        hx(i,1) = hx(i,2);
        hy(i,1) = hy(i,2);
        
        ez(i,sizeY) = ez(i,sizeY-1);
        hx(i,sizeY) = hx(i,sizeY-1);
        hy(i,sizeY) = hy(i,sizeY-1);
       end
       
       % apply boundary condition for end
       for j=1:sizeY
           %ez(sizeX,j) = ez(sizeX-1,j);  %PMC
           hx(sizeX,j) = hx(sizeX-1,j);  %PEC
           
           hy(sizeX,j) = hy(sizeX-1,j);
           
       end
    
       % put a PEC vertical plate 41 cells from the left edge and
       % 60 cells high
    
       if 1 == 1
        for j= 21:40; %20:40
            ez(80, j) = 0;
            %hx(20, j) = 0;
            %hy(20, j) = 0;
             hx(80, j) = hx(79, j);
             hy(80, j) = hy(79, j);
             ez(81, j) = 0;
             hx(81, j) = hx(82, j);
             hy(81, j) = hy(82, j);
        end
       end
 
       % generate wave 
       arg = pi * Cdtds*n / 40 - 1;
       arg = arg * arg;
       for j= 1:sizeY
           ez(1,j) = (1 - 2*arg)*exp(-arg);
           hy(1,j) = -ez(1,j) / imp0;
           
       end
    
    
     lc=lc+1;
     if mod(lc,10)==0     % animate graph
        
        surf(ez')   
        shg

    end
    
   
    
    
end

figure(2)
surf(ez')
%ezyee
%hxyee
%hyyee

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