%% Matlab version of typical Mathcad calculations
%% Basic calculations - solution to quadratic equation:
%%  {\mathrm{ax}}^2 +\mathrm{bx}+c=0

a=1; b=2; c=3;

% x_{1\mathrm{st}} =\frac{-b+\sqrt{b^2 -4\mathrm{ac}}}{2a}        x_{2\mathrm{nd}} =\frac{-b-\sqrt{b^2 -4\text{ac}}}{2a}

x_1st = (-b + sqrt(b^2 - 4*a*c)) / (2*a);
x_2nd = (-b - sqrt(b^2 - 4*a*c)) / (2*a);
disp (['x_1st = ' num2str(x_1st) '    x_2nd =' num2str(x_2nd)]);
% checking results in my_quadratic function:

% function [f] = my_quadratic(x, a, b, c)
% % Function to evaluate the quadratic function with predefined a, b, c
% f = a*x.^2 + b*x + c;
% end

fprintf ('f(x_1st) = %i     f(x_2nd) = %i\n', my_quadratic(x_1st, a, b, c), my_quadratic(x_2nd, a, b, c));

%
%% Plotting a function with automated ranges and number of points

ezplot('1*x^2 + 2*x + 3')

%
%% Plotting a function using a vector of values, with custom display

x = -5 : 0.05 : 3;
y = my_quadratic(x, a, b, c);
plot (x, y)
title ('Custom plot of quadratic function');
xlabel('x');
ylabel('f(x)');
grid on
% Display both ends of x vector
x_length = length(x);
for i = 1:15   
    x_lower(i) = x(i);
    x_upper(i) = x(x_length - 15 + i);
end
x_lower
x_upper

%
%% Using units and formatted display
%   (unit conversion functions available in Aerospace Toolbox only)

% m = convmass (100, 'lbm', 'kg');
m = 100 / 2.204622622;    % conver lbm to kg
% v = convvel (60, 'mph', 'm/s');
v = 60 * 0.44704;     % convert mph to mps
% a = convacc (20, 'ft/s^2', 'm/s^2');
a = 20 * 0.3048;     % convert fps2 to mps2
p = m*v
F = m*a;
% convforce(F, 'N', 'lbf')
F = F / 4.448       % conver N to lbf

%
%% Symbolic algebra
% \frac{x}{\left(2x-3\mathrm{xy}\right)}=\frac{{\left(x-2\right)}^2 }{y+2}
% solve equation:

syms x y
eqn = x / (2*x - 3*x*y) == (x-2)^2/(y+2);
x_ans = solve (eqn)

% evaluate solutions:

subs(x_ans, 'y', 5)

%
%% Symbolic calculus

syms x a
fx = (x - a)^2 + 10*sin(2*x)/x
dfx = diff (fx)
x = -3:0.1:5;
a = 6.096;
y = eval(fx);
dy = eval(dfx);
plot (x, y)
title ('f(x)')
plot (x, dy)
title ('df(x)')

%
%% Vector and matrix calculations
% v=\left\lbrack \begin{array}{c}
% -1\\
% -2
% \end{array}\right\rbrack

vx = -1; vy = -2;
v = [vx; vy];
v = vx + j*vy

% \left|v\right|

abs(v)

% v\cdot v

dot(v,v)
v_ang = angle(v)*180/pi;
display (['angle of v = ' num2str(v_ang) ' deg']);
display (['polar form of v = ' num2str(v_mag) ' < ' num2str(v_ang)]);

% A=\left\lbrack \begin{array}{c}
% 1 & 2 & 3\\
% 2 & 1 & 5\\
% 0 & -2 & 3
% \end{array}\right\rbrack


A = [1 2 3; 2 1 5; 0 -2 3];

% A^{-1}

inv(A)

% I={\mathrm{AA}}^{-1}

I = A * A_inv

%
%% Programming a piecewise function
% f\left(x\right)=\left\lbrace \begin{array}{c}
% x & x<1\\
% -{\left(x-1\right)}^2 +1 & 1\le x\le 3\\
% -3 & x>3
% \end{array}\right.

