// Pb 7.13-14-15-16
clear; //deletes all data from  memory
mode(0); // prints each line like MATLAB unless suppressed by ;
funcprot(0); // to avoid warnings
exec('C:/Users/EJB/OneDrive/Scilab/CLT.sci',0);

orientation = [0,45,-45,45,0]; // deg
thickness   = [1, 1, 2, 1, 1]*2.5; // mm
N = length(thickness);  // # of laminas
// properties E-glass/polyester
E1 = 204.6E3;   // MPa
E2 = 20.46E3;     // MPa
v12 = 0.3;
v23 = 0.5;      // not given
G12 = 6.89E3;    // MPa
F1t = 826.8;    // MPa
F1c = 482.3;
F2t = 68.9;
F2c = F2t;
F6 = 241.15;
// From Pb. 7.8
Q = buildQ(E1,E2,G12,v12);    
[A,B,D] = buildABD(Q,thickness,orientation)
// From Pb. 7.9
N = [1000,1000,0]'; // ' makes it a column array
eps0 = A\N;
// From Pb. 7.10-11-12
// 0 deg lamina
T=transf(0); k=1;
eps(:,k) = Reuter*T*Reuter^-1*eps0;  // 7.11
sig(:,k) = Q*eps(:,k);        // 7.12 stress in lamina c.s
// 45 deg lamina
T=transf(45); k=2;
eps(:,k) = Reuter*T*Reuter^-1*eps0;  // 7.11
sig(:,k) = Q*eps(:,k);        // 7.12 stress in lamina c.s
// -45 deg lamina
T=transf(-45); k=3; 
eps(:,k) = Reuter*T*Reuter^-1*eps0;  // 7.11
sig(:,k) = Q*eps(:,k);        // 7.12 stress in lamina c.s
sig,    // display stresses in lamina c.s. by columns
// Max. stress Pb. 7.13
FT = [F1t, F2t, F6]';
FC = [F1c, F2c, F6]';
Rk(1) = maxStress(FT,FC,sig(:,1))
Rk(2) = maxStress(FT,FC,sig(:,2))
Rk(3) = maxStress(FT,FC,sig(:,3))
Rlaminate = min(Rk)
// Max. Strain Pb. 7.14 
FT = [F1t/E1, F2t/E2, F6/G12]';
FC = [F1c/E1, F2c/E2, F6/G12]';
Rk(1) = maxStress(FT,FC,eps(:,1))
Rk(2) = maxStress(FT,FC,eps(:,2))
Rk(3) = maxStress(FT,FC,eps(:,3))
Rlaminate = min(Rk)
// Interacting coefficients Pb. 7.15
f2  = F2t^-1 - F2c^-1;  
f22 = (F2t*F2c)^-1;     
f66 = F6^-2;            
f2*1E6
f22*1E6
f66*1E6
// Interacting
for k=1:3
// fiber mode
    R_1(k) = 1 ./max(Macaulay(sig(1,k)./F1t) + Macaulay(-sig(1,k)./F1c));
    // matrix mode
    a = f22*sig(2,k)^2 + f66*sig(3,k)^2;
    b = f2*sig(2,k);
    R_2(k) = -b+sqrt(b^2+4*a)/(2*a);
    Rk(k) = min(R_1(k),R_2(k));
end
Rk
// truncated max. strain Pb. 7.16
format('v',8)   // to see small values of strain 
// calculate strains to failure if no exp. data available
eps1t = F1t/E1// longitudinal strain to failure tensile
eps1c = F1c/E1// longitudinal strain to failure compressive
eps2t = F2t/E2// transverse strain to failure tensile
eps2c = F2c/E2// transverse strain to failure compressive
// CADEC's implementation TMS
// does not check complementary lamina
for k=1:3
    if eps(1,k)>0 then
        R_1(k) = min(eps1t,eps2t)/eps(1,k);
    else
        R_1(k) = -eps1c/eps(1,k);
    end
    if eps(2,k)>0 then
        R_2(k) = min(eps1t,eps2t)/eps(2,k);
    else
        R_2(k) = -eps1c/eps(2,k);
    end
    R_6(k) = (1+v12)*max(eps1t,eps1c)/abs(eps(1,k)-eps(2,k));
    Rk(k) = min(R_1(k),R_2(k),R_6(k));
end
Rk

//table = sig' 
//table = eps'*1E6
//table(:,4) = Rk
//prettyprint(table)

