// Pb. 7-27: optimize Ex. 7.14 with TMS criterion 
// Pb. 7-28: optimize Ex. 7.14 with Interacting criterion
mode(0);clc;clear;
exec('C:\Users\EJB\OneDrive\Scilab\CLT.sci');
// Nex few lines brought from Example 7.14
d0 = 26.10;// 
alpha = 1.;//load factor
phi = 1/2;//resistance factor
c = 113000//N mm
M = c * alpha;// 
T = c * alpha;// 
t_k = 0.125;//mm 
// change next 2 lines:
//orientation =[0,45,-45,90,0,45,-45,90];// 
//thickness   =[1,1,1,1,1,1,1,1]*t_k;
orientation =[0,45,-45,0,45,-45,90];// [0_2/pm45/0/pm45/90]  
thickness   =[2,1,1,1,1,1,1]*t_k;//keep 10% rule
orientation =[0,45,-45,0,45,-45];// [0_2/pm45]_2  
thickness   =[2,1,1,2,1,1]*t_k;//violates 10% rule
// end change
n = length(thickness);
t = sum(thickness)
di = d0 - 2*t;
// New calculations begin here................................
E_1 = 171400; E_2 = 9080; G_12 = 5290; v_12 = 0.32;//Table 1.4
F_1t = 2326.2;// MPa Table 1.4
F_1c = 1200.1;// MPa Table 1.4
F_1t_all = F_1t * phi;
F_1c_all = F_1c * phi;
F_2t = 62.3;
F_2c = 199.8;
F_6 = 120;
F_2t_all = F_2t * phi;
F_2c_all = F_2c * phi;
F_6_all = F_6 * phi;
f2 = F_2t^-1 - F_2c^-1;
f22 = (F_2t*F_2c)^-1;
f66 = F_6^-2;
// e_1t = F_1t/E_1 // use the experimental value instead
e_1t = 1.62/100//Table 1.4 experimental
e_1c = F_1c/E_1// calculated, no experimental value available
e_1t_all = e_1t * phi;//allowable tensile
e_1c_all = e_1c * phi;//allowable compressive
//
Q = buildQ(E_1,E_2,G_12,v_12);
A = buildABD(Q,thickness,orientation);
sigma = - 32*M*d0/%pi/(d0^4-di^4);// compression side of the beam
tau   = 16*T*d0/%pi/(d0^4-di^4);
N = [sigma*t; 0; tau*t];
eps_ = A\N;
// Truncated Maximum Strain
R = 1E12;//initialize with a large value
for k = 1:n 
    eps = rotateEps_(eps_,orientation(k));//==> lamina cs
    R1 = Macaulay(e_1t_all/eps(1)) + Macaulay(-e_1c_all/eps(1));//fiber t/c
    R6 = (1+v_12)*max(e_1t_all,e_1c_all)/abs(eps(1)-eps(2));//fiber shear
    R = min(R1,R6,R);
end
R
// Interacting Criterion
R = 1E12; R1min=R; R2min=R;//initialize with large value
for k = 1:4
    eps = rotateEps_(eps_,orientation(k));//==> lamina cs
    sigma = Q * eps;//lamina cs
    R1 = Macaulay(F_1t_all/sigma(1)) + Macaulay(-F_1c_all/sigma(1));//fiber t/c
    a = f22*sigma(2)^2 + f66*sigma(3)^2;
    b = f2*sigma(2);
    R2 = -b+sqrt(b^2+4*a)/2/a;
    R1min = min(R1min,R1);//track fiber failure
    R2min = min(R2min,R2);//track matrix-dominated failure
    R = min(R1,R2,R);//track overal strength ratio
end
R, R1min, R2min
