// Pb. 4.19
mode(0)
function g = g(alpha,s2,s6,F4,F6,nT,nL)
    sn = s2*cosd(alpha).^2
    tauT = -s2*sind(alpha).*cosd(alpha)
    tauL = s6*cosd(alpha)
    g = (tauT./(F4-nT*sn)).^2 + (tauL./(F6-nL*sn)).^2 - 1
    g = -g // for fminsearch to find the maximum
endfunction
// Glass/LY556 Table 1.3
alp_0 = 53.0;   // degrees alpha zero
F_2c  = 138.0;    // [MPa] at alpha_zero
F_6   = 61.0;     // [MPa]
//
F_4 = 138 * cosd(alp_0) * ( sind(alp_0) + cosd(alp_0)*cotd(2*alp_0) )
eta_T = -cotd(2*alp_0)
eta_L = -(F_6*cosd(2*alp_0)/(F_2c*cosd(alp_0)^2))
// choose the loading
s_2 = -60   //[MPa]
s_6 =  80.255   // [MPa]
// find the minimum of the negative of g
options = optimset("TolFun",1E-12)
[alpha,min_g] = fminsearch(list(g,s_2,s_6,F_4,F_6,eta_T,eta_L),53)
max_g = - min_g  // minus to recover max(g)
// plot the actual g_ = minus g
close()			// in case a figure is open
x = [0:1:90]    // alpha
g_ = - g(x,s_2,s_6,F_4,F_6,eta_T,eta_L) // minus to recover the actual g
plot(x, g_, 'thickness',2)
xlabel('$\text{Fracture plane angle }\alpha$','font_size',4)
ylabel('Activation function g','font_size',3)
legend('g','I_F')
a=gca()
a.font_size=2
l=a.children(1) // legend handle
l.font_size=3
//a.data_bounds=[0,90,0,0.3]	// set axis range
xs2png(gcf(),'Pb417.png') // save plot as a PNG file
I_F = max_g + 1
disp(['alpha, g, I_F, R = ',string(alpha),string(max_g),string(I_F),string(1/I_F)])
