
%% CALCULATING COMPRESSION MODULUS %%

% givens
E1 = 3; E2 = 5; %GPa
h = [1,1,1,1]; % h should be N-length array based on N-layers

% for 2 layers, effective modulus:

E_eff_2 = (h(1) + h(2))/(h(1)/E1 + h(2)/E2)

% for N-Layers, effective modulus: 

E_eff_N = (E1*E2*sum(h))/((sum(h(1:2:end))*E2)+(sum(h(2:2:end))*E1))

% denominator sums the even and odd elements of h 
% corresponding to the different material layers


%% EXPLORING PARAMETERS %%

% how does the number of layers impact E_eff?
N = 100; % modify this

h = ones(N,1);
E_eff = [E1];
for i = [2:length(h)]
    E_eff(i) = (E1*E2*sum(h(1:i))/((sum(h(1:2:i))*E2)+(sum(h(2:2:i))*E1)));
end

plot(1:N,E_eff)

% observation: as number of layers increases, E_eff oscillates around
% an constant value and amplitude of oscillations diminishes

%% how does changing layer thickness uniformly impact E_eff?

t1 = 2;
t2 = 10;
h1 = t1*ones(N,1);
h2 = t2*ones(N,1);
E_eff = [E1];
for i = [2:length(h1)]
    E_eff1(i) = (E1*E2*sum(h1(1:i))/((sum(h1(1:2:i))*E2)+(sum(h1(2:2:i))*E1)));
end
  
for i = [2:length(h2)]
    E_eff2(i) = (E1*E2*sum(h2(1:i))/((sum(h2(1:2:i))*E2)+(sum(h2(2:2:i))*E1)));
end

plot(1:N,E_eff1); hold on
plot(1:N,E_eff2)

% observation: when thicknesses are kept constant for the 2 materials, no
% change in E_eff is observed from the previous case

%% how does changing layer thickness for one material impact E_eff?

N = 100;
t = 10;
h = ones(N,1);
h(2:2:end) = t; % changing only the thickness of all layers with E2

E_eff = [E1];
for i = [2:length(h)]
    E_eff(i) = (E1*E2*sum(h(1:i))/((sum(h(1:2:i))*E2)+(sum(h(2:2:i))*E1)));
end

plot(1:N,E_eff)

% observation: the E_eff value shifts drastically towards the E value of 
% the thicker material, as expected

%% CALCULATING BENDING MODULUS %%

% givens
E1 = 3; E2 = 5; % GPa
z = [3, 5, 6, 8]; 
N = length(z);
% z should be N-length array based on N-layers
% z(N) = total height of the stack

% for N-Layers, calculate effective modulus:
zc_num = [E1*(z(1)^2)];
zc_den = [E1*(z(1))];

for i = [2:length(z)]
    if mod(i,2) == 0 % for even items, using E2
        zc_num(i) = E2*(z(i)^2-z(i-1)^2);
        zc_den(i) = E2*(z(i)-z(i-1));
    else % for odd items, using E1
        zc_num(i) = E1*(z(i)^2-z(i-1)^2);
        zc_den(i) = E1*(z(i)-z(i-1));
    end
end

zc = sum(zc_num)/(2*sum(zc_den));

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

% compute final E_eff bending value !
E_eff = 4 * sum(E)

%% EXPLORING PARAMETERS %%

% how does the number of layers impact E_eff?
N = 100; % modify this

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

for i = [2:length(z)]
    zc_num = [E1*(z(i)^2)];
    zc_den = [E1*(z(i))];

    for i = [2:length(z)]
        if mod(i,2) == 0 % for even items, using E2
            zc_num(i) = E2*(z(i)^2-z(i-1)^2);
            zc_den(i) = E2*(z(i)-z(i-1));
        else % for odd items, using E1
            zc_num(i) = E1*(z(i)^2-z(i-1)^2);
            zc_den(i) = E1*(z(i)-z(i-1));
        end
    end

    zc = sum(zc_num)/(2*sum(zc_den));

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

    % compute final E_eff bending value !
    E_eff(i) = 4 * sum(E); 
end

plot(1:N,E_eff)

