% script for excercise 2 of case study 1 for EMATH171
% in this script we solve for power using newtons method repetating untill we
% have a tolerence of 0.0001
clear
clc
close all

% Coefficients for the function and its deriviative

m = 1500; % mass
g = 9.81; % gravity
R = 0.205; % wheel raduis
beta = 0.44;
acell_angular = 420;
rgb = 0.8;
rfd = 3.5;
drag_coff = 0.3;
area = 2;
fluid_density = 1.2;

crr = 0.01;

acceleration = 0;
theta = 10;

% functions
k = R /(rgb * rfd);

pwr = @(v) ((acell_angular .* v)/k) - ((beta .* v.^2) / (k.^2));
% power output

roll = @(v) crr .* m .* g .* (cos(theta)) .* v;
% power dissipated by rolling resistance

drag = @(v) drag_coff .* area .* 0.5 .* fluid_density * v.^3;
% power dissipated by aerodynamic drag

acell = @(v) m .* acceleration .* v;
% power required for acceleration at rate of dv/dt

grad = @(v) m .* g .* (sin(theta)) .* v;
% Gradient climbing resistance

f = @(v) drag(v) + grad(v) + roll(v) + acell(v) - pwr(v);
% origional function

% derrivative functions
drag_derivative = @(v) 3 .* drag_coff .* area .* 0.5 .* fluid_density .* v^2;

grad_derivative = m * g * sin(theta);

roll_derivative = crr * m * g * cos(theta);

pwr_derivative = @(v) (acell_angular / k) - ((2 .* beta .* v)/(k.^2));

g = @(v) drag_derivative(v) + grad_derivative + roll_derivative + m * acceleration - pwr_derivative(v);
% deriviative

h = @(old_x) old_x - ((f(old_x))/(g(old_x)));
% newtons method propper

% Coefficients for newtons method

N = 200;
% maximuim amount of repetitions, unlikley to actualy reach this value

x = 200;
% initial guess

tol = 0.0001;
%tolerence of diffrence between iterations

out_array = [x];
% output array of the results of newtons methods approximations

error_array = [abs(f(x))];
% array of the error

delta_array = [];
% array of change between iterations

counter = 0;
% tracks number of iterations

% Newtons method implementaion

for item = 1:N
    old_x = x; % saving old x
    delta_array(item) = (x);
    x = h(x); % creating new x
    counter = counter + 1;
    out_array(item+1) = x;
    error_array(item+1) = abs(f(x));
    delta_array(item) = abs(delta_array(item)- x);
    if and((abs(f(x))<tol),(abs(x-old_x)<tol))
        counter = counter - 1;
        break
    end
end

% display
figure
hold on
x_actual_array = [-200:1:1500]; % road speed v
plot(x_actual_array, pwr(x_actual_array))
plot(x_actual_array, roll(x_actual_array))
plot(x_actual_array, drag(x_actual_array))
plot(x_actual_array, acell(x_actual_array))
plot(x_actual_array, grad(x_actual_array))
plot(x_actual_array, f(x_actual_array))
legend('power', 'rolling resistance', 'drag resistance', 'aceleration', 'gradient resistance', 'total')
hold off
x
counter
out_array;
error_array;
delta_array;

%
% The xlabel and ylabel commands generate labels along x-axis and y-axis.
% The title command allows you to put a title on the graph.
% The grid on command allows you to put the grid lines on the graph.
% The axis equal command allows generating the plot with the same scale factors and the spaces on both axes.
% The axis square command generates a square plot.