%% Function name
%       moduli_elasticity

%% Revised
%       28 January 2014

%% Authors
%       Mateusz Malinowski, Trey Moore, & Autar Kaw
%       Section: All
%       Semester: Fall 2013

%% Purpose
%       Given the elastic modulii, Poisson's ratios, of a fiber and matrix,
%       as well as the fiber volume fraction, output the elastic modulii
%       and Poisson's ratios of the lamina using the elasticity method

%% Usage 
%       function [moduli] = moduli_elasticity(Ef,Em,nuf,numm,Vf)
% Input Variables
%       Ef=fiber elastic modulus
%       Em=matrix elastic modulus
%       nuf=fiber Poisson's ratio
%       numm=matrix Poisson's ratio
%       Vf=fiber volume fraction
% Output Variables 
%       moduli=vector of elastic modulii and Poisson's ratio of unidirectional lamina
%       [E1 E2 nu12 G12]
% Keyword
%       elastc modulii
%       Poisson's ratio
%       elasticity method

%% License Agreement
%       http://www.eng.usf.edu/~kaw/OCW/composites/license/limiteduse.pdf

%% Code 
function [moduli] = moduli_elasticity(Ef,Em,nuf,numm,Vf)
% Longitudianl elastic modulus
E1=Ef*Vf+Em*(1-Vf)-((2*Em*Ef*Vf*(1-Vf)*(nuf-numm)^2)/(Ef*(2*numm^2*Vf-...
    numm+Vf*numm-Vf-1)+Em*(-1-2*Vf*nuf^2+nuf-Vf*nuf+2*nuf^2+Vf)));
% Major Poisson's ratio
nu12=nuf*Vf+numm*(1-Vf)+((Vf*(1-Vf)*(nuf-numm)*(2*Ef*numm^2+numm*Ef-Ef+...
    Em-Em*nuf-2*Em*nuf^2))/((2*numm^2*Vf-numm+numm*Vf-1-Vf)*Ef+...
    (2*nuf^2-Vf*nuf-2*Vf*nuf^2+Vf+nuf-1)*Em));
% Transverse elastic modulus
Nf=3-4*nuf;
Nm=3-4*numm;
% Bulk modulus of the fiber under longitudinal plane strain
Kf=((Ef)/(2*(1+nuf)*(1-2*nuf)));
% Bulk modulus of the matrix under longitudinal plane strain
Km=((Em)/(2*(1+numm)*(1-2*numm)));
% Shear elastic modulus of matrix
Gm=Em/(2*(1+numm));
% Shear elastic modulus of fiber
Gf=Ef/(2*(1+nuf));
% Bulk modulus of lamina under longitudinal plane strain
K_s=((Km*(Kf+Gm)*(1-Vf)+Kf*(Km+Gm)*Vf)/((Kf+Gm)*(1-Vf)+(Km+Gm)*Vf));
% Constants of quadratic equation to solve for shear modulus
A=3*Vf*(1-Vf)^2*((Gf/Gm)-1)*((Gf/Gm)+Nf)+((Gf/Gm)*Nm+Nf*Nm-...
    ((Gf/Gm)*Nm-Nf)*Vf^3)*(Vf*Nm*((Gf/Gm)-1)-((Gf/Gm)*Nm+1));
B=-3*Vf*(1-Vf)^2*((Gf/Gm)-1)*((Gf/Gm)+Nf)+(1/2)*(Nm*(Gf/Gm)+...
    ((Gf/Gm)-1)*Vf+1)*((Nm-1)*((Gf/Gm)+Nf)-2*((Gf/Gm)*Nm-Nf)*Vf^3)+...
    (Vf/2)*(Nm+1)*((Gf/Gm)-1)*((Gf/Gm)+Nf+((Gf/Gm)*Nm-Nf)*Vf^3);
C=3*Vf*(1-Vf)^2*((Gf/Gm)-1)*((Gf/Gm)+Nf)+(Nm*(Gf/Gm)+...
    ((Gf/Gm)-1)*Vf+1)*((Gf/Gm)+Nf+((Gf/Gm)*Nm-Nf)*Vf^3);
% Shear modulus quadratic equation
G23_1=((((-2*B)/Gm)+sqrt((((2*B)/(Gm))^2)-4*(A/Gm^2)*C))/(2*(A/Gm^2)));
G23_2=((((-2*B)/Gm)-sqrt((((2*B)/(Gm))^2)-4*(A/Gm^2)*C))/(2*(A/Gm^2)));
% Solving for m
m=1+4*K_s*((nu12^2)/E1);
% Conditional statement for elastic modulii and Poisson's ratio selection
if G23_1>=0
    nu23=((K_s-m*G23_1)/(K_s+m*G23_1));
    E2=2*(1+nu23)*G23_1;
else
    nu23=((K_s-m*G23_2)/(K_s+m*G23_2));
    E2=2*(1+nu23)*G23_2;    
end
% Shear modulus
G12=Gm*((Gf*(1+Vf)+Gm*(1-Vf))/(Gf*(1-Vf)+Gm*(1+Vf)));
% Vector of elastic modulii and Poisson's ratio of unidirectional lamina
[moduli]=[E1;E2;nu12;G12];
end