%% COMPRESSION MODULUS

clc
clear
close all

% Step 1: defination of variables
h_eff = (0: .01: 1); % total h_eff values of material 1
h_eff = h_eff'; 
[r, c] = size(h_eff); % define the size of the h_eff matrix
h2_eff = (ones([r, c]))-h_eff; % create a matrix of h2_eff values
Eratio = [0.01; 0.1; 1; 10; 100]; % define possible values of E_ratio (can do this since Eratio is dimensionless)
B_bar = zeros(r, length(Eratio)); % define initial Matrix

% iterate through the columns of B_bar to update it to contain the relevent
% value
for i = 1:5
    B_bar(:,i) = 1./(h_eff + Eratio(i)*h2_eff);
    
end

% plot B_bar
figure()
hold on
plot(h_eff, B_bar(:,1:5), "Linewidth", 2)
xlim([0 1])
ylim([0 10])
xlabel("Effective 1/h")
ylabel("Normalized compression modulus")
legend('E{ratio} = 0.01','E{ratio} = 0.1', 'E{ratio}=1', 'E{ratio} = 10', 'E{ratio} = 100')
title("Normalized compression modulus as a funciton of Effective 1/h")


%% Bending Modulus

% initialize values
N = round(rand()*10); % generate a number for the number of layers
z_ratio = linspace(0, 1, N); % ratio of distance at the neutral axis to the distance in question (represented as a ratio here, but calculated for a final neutral axis location later)
Eratio = [0.01; 0.1; 1; 10; 100]; % possible values of the E ratio
h = ones(size(N)); % an array of ones representing the total length
z = [round(rand(5))*10]; % generate many random numbers for z
E1 = 5; % randomly assign a value to E1 and E2
E2 = 7;

% initialize expressions
z_numerator = [E1*(z(1)^2)];
z_denominator = [E1*(z(1))];

% calculate what the distance of the neutral axis is
for i = [1:length(z)]
    if mod(i,2) == 0 % check to see if even, and if so, use E2
        z_numerator(i) = E2*(z(i)^2-z(i-1)^2);
        z_denominator(i) = E2*(z(i)-z(i-1));
    if mod(i,2) ~= 0 % check to see if odd, and if so, use E1
        z_numerator(i) = E1*(z(i)^2-z(i-1)^2);
        z_denominator(i) = E1*(z(i)-z(i-1));
    end
    end
end

% Using the information from what is above, calculate the distance to the
% neutral axis
zc = sum(z_numerator)/(2*sum(z_denominator));


E = [E1 * ((z(1)/z(N) - zc/z(N))^3)];

for i = [2:length(z)]
    if mod(i,2) == 0 % calculate for even layers (E2)
        E(i) = E2 * ((z(i)/z(N) - zc/z(N))^3 - (z(i-1)/z(N) - zc/z(N))^3);
    if mod(i, 2) ~=0 % calculate for odd layers (E1)
        E(i) = E1 * ((z(i)/z(N) - zc/z(N))^3 - (z(i-1)/z(N) - zc/z(N))^3);
    end
    end
end

% calculate the effective final values
E_eff_final=4 * sum(E)