clear;

%% Data

rng(2022);

T = 4;
nSteps = 10.^(2);
nTrajectories_vec = [1];


%% Brownian motion

dt = T/nSteps;
disp(dt)
tVec = 0:dt:T;
for j = 1:length(nTrajectories_vec)       
    nTrajectories = nTrajectories_vec(j); 
    W = BM(T, nSteps, nTrajectories);
    figure;
    plot(tVec, W, 'color', 0.75*ones(1, 3));
    xlabel('t', 'Fontsize', 20);
    title(['{\Delta}t = ', num2str(dt), ', ', num2str(nTrajectories), ' sample paths'], 'Fontsize', 20);
    hold on;
    plot(tVec, mean(W), 'r', 'Linewidth', 2);
    plot(tVec, zeros(size(tVec)), 'g', 'Linewidth', 2);
    set(gca, 'FontSize', 20);
    saveas(gcf, ['brownian_motion_14/figures', num2str(j), '.png']);
end

   
%% Refinement

alpha_vec = [0.25,0.5,0.75];
for k = 1:length(alpha_vec)
    alpha = alpha_vec(k);
    dt = T/nSteps;
    nSteps_r = nSteps*2;
    dt_r = T/nSteps_r;
    tVec_r = 0:dt_r:T;
    W_r = [zeros(nTrajectories_vec, 1),];
    
    for j = 1:length(nTrajectories_vec)       
        nTrajectories = nTrajectories_vec(j); 
        Y = sqrt(dt*(alpha-alpha^2))*randn(nTrajectories, (nSteps_r-nSteps));
        for t =1:nSteps_r
            if mod(t, 2) == 0
                disp("yes");
                W_r(end+1) = W((t/2)+1);
            else
                disp(floor(t/2)+1);
                disp(ceil(t/2)+1);
                W_r(end+1) = (1-alpha)*W(floor(t/2)+1) + alpha*W(ceil(t/2)+1) + Y(floor(t/2)+1);
            end
        end
    
        figure;
        xlabel('t', 'Fontsize', 20);
        title(['{\Delta}t = ', num2str(dt), ', alpha=', num2str(alpha)], 'Fontsize', 20);
        hold on;
        plot(tVec, W, 'r', 'Linewidth', 2);
        plot(tVec, W, '*');
        plot(tVec_r, W_r, 'g', 'Linewidth', 2);
        plot(tVec_r, W_r, '.');
        set(gca, 'FontSize', 20);
        saveas(gcf, ['brownian_motion_14/figures_alpha', num2str(alpha), '.png']);
    end
end