clear;clc; close all

%----------------- datos de entrada ---------------------
Em = 0.675e9                    %Pa
Eg = 0.1013e6                   %Pa
tg = 1.5e-5                     %m espesor del rve
lrve = 1.5e-5               %m longitud RVE
Ncells = 10                   %numero de celdas que 
                              %conforman el RVE
%--------------------------------------------------------

wgi = 0.01;
dwg = 0.001;
wgf = 0.99;
hrve = lrve                     %m altura RVE
l = lrve/Ncells                 %m logitud de celda que rodea a grafeno l*Ncells = long RVE
h = l;                          %m altura de celda que rodea a grafeno h*Ncells = altura RVE

ki=0.8;
dk=0.05;
kf = 1.1;

w = 0;
Exxlyr = [];
for k = ki:dk:kf
  w = w + 1;  
  s=0;
  for wg = wgi:dwg:wgf    
    s = s + 1;    
    agrafcell = l*h*wg;                       %area del grafeno      
    lg = sqrt(k*agrafcell);
    hg = lg/k;
    Exxlyr(w,s) = CalcYoungMod_4(Ncells,l,h,Em,Eg,tg,wg,lg,hg);
  end
end

k = [ki:dk:kf];
for i = 1 : w
  id{i} = strcat("l_g/h_g = ",num2str(k(i)));
end

figure(1)
plot([wgi:dwg:wgf],Exxlyr/1e9,'-','linewidth',1.5)
legend(id(1,:),"location",'northwest','numcolumns',1)
xlabel('GRAPHENE WEIGHT - Wg')
ylabel('YOUNG MODULUS - [GPa]')
#grid minor on

w = [wgi:dwg:wgf];
k = 0
for i = 1 : 200 : s
  k = k + 1;
  idw{k} = strcat("w_g = ",num2str(w(i)));
end

figure(2)
plot([ki:dk:kf],Exxlyr(:,[1:100:end])'/1e9,'-','linewidth',1.5)
legend(idw(1,:),"location",'northwest','numcolumns',1)
xlabel('GRAPHENE RATIO - l_g/h_g')
ylabel('YOUNG MODULUS - [GPa]')
#grid minor on