%*************************************************************************
%FDTD-1D
%Editer: Chun-Feng Lai
%*************************************************************************
%description: soft source, 1D, absorbing boundary condition (which really
%only work for epsrel=murel=1), reflection at an interface
%type ward_fdtd1d to run. 
clear
%things to specify before running
w=10.0; %(THz) driving frequency
phiE=0*pi; %phase of the E field source
ie=600.0; %(problem space size in microns) note-speed of light in vacuum is 300 microns/ps
sx=100; %(microns) source location 
st=0.5; %(ps) source temporal duration
T=3.0; %(picoseconds) duration of the sim
upnum=20; %display output every upnum timesteps
startd=250; %(meshpoint) where the dielectric starts (epsrel1 on left, epsrel2 on right)
epsrel1=1; 
epsrel2=10;
murel1=1;
murel2=1;


%Past here, change at your own risk***************************************
%constants and convert to simulation units (SI)
w=2*pi*w*1e12; %adjusting to radial units and SI units
T=T*1e-12;
epso=8.854187817620389850536563e-12;
muo=1.2566370614359172953850573e-6;
c0=1/sqrt(epso*muo);
if epsrel1>epsrel2
    lambda=2*pi*c0/(epsrel2^2*w); %wavelength of excitation
else 
    lambda=2*pi*c0/(epsrel1^2*w); %wavelength of excitation
end
dx=lambda/20; % lambda over 20 is a good rule of thumb
ie=round((ie*1e-6)/dx); %convert proble space to discrete units
dt=dx/(2*c0); %Courant Condition for stability, which is related to Nyquist theorem
nsteps=round(T/dt);
sx=round((sx*1e-6)/dx);
st=(st*1e-12)/dt;

%initialize fields

Ex(1:ie+1)=0.0;
Hy(1:ie)=0.0;

%initialize boundary conditions
Ex_lm2=0;
Ex_lm1=0;
Ex_hm2=0;
Ex_hm1=0;

%make dielectric array
eps(1:ie+1)=epso;
mu(1:ie+1)=muo;
eps(1:startd)=eps(1:startd)*epsrel1;
mu(1:startd)=mu(1:startd)*murel1;
eps(startd+1:end)=eps(startd+1:end)*epsrel2;
mu(startd+1:end)=mu(startd+1:end)*murel2;

%initialize graphics
x=linspace(dx,ie*dx*1e6,ie);
subplot(2,1,1),plot(x,Ex(1:ie)),xlabel('Position (microns)'),ylabel('Ex');
subplot(2,1,2),plot(x,Hy),xlabel('Position (microns)'),ylabel('Hy');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nsteps/upnum,gcf,rect);

%start simulation
t=0;
for n=1:nsteps
    %source for E
    if n<st
        sourceE=sin(w*t+phiE);
        Ex(sx)=Ex(sx)+sourceE;
    end
    %update time by half time step
    t=t+dt/2;
    
    %update E field***************
    for i=2:ie
        Ex(i)=Ex(i)+(dt/(eps(i)*dx))*(Hy(i-1)-Hy(i));
    end
    %boundary condition*******************************************
    Ex(1)=Ex_lm2;
    Ex_lm2=Ex_lm1;
    Ex_lm1=Ex(2);
    Ex(ie-1)=Ex_hm2;
    Ex_hm2=Ex_hm1;
    Ex_hm1=Ex(ie-2);
    %boundary condition*******************************************
    %update time by half time step
    t=t+dt/2;
    
    %update H field***************
    for i=1:ie
        Hy(i)=Hy(i)+(dt/(dx*mu(i)))*(Ex(i)-Ex(i+1));
    end

    %update graph every upnum time steps
    if mod(n,upnum)==0
        subplot(2,1,1),plot(x,Ex(1:ie)),xlabel('Position (microns)'),ylabel('Ex');
        subplot(2,1,2),plot(x,Hy),xlabel('Position (microns)'),ylabel('Hy');        
        M(:,n/upnum)=getframe(gcf,rect);
    end
    
end