clear;
%Boundary Conditions
%bc=0 for open, bc=1 for reflecting
bc=1;
%Permeability and Permittivity
mu=pi*4e-7;
epsilon=8.854e-12;
%Calculate Impedance
eta=sqrt(mu/epsilon)
%Calculate Wave Speed C
C=sqrt(1/(mu*epsilon))
%Spatial step size and temporal step duration
dx=0.001875;  %Grid step distance in meters
%Set CFL Number
CFL=1;                   %Any more than one -> Unstable!
%Calculate Time Step dt
dt=CFL.*dx./C            
%Max domain size in space, and evolution duration time
Imin=0;   %First Grid Point
Imax=99;  %Last Grid Point
Xmin=(Imin-1)*dx;  %One Grid Point added for Boundary Conditions
Xmax=(Imax+1)*dx;  %Domain spans from -1 to 101 meters
Xdomain=(Imax-Imin+1)*dx;

%Number of waves that fit in box (determines wavelength)
%Half-integer multiples preferred, others accepted for fun examples
M=1.5;
%Calculate Resulting Wavelength
lambda=Xdomain/M
%Calculate Resulting Frequency
frequency=C/lambda
%Number of Full-Domain Propagations Per Wavelength to Run
NT=10;
%Wave Amplitude (H-field in [A/m])
A=1;
%Gaussian envelope "width" (standard dev.) as fraction of total runtime
G=1/5;
%Gaussian packet "peak" time as fraction of total runtime
%Note that 0 will result in a packet that begins at full-strength and
%decays over time
T=1/2;


Nmin=0;                   %First time step to take
Nmax=Xdomain*NT/(C*dt);   %Last time step to take
Tmin=Nmin*dt;             %Time of first time step
Tmax=Nmax*dt;             %Time of last time step

%Domain X-axis Reference
x=[Xmin:dx:Xmax];

%Calculate later times t...
for n=Nmin:1:Nmax   
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        %Define initial zero state
        Hya=x.*0;
        Exa=x.*0;
        Hy=Hya; %Lax-Wendroff Init.
        Ex=Exa;       
        %Initialize temporary U(n+1) 
        H1ya=Hy;
        E1xa=Ex;
    else        %Evolve solution at later timesteps
               
        for im=Imin-1:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1/2) for Lax-Wendroff 2-Step Method
            H1ya(i)=(1/2)*(Hy(i+1)+Hy(i))-((1/mu)*Ex(i+1)-(1/mu)*Ex(i))*dt/(2*dx);       
            E1xa(i)=(1/2)*(Ex(i+1)+Ex(i))-((1/epsilon)*Hy(i+1)-(1/epsilon)*Hy(i))*dt/(2*dx);     
        end;

        %Perform second Lax-Wendroff step:
        %Update results for n+1/2 (Lax-Wendroff)
        Hya=H1ya; 
        Exa=E1xa; 
        %Initialize Array
        H1ya=Hya;
        E1xa=Exa;
        
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1) for Lax-Wendroff 2-Step Method
            H1y(i)=Hy(i)-((1/mu)*(Exa(i)-Exa(i-1)))*(dt/dx);   
            E1x(i)=Ex(i)-((1/epsilon)*(Hya(i)-Hya(i-1)))*(dt/dx);            
        end;
        
        %Set solution to updated values
        Hy=H1y;
        Ex=E1x;
       
        %Update boundaries for Lax-Wendroff
        Hy(Imin+1)=Hy(Imin+2)+A*sin(2*pi*((n*dt)/(Xdomain/(C*M))))*exp(-.5*((n*dt-Tmax*T)^2)/(Tmax*G)^2);
        Hy(Imax+3)=Hy(Imax+2);
        Ex(Imin+1)=(1-2*(bc==1))*Ex(Imin+2);
        Ex(Imax+3)=(1-2*(bc==1))*Ex(Imax+2);
    end;
    
    %Show Results    
    figure(1);
    subplot(2,1,1);
    plot(x,Hy,'r-',x,Ex/eta,'b-');
    axis([Xmin Xmax -3*A 3*A]);
    xlabel('Horizontal Position (m)');
    ylabel('Amplitude (A/m)');
    title(['$E_{x}/\eta$ (Blue) and $H_{y}$ (Red) @ t=',num2str(n*dt), ', n=',num2str(n),''],'Interpreter','latex');
    pause(0.005); 
    subplot(2,1,2);
    plot(x,Ex.*Hy,'g-');
    axis([Xmin Xmax -A*eta A*eta]);
    xlabel('Horizontal Position (m)');
    ylabel('Amplitude (W/m^2)');
    title(['$S_{z}$ @ t=',num2str(n*dt), ', n=',num2str(n),''],'Interpreter','latex');
    pause(0.005); 
    
end;
