%Instructors solution to HW1_16, part 2
%Bare-bones...

close all;

% choose time step:
dtoverdx = 0.5;

%define anlytic target
Disp = @(k) 2/dtoverdx*asin(dtoverdx*sin(k/2));

% define initial electric field fucntion
InitialE = @(x,k) sin(k*x);

% computational domain
NX = 2048;
cx = linspace(0,NX-1,NX);

% parameters
nl = NX/2;
test=0;

%scanning wavenumbers: nl will stand for inverse number of points per wavelength
for nl=-NX/2:50:NX/2

%avoid zero wavenumber
if abs(nl) > 0

% this is the selected wavelength
la = NX/nl;

% and this is the selected wavenumber
k0 = 2.0*pi/la;

% set an initial condition with a sharp value of wavenumber
EF = +InitialE(cx,k0);
HF = -InitialE(cx - 0.5 - 0.5*dtoverdx,k0);

%num of simulation steps
%multiply this by 8 to obtain smaller error
NS = 1024;

% auxiliary = holder for field history at a selected point
history = zeros(NS,1);

% execute some steps
for s=1:NS
EF = EF + dtoverdx*(circshift(HF,[1,-1]) - HF);
HF = HF - dtoverdx*(circshift(EF,[1,+1]) - EF);

%collet data during simulation
history(s,1) = EF(NX/2);
end

%take spectrum at the selected point
hspectrum=abs(fft(history));

%find peaks
peak = find(hspectrum>=max(hspectrum));

%convert peaks to omegas (angular frequency)
oplu = (peak(1)-1)*2*pi/(NS*dtoverdx);
omin = -(NS+1-peak(2))*2*pi/(NS*dtoverdx);

display([k0 oplu omin])

%while scanning k0, build lists to plot results later
if test==0
wavenumbers= k0;
omegaspos = oplu;
omegasneg = omin;
test=1;
else
wavenumbers=cat(1,wavenumbers, k0);
omegaspos =cat(1,omegaspos, oplu);
omegasneg =cat(1,omegasneg, omin);
end

end

end

figure(1)
plot(wavenumbers,omegaspos);
hold on;
plot(wavenumbers,omegasneg);
target = Disp(wavenumbers);
title('Four Branches of Dispersion Curve');

figure(2)
plot(wavenumbers,abs(target) - omegaspos,'r');
title('Dieprsion Measurement Error')