%% Function name
%       QbarandSbar

%% Revised
%       September 29, 2013

%% Author
%       Britton Steele, Christian Muneton, Trey Moore, & Autar Kaw
%       Section: All
%       Semester: Fall 2012

%% Purpose
%       Given the angle of the ply, the the elastic modulii, 
%       Poisson's ratio, and shear modulus of a unidirectional lamina,
%       output the Qbar and Sbar matricies

%% Usage
%       function [Qbar,Sbar] = QbarandSbar(angle,moduli)
% Input variables
%       moduli=vector with four elastic modulii of unidirectional lamina
%       [moduli]=[E1 E2 nu12 G12]
%       E1=longitudinal elastic modulus
%       E2=transverse elastic modulus
%       nu12=major Poisson's ratio
%       G12=in-plane shear modulus
%       angle=angle of ply given in degrees.
% Output variables
%       [Qbar]=vector of transformed reduced stiffness matrix
%       [Sbar]=transformed compliance matrix
% Keyword
%       transformed compliance matrix 
%       transformed reduced stiffness matrix

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

%% Testing code
clc
clear all
%% Inputs
% Material: Graphite/Epoxy
moduli=[181 10.3 0.28 7.17];
fprintf('\nLongitudinal Elastic Modulus: %g ',moduli(1))
fprintf('\nTransverse Elastic Modulus: %g ',moduli(2))
fprintf('\nPoisson''s Ratio: %g ',moduli(3))
fprintf('\nShear Modulus: %g \n\n\n',moduli(4))
%% Test 1
% Testing for individual angle
% Input desired angle in degrees
angle=60;
% Call function: Qbar,Sbar
[Qbar,Sbar] = QbarandSbar(angle,moduli);
% Results for a particular angle
fprintf('\n\nAngle of Lamina: %G\n\n',angle)
disp('      Transformed Reduced Stiffness Matrix')
disp('__________________________________________________')
disp(Qbar)
disp(' ')
disp('       Transformed Compliance Matrix')
disp('__________________________________________________')
disp(Sbar)
%% Test 2
% Table of Transformed Reduced Stiffness and Transformed Compliance Matrices
% of an Angle Lamina as a Function of Fiber Angle
for i=1:5:91
    angle(i)=i-1;
    [Qbar,Sbar] = QbarandSbar(angle(i),moduli);
    Qbar11(i)=Qbar(1,1);
    Qbar12(i)=Qbar(1,2);
    Qbar16(i)=Qbar(1,3);
    Qbar22(i)=Qbar(2,2);
    Qbar26(i)=Qbar(2,3);
    Qbar66(i)=Qbar(3,3);
    Sbar11(i)=Sbar(1,1);
    Sbar12(i)=Sbar(1,2);
    Sbar16(i)=Sbar(1,3);
    Sbar22(i)=Sbar(2,2);
    Sbar26(i)=Sbar(2,3);
    Sbar66(i)=Sbar(3,3);
end
disp('                                    Transformed Reduced Stiffness Matrix')
fprintf('______________________________________________________________________')
fprintf('_____________________________________________\n')
fprintf('  angle          Qbar11             Qbar12          Qbar16')
fprintf('          Qbar22          Qbar26          Qbar66\n')
for i=1:5:91
   fprintf('\t%2.0f |\t  %E |\t %E |\t %E |\t %E |\t %E |\t %E \n',...
       angle(i),Qbar11(i),Qbar12(i),Qbar16(i),Qbar22(i),Qbar26(i),Qbar66(i))
end
%% Test 3
% Graph of Transformed Reduced Stiffness and Transformed Compliance Matrices
% of an Angle Lamina as a Function of Fiber Angle
for i=1:1:91
    angle(i)=i-1;
    [Qbar,Sbar] = QbarandSbar(angle(i),moduli);
    Qbar11(i)=Qbar(1,1);
    Qbar12(i)=Qbar(1,2);
    Qbar16(i)=Qbar(1,3);
    Qbar22(i)=Qbar(2,2);
    Qbar26(i)=Qbar(2,3);
    Qbar66(i)=Qbar(3,3);
    Sbar11(i)=Sbar(1,1);
    Sbar12(i)=Sbar(1,2);
    Sbar16(i)=Sbar(1,3);
    Sbar22(i)=Sbar(2,2);
    Sbar26(i)=Sbar(2,3);
    Sbar66(i)=Sbar(3,3);
end
figure
hold on
title('Qbars vs Fiber Angle')
xlabel('Angle (degrees)')
ylabel('Qbars')
grid
plot(angle,Qbar11,'b','LineWidth',2)
plot(angle,Qbar12,'r','LineWidth',2)
plot(angle,Qbar16,'c','LineWidth',2)
plot(angle,Qbar22,'g','LineWidth',2)
plot(angle,Qbar26,'y','LineWidth',2)
plot(angle,Qbar66,'k','LineWidth',2)
legend('Q11','Q12','Q16','Q22','Q26','Q66')
hold off
figure
hold on
title('Sbars vs Fiber Angle')
xlabel('Angle (degrees)')
ylabel('Sbars')
grid
plot(angle,Sbar11,'b','LineWidth',2)
plot(angle,Sbar12,'r','LineWidth',2)
plot(angle,Sbar16,'c','LineWidth',2)
plot(angle,Sbar22,'g','LineWidth',2)
plot(angle,Sbar26,'y','LineWidth',2)
plot(angle,Sbar66,'k','LineWidth',2)
legend('S11','S12','S16','S22','S26','S66')