%% ENGN0040 2021 Homework 4 Solutions
%
%% Problem 1
%

syms r real
syms U(r)
Jtoev = 6.2418e18; % Conversion from Joules to eV
nmtoA = 10;        % Conversion from nanometers to Angstroms
Atom = 1e-10;       % Conversion from Angstroms to m
e = 1.6022e-19;     % Electron charge, Coulombs
e0 = 8.8542e-12;    % Permittivity of free space, Farads/m
B = 3.38e-20*exp(2.96e09*(1.463+1.585)*10^(-10))*Jtoev; % B, in eV - this is the 8.33e-20 given
C = 48e-25*Jtoev*(nmtoA)^6;        % C, in eV (A)^6
D = 73e-27*Jtoev*(nmtoA)^8;        % D, in ev (A)^8
alpha = 2.96/nmtoA ;               % alpha, in A^-1
% Energy plot
U(r) = -e^2*Jtoev/(4*pi*e0*r*Atom)+B*exp(-alpha*r)-C/r^6+D/r^8;
fplot(U(r),[1,5],'Color',[1,0,0],'LineWidth',2)
xlabel('Separation (Angstroms)')
ylabel('Energy (eV)')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
% Force plot
F = diff(U(r),r);
figure
fplot(F,[1.2,5],'Color',[1,0,0],'LineWidth',2)
xlabel('Separation (Angstroms)')
ylabel('Force (eV/Angstrom)')
grid on
axes1 = gca;
set(axes1,'FontSize',16,'XGrid','on','YGrid','on');
%
% Equilibrium separation
req = vpasolve(F==0,1.1)
%
%Stiffness - convert the units from eV/A^2 to J/m^2 (SI, therefore Newtons/m)
stif = subs(diff(F,r),r,req)/Jtoev/Atom^2
%
%Bond strength; convert from eV/A to J/m (SI, therefore Newtons)
rbreak = vpasolve(diff(F,r)==0,1.5)
Fbreak = subs(F,r,rbreak)

%
%% Problem 3
%

clear all
syms x
FL = 67.701*x^3 - 546.449*x^2+1735*x
FU = 81.032*x^3-617.95*x^2 + 1609*x
xmax = vpasolve(FL==FU,5.5)
xmax(3)
figure
p1 = fplot(FL,[0,double(xmax(3))],'LineWidth',2)
hold on
p2 = fplot(FU,[0,double(xmax(3))],'LineWidth',2)
xlabel('Stretch (m)')
ylabel('Force (N)')
title('Force-Extension curve for bungee cord')
axes1 = gca;
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
legend1 = legend([p1,p2],'Loading','Unloading');
set(legend1,...
    'Position',[0.2665983606557377 0.608333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);
% Work calculations
Wloading = double(int(FL,x,[0,6.71]))
WUnloading = -double(int(FU,x,[6.71,0]))
% Numbers
V=sqrt((20*WUnloading-195*9.81*6.71*sind(15))/195)
Efficiency = (1/2)*175*V^2/(20*Wloading)

%
%