function matlab_intro_part_2_assignment
%% MATLAB: Introduction - Part 2: Assignment
clear all
close all
clc


%% Problem 1

y = randn(2,3,40,12);

% Part (a)
ya = y(:,:,1,1);
size(ya)

% Part (b)
yb = y(:,:,1:end/2,5);
size(yb)

% Part (c)
y(1,1:2,1:10,3:5) = exp(1);

% Part (d)
try
    y(2,:,:,12) = ones(3,40);
    disp('Problem 1 (d): No exception caught');
catch EX
    disp(['Problem 1 (d): Exception caught: ' EX.message]);
end

% Part (e)
try
    y(:,:,15,1) = eye(3);
    disp('Problem 1 (e): No exception caught');
catch EX
    disp(['Problem 1 (e): Exception caught: ' EX.message]);
end


%% Problem 2

N = 50;
s = 0;
for i = 1:N
    s = s + 1/i + 1/((i + 2)*(i + 3));
end
sprintf('Sum = %g',s)


%% Problem 3

n = 8;
H = zeros(8); % Preallocate H

for i = 1:n
    for j = 1:n
        H(i,j) = 1/(i + j - 1);
    end
end
H


%% Problem 4

x = linspace(-3,3,200);
y = zeros(size(x)); % Preallocate y

for i = 1:length(y)
    y(i) = prob4f(x(i));
end

figure('Name','Problem 4','Color','white');
plot(x,y);
xlabel('x');
ylabel('f(x)');


%% Problem 5

u = cmu.units;

% Methane component
methane.Pc = 46.4*u.bar;
methane.Vc = 99E-6*u.m^3/u.mol;
methane.a = 3*methane.Pc*methane.Vc^2;
methane.b = methane.Vc/3;

% Ranges for T and V
Tmin = 50*u.K;
Tmax = 600*u.K;
Vmin = 10E-4*u.m^3/u.mol;
Vmax = 200E-4*u.m^3/u.mol;

% Calculate P based on a meshgrid of T and V
T = linspace(Tmin,Tmax,300).';
V = linspace(Vmin,Vmax,300).';
[TT VV] = meshgrid(T,V);
P = vdWP(methane,VV,TT);

% Plot P with contours
figure('Name','Problem 5: Surface Plot','Color','white');
meshc(TT,VV,P/u.bar);
title('PVT Diagram for Methane');
xlabel('T [K]');
ylabel('V [m^3/mol]');
zlabel('P [bar]');

% Isotherms
Tiso = [100; 190; 250; 400]*u.K;
Piso = zeros(length(V),length(Tiso))*u.bar;

% Calculate P for each Tiso and plot the isotherm
figure('Name','Problem 5: Isotherms Plot','Color','white');
legendcell = cell(size(Tiso));
colors = lines(length(Tiso));
for i = 1:size(Piso,2)
    Piso(:,i) = vdWP(methane,V,Tiso(i));
    
    plot(V,Piso(:,i)/u.bar,'Color',colors(i,:));
    hold on;
    legendcell{i} = sprintf('T = %g',Tiso(i));
end
hold off;

title('Isotherms for Methane');
xlabel('V [m^3/mol]');
ylabel('P [bar]');
legend(legendcell{:});


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fx = prob4f(x)
if (x < 0)
    fx = 0;
elseif (x >= 0 && x <= 1)
    fx = x;
elseif (x >= 1 && x <= 2)
    fx = 2 - x;
else
    fx = 0;
end


function P = vdWP(component,V,T)
u = cmu.units;
R = 8.314472E-5*u.m^3*u.bar/(u.K*u.mol);

P = R*T./(V - component.b) - component.a./V.^2;