%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% batch.m
% Batch execution script. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function z = batch(n) % Remember to increase n by powers of 2

fprintf(1,'\n\nStarting...\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 0 -  Finding the "correct" answer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nFinding the correct answer...\n');
y=[pi/2 0]; y2 = y; t=0;

n2 = 100000; h2 = 10/n2;
for i = 1:n2
    t = t + h2;
    y = trapstep(t, y, h2);
end
correct_ans = y(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 1 - Characterize the function E(h)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nCharacterizing E(h)...\n');
q = 0;
while(n > 10)
    q = q + 1; h = 10 / n;
    
    y = [pi/2 0]; y2 = y; t = 0;

    % Trap Method / Euler Method
    for i = 1:n
        t = t + h;
    
        y = trapstep(t, y, h);
        y2 = eulerstep(t, y2, h);
    end
    trap_ans = y(1)
    euler_ans = y2(1)

    trap_error(q) = correct_ans - trap_ans;
    euler_error(q) = correct_ans - euler_ans;
    
    n = n * .5;
end

trap_error'
euler_error'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 2 - Damped Pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nAnimating the damped pendulum...\n');
d = 0.1; n = 100; h = 10 / n;
while ( d <= 5.1 )
    y = [pi/2 0]; t = 0;

    d
    damppend([0 10], y, h, 1, d);
    d = d + 1.0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 3A - Forced Pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nAnimating the forced pendulum...\n');
d = 0.5; a = 60; n = 100; h = 10 / n;
while (a <= 150)
    y = [3.1 0]; t = 0;

    a
    forcepend([0 10], y, h, 1, d, a);
    a = a + 30;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 3B - Forced Pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nAnimating the forced pendulum (again)...\n');
d = 0.5; a = 150; n = 100; h = 10 / n; 
yi = 3.1;
while ( yi >= 2.7 )
    y = [yi 0]; t = 0;

    yi
    forcepend([0 10], y, h, 1, d, a);
    yi = yi - 0.1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 4 - Double Pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\nAnimating the double pendulum...\n');
n = 1000; h = 10 / n; y = [pi/2 0 pi/2 0];
friction = 0.1;
doublepend([0 n], y, h, 1, friction);

fprintf(1,'\n\nFinished.\n');