clear;clc; close all
Emxx = 3.5e9                    %Pa
Egxx = 1.0e12                   %Pa
tg = 1e-9                       %m espesor del rve
lrve = 1000e-9                  %m longitud RVE
hrve = lrve                     %m altura RVE
Ncells = 10                     %numero de celdas que conforman el 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

%calcula la funcion que determinar el modulo elasticos para 
%distintas concentraciones de grafeno
k = 0;
wgi = 1e-4
wgo = 0.9999
dwg = 1e-2
for wg = wgi : dwg : wgo
  k = k + 1;
  Exx(k) = CalcYoungMod(Ncells,l,h,wg,Emxx,Egxx,tg);
  wgg(k) = wg;
end

semilogy(wgg,Exx/1e9,'b-','linewidth',1.5)
xlabel('PORCENTAJE DE GRAFENO')
ylabel('MODULO DE YOUNG (GPA)')
grid minor on
  
%calcula modulo elástico equivalente para cada capa, promediando valores en las capas
thk = 0.1                   %espesor de la placa
nlayers = 15                %numero de capas
thk_ly = thk/nlayers        %espesor de cada capa

zo = 0.0;
for i = 1 : nlayers
  zf = zo + thk_ly;
  %la relacion entre la concentracion wg del grafeno y 
  %la distancia desde el borde inferior de la placa esta
  %dada por una relacion lineal del tipo: wgz = (wgo-wgi)/thk*z + wgi
  %donde wgi es la concentracion en la superficie inferior
  %wgo es la concentracion en la superficie exterior y wgz es el valor
  %de la concentracion a una distancia z desde el borde inferior.
  %En la siguiente linea se calcula el promedio Exx_lyr para 1 capa
  wgzo = (wgo-wgi)/thk*zo + wgi;
  wgzf = (wgo-wgi)/thk*zf + wgi; 
  ExxPoly = spline(wgg,Exx,[wgzo:(wgzf-wgzo)/100:wgzf]);
  Exx_lyr(i) = 1/(wgzf-wgzo)*trapz([wgzo:(wgzf-wgzo)/100:wgzf],ExxPoly);
  zo = zf;
end

Exx_lyr'

%calculo del modulo de rigidez para cada capa


%calculo de modulo de poisson para cada capa
nu_lyr = Exx_lyr/(2*G)-1


