function HW2P3F1()
% solving a periodic domain with IC between 0 period
global period;
period = 1;

c = 1.0;
rho = 1;
E = c * c * rho;
Z = sqrt(E * rho);
kn = 0.50;
numSegments = 100;
finalTime = period / c * 0.35;

numSegs = [10, 100, 1000];
kns = [0.10, 0.49, 0.50, 0.51, 0.99, 1.01, 1.99, 2.01];

ns = length(numSegs);
knsz = length(kns);
for nsi = 1:ns
    numSegments = numSegs(nsi);
    for kni = 1:knsz
        kn = kns(kni);
        FDFV_comparisonED(numSegments, finalTime, kn, c, E, Z);
    end
end
close('all');


function FDFV_comparisonED(numSegments, finalTime, kn, c, E, Z)
global period;
addRiemannCorrection = 0;
[x, ysFD, numberOfTimeStep, finalTimeComputed, h] = FVED(c, E, Z, kn, numSegments, finalTime, addRiemannCorrection);
yFinalFD = ysFD(numberOfTimeStep + 1, :);
addRiemannCorrection = 1;
[x, ysFV, numberOfTimeStep, finalTimeComputed, h] = FVED(c, E, Z, kn, numSegments, finalTime, addRiemannCorrection);
yFinalFV = ysFV(numberOfTimeStep + 1, :);
xExact = 0:(period/1000):period;
yFinalExact = exactSolution(xExact, finalTimeComputed, c, E);

base{1} = 'HW2_FVFD_comparison_1Field';
[status,message,messageid] = mkdir(base{1});
base{2} = [base{1}, '/FD'];
[status,message,messageid] = mkdir(base{2});
base{3} = [base{1}, '/FV'];
[status,message,messageid] = mkdir(base{3});
leg{1} = {'FD', 'FV', 'exact'};
leg{2} = {'FD', 'exact'};
leg{3} = {'FV', 'exact'};

for i = 1:3
    base{i} = [base{i}, '/F1-'];
end

figure(1);
clf
plot(x, yFinalFD, '-b', 'Linewidth', 2);
hold on;
plot(x, yFinalFV, '-r', 'Linewidth', 2);
hold on;
plot(xExact, yFinalExact, '-k', 'Linewidth', 2);
hold on;

figure(2);
clf
plot(x, yFinalFD, '-b', 'Linewidth', 2);
hold on;
plot(xExact, yFinalExact, '-k', 'Linewidth', 2);
hold on;
figure(3);
clf
plot(x, yFinalFV, '-r', 'Linewidth', 2);
hold on;
plot(xExact, yFinalExact, '-k', 'Linewidth', 2);


for pl = 1:3
    figure(pl);
    legend(leg{pl}, 'FontSize', 18, 'Location', 'North');
    legend('boxoff');

    titleText = ['h = ', num2str(h), 'kn = ', num2str(kn), ' T = ', num2str(finalTimeComputed), ' c = ', num2str(c)];
    title(titleText, 'FontSize', 18);
    xlabel('x', 'FontSize', 18);
    ylabel('u', 'FontSize', 18);
    printName = [changeName4Output([base{pl}, titleText]), '.png'];
    print('-dpng', printName);
end
close('all')

function [x, ys, numberOfTimeStep, finalTimeComputed, h] = FVED(c, E, Z, kn, numSegments, finalTime, addRiemannCorrection)
global period;
h = period / numSegments;
k = kn * h / c;

numberOfTimeStep = ceil(finalTime / k);
finalTimeComputed = numberOfTimeStep * k;
ys = zeros(numberOfTimeStep + 1, numSegments);

for m = 1:numSegments
    x(m) = (m - 1) * h;
    [ys(1, m), ydot(m)] = IC(x(m), E);
    ym1(m) = ys(1, m) - k * ydot(m);
end


for n = 1:numberOfTimeStep
    for m = 1:numSegments
        mm1 = m - 1;
        mm2 = m - 2;
        mp1 = m + 1;
        mp2 = m + 2;
        if (m == 1)
            mm1 = numSegments; mm2 = mm1 - 1;
        elseif (m == 2)
            mm2 = numSegments;
        elseif (m == numSegments)
            mp1 = 1;    mp2 = 2;
        elseif (m == numSegments - 1)
            mp2 = 1;
        end
        if (n == 1)
            y_mm1_nm1 = ym1(mm1);
            y_m_nm1 = ym1(m);
            y_mp1_nm1 = ym1(mp1);
        else
            y_mm1_nm1 = ys(n - 1, mm1);
            y_m_nm1 = ys(n - 1, m);
            y_mp1_nm1 = ys(n - 1, mp1);
        end
        ys(n + 1, m) = 2 * ys(n, m) - y_m_nm1 + kn * kn / 4.0 * (ys(n, mp2) + ys(n, mm2) - 2 * ys(n, m));
        if (addRiemannCorrection == 1)
            ys(n + 1, m) = ys(n + 1, m) + kn * ((ys(n, mp1) + ys(n, mm1) - 2 * ys(n, m)) - (y_mp1_nm1 + y_mm1_nm1 - 2 * y_m_nm1));
        end
    end
end
ys(:, numSegments + 1) = ys(:, 1);
x(numSegments + 1) = period;



function [y, ydot, yStrs] = IC(x, E)
global period;
xp = mod(x, period);

a = (.5 - 0.1) * period;
b = 0.5 * period;
c = (0.5 + 0.1) * period;
if (xp < a)
    y = 0.0;
    yStrs = 0.0;
elseif (xp < b)
    xn = (xp - a) / (b - a);
    yStrs = E / (b - a);
    y = xn;
elseif (xp < c)
    xn = (xp - b) / (c - b);
    y = 1 - xn;
    yStrs = E / (b - c);
else
    y = 0.0;
    yStrs = 0.0;
end
ydot = 0.0;


% assuming initial velocity is zero
function y = exactSolution(x, time, c, E)
global period;
for m = 1:length(x)
    y(m) = 0.5 * (IC(x(m) - c * time, E) + IC(x(m) + c * time, E));
end
