// Pb 4.20 - just plot it
clear;clc;
mode(0); // prints each line like MATLAB unless suppressed by ;
exec('C:/Users/EJB/OneDrive/Scilab/CLT.sci',0);//See Appendix

function g = MohrCoulomb(alpha,sigma2,sigma6,F4,F6,etaT,etaL)
    sigmaN = sigma2*cosd(alpha).^2;
    tauT = -sigma2*sind(alpha).*cosd(alpha);
    tauL = sigma6*cosd(alpha);
    // minus for fminsearch to find the maximum
    g = - ( (tauT./(F4-etaT*sigmaN)).^2 + (tauL./(F6-etaL*sigmaN)).^2 - 1 );
endfunction

F6 = 61;//MPa
F2c=138;//MPa
alpha0=53;//deg
sigma2 = -90;
sigma6 = 90;

F4 = F2c*cosd(alpha0)*(sind(alpha0)+cosd(alpha0)*cotd(2*alpha0))
etaT = -cotd(2*alpha0)
etaL = -F6*cosd(2*alpha0)/F2c/cosd(alpha0)^2

//plot
alpha = [0:1:53]
plot(alpha,-MohrCoulomb(alpha,sigma2,sigma6,F4,F6,etaT,etaL))
xtitle('','$\alpha$','$g(\alpha)$')
xs2png(gcf(),'Pb420.png');

// serch for the maximum
[alpha,g] = fminsearch(list(MohrCoulomb,sigma2,sigma6,F4,F6,etaT,etaL),40)
sigmaN = sigma2*cosd(alpha).^2
tauT = -sigma2*sind(alpha).*cosd(alpha)
tauL = sigma6*cosd(alpha)