% function [f] = my_piece_wise(x, a, b, c)
% % Function to evaluate the quadratic function with predefined a, b, c
% if (x < 1)
%     f = x;
% elseif ((x >= 1) && (x <= 3))
%     f = -(x-1)^2 + 1
% else
%     f = -3
% end

x = -2 : 0.1 : 5;
clear y;
for (i = 1 : length(x))
    y(i) = my_piece_wise(x(i));
end
plot (x, y);
title ('Plot of piece-wise function');
xlabel('x');
ylabel('f(x)');
axis([-2.5 5.5 -3.5 1.5]);

%
%% General programming problem example
% Find the sum of the first N numbers divisible by 3

% function [i total] = my_program(N)
% % Function to calculate the sum of the first N numbers divisible by 3
% i = 0;
% n = 0;
% total = 0;
% 
% while (n < N)
%     i = i + 1;
%     remainder = mod (i, 3);
%     if (remainder == 0)
%         total = total + i;
%         n = n + 1;
%     end
% end

N = 10000;
[i, total] = my_program(N)

%
%% Finding roots
% f\left(x\right)=2x^2 -4\mathrm{sin}\left(x\right)-2

% function [f] = my_root_func(x)
% % Function to evaluate the quadratic function with predefined a, b, c
% f = 2*x^2 - 4*sin(x) - 2;
% end

% roots for different guesses:

x0 = 1;
fzero (@my_root_func, x0)
x0 = -1;
fzero (@my_root_func, x0)
x = -1 : 0.1 : 2;
y = 2*x.^2 - 4*sin(x) - 2;
plot (x, y);
title ('Plot of root function');
xlabel('x');
ylabel('f(x)');
hold on;
x = [-1 2];
y = [0 0];
plot (x,y,'LineStyle',':','Color',[1 0 0]);
axis([-1.5 2.5 -4 4]);
hold off;

%
%% Solving a set of nonlinear equations
% x=2-y^2
% y=\frac{\mathrm{sin}\left(x\right)}{x}+\mathrm{xy}

% function [ F ] = my_non_linear_equations( x )
% % Define set of nonlinear equations to be solved numerically
% F = [x(1) - 2 + x(2)^2;  x(2) - sin(x(1))/x(1) + x(1)*x(2)];
% end

x0 = [1; 1];  % initial guesses
[x_sol,fval] = fsolve(@my_non_linear_equations,x0);
x_sol

% Checking results (solving symbolically and plotting)
syms x y;
fa = solve(x == 2 - y^2, y)
fb = solve(y == sin(x)/x + x*y, y)
x = x_sol(1)
eval(fa(1))
eval(fb)
x = 0.01 : 0.04 : 0.5;
ya = eval (fa(1));
yb = eval (fb);
plot (x, ya);
hold on;
plot (x, yb, 'Color',[1 0 0]);
legend ('fa(x)', 'fb(x)');
hold off;

%
%% Iterative calculations

clear x y;
x(1) = 1, y(1) = 1
for i = 1 : 6
    x(i+1) = x(i) + 2;
    y(i+1) = (x(i) + x(i+1)) / 2;
end
x
y

%
%% Finding an optimal solution given constraints
% (requires Optimization Toolbox)
% \mathrm{maximize}
% z\left(x,y\right)={\left(x-1\right)}^2 -\mathrm{xsin}\left(y\right)
% \mathrm{subject}\text{ }\mathrm{to}
% x>-2
% x<2y^2 +3
% -3<y<5

% function [ F ] = my_objfun( x )
% % Function definition for constrained optimization problem
% % (minus sign in front for max vs. min)
% F = - ((x(1)-1)^2 - x(1)*sin(x(2)));
% end

% function [c, ceq] = my_confun(x)
% % Nonlinear inequality constraints
% c = [-x(1) - 2; x(1) - 2*x(2)^2 - 3; x(2) - 5; -x(2) - 3];
% % Nonlinear equality constraints
% ceq = [];

x0 = [1; 1];
[x,fval] = fmincon(@my_objfun,x0,[],[],[],[],[],[],@my_confun);
x
-fval   % minus for max vs. min