%% Solutions to HW1 ENGN0040 Dynamics and Vibrations, Spring 2017
% 12 POINTS TOTAL

clear all;

%
% Problem #1 [2 POINTS]

syms x y
eq1 = x^2 + 2 * y^2 == sym(18);
eq2 = x + y - x * y == sym(1);
[x, y] = solve([eq1, eq2], [x, y])

%
% Problem #2 [2 POINTS]

syms P(x,s) s
P(x, s) = 1/(x * (1 - x)) * exp(-(log(x/(1 - x))/s)^2)

%
%

fplot(P(x,0.5), [0.01, 0.99],'DisplayName','P(x,0.5)','Color',[1 0 0],...
     'LineWidth',2);
hold on
fplot(P(x,1), [0.01, 0.99],'DisplayName','P(x,1)','Color',[0.5 0.5 0],...
     'LineWidth',2);
fplot(P(x,2), [0.01, 0.99],'DisplayName','P(x,2)','Color',[0 1 0],...
     'LineWidth',2);
fplot(P(x,2.5), [0.01, 0.99],'DisplayName','P(x,2.5)','Color',[0 0 1],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x'},'FontAngle','italic','FontName','Times New Roman');
title('The Logistic Normal Distribution');
ylabel({'P(x,m)'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

%
% Problem #3 [1 POINT]

syms s
assume(0 < s);
int(P(x,s), x, sym(0), sym(1))

%
% Problem #4 [2 POINTS]

syms a
f = x^2/(x^3 + a^3);
xatmax = solve(diff(f, x) == sym(0), x)

%
% (this is all the roots. The first cant be a max (f=0); the third and fourth are complex, so we want the second one)

fmax = simplify(subs(f, x, xatmax(2)))

%
% Problem #5:

clear all
syms R Ts g lambda ps
syms p(z)
diffeq = diff(p(z), z) == (g * p(z))/(R * (Ts - lambda * z))
IC = p(sym(0)) == ps;
press = simplify(dsolve(diffeq, IC,symvar(p(z))))
 

% [2 POINTS]
% To simplify:
% 1. Combine the logs: p = ps exp( (g/lambda R) log (Ts/(Ts-lambda z))
% 2. Put the g/lambda R inside the log: p = ps exp( log ( (Ts/(Ts-lambda z)) ^ (g/lambda R) )
% 3. Get rid of the exp(log) p = ps (Ts/(Ts-lambda z)) ^ (g/lambda R)
% 4. Alternatively p = ps ( 1 - lambda z/Ts ) ^(-g/lambda R)
% [1 POINT]
% Calculating the % change in pressure over 160 feet:

preshi = vpa(subs(press, {R, Ts, ps, lambda, g, z}, {287,298,101325,6, 9.81, 3072.384/1000}));

preslo = vpa(subs(press, {R, Ts, ps, lambda, g, z}, {287,298,101325,6,9.81,3023.616/1000}));

percenterr = (100 * (preshi - preslo))/preshi

%
% [2 POINTS]
%


