// ASYMPTOTIC STATE BOUNDARY SURFACE (ASBS) OF CLAY HYPOPLASTICITY (Masin, 2013)
// stress space
// Project: Animating Soil Models https://soilmodels.com/soilanim
// by: Gertraud Medicus, University of Innbruck, 06/2023
//
// Software: asymptote.sourceforge.io
// Asymptote Web Application: http://asymptote.ualberta.ca/

import graph3;
import three;
import patterns;
import texcolors;

//
size3(200,200,200);
currentprojection=orthographic(-40,100,70); 

//
scale(Linear,Linear,Linear);
defaultpen(Symbol(series="l",shape="n"));
defaultpen(linewidth(.5));


// hypo clay (Masin, 2013) parameters
real phi=25*pi/180;
real M=6*sin(phi)/(3-sin(phi));
real N= 1;
real lambda=0.1;
real kappa=0.01;
real a  =0.3;
//
real p=110;
real T22=-100;
real T33=T22;
real pe=100;
real v_plot=exp(N-lambda*log(pe)); 

////////////////////////
triple f(pair theta) {
real phi_m=theta.y*pi/180; 
real Kc=(1+sin(phi_m))/(1-sin(phi_m));
real T11=Kc*T22;
real I1=T11+T22+T33;
//
real I2= 1/2*(T11*T11+T22*T33*2-I1^2);
real I3=T11*T22*T33;
real Fm=(9*I3+I1*I2)/(I3+I1*I2);
real omega=-log(cos(phi)^2)/log(2)+a*(Fm-sin(phi)^2);
real p_zu_pe=(1-Fm)^(1/omega);
real q_zu_p=3*(T11-T22)/(T11+T22+T33);
real q_zu_pe=p_zu_pe*q_zu_p;
real p_p=pe*p_zu_pe;

real Kmn=(9-sin(phi_m)^2)/(1-sin(phi_m)^2);

// MN Griffith f=0 cubic equation
// ax3 + bx2 + cx + d = 0.
real a= -sqrt(2)*Kmn*sin(3*theta.x); //a3
real b=(9-3*Kmn)*p_p*sqrt(3); //a2
real c=0; //a1
real d=(2*Kmn-18)*(p_p*sqrt(3))^3; //a0

real [] t=cubicroots(a,b,c,d);
//write(t);
// Remove negative or zero values of t:
t=t>= 0 ? t : null;

real t=min(t);
write(t);

////// MATSUOKA NAKA 
real L=2*sqrt(2)*sin(phi_m)/(3-sin(phi_m)); //DP
return (
//MN
(p_p+sqrt(2/3)*t*sin(theta.x-2*pi/3)),
(p_p+sqrt(2/3)*t*sin(theta.x)),
(p_p+sqrt(2/3)*t*sin(theta.x+2*pi/3))
//   Drucker-Prager:
// (p_p+sqrt(2)*L*p_p*sin(theta.x-2*pi/3)),
//(p_p+sqrt(2)*L*p_p*sin(theta.x)),
//(p_p+sqrt(2)*L*p_p*sin(theta.x+2*pi/3))
);
}
surface asbs=surface(f,(-179.9*pi/180,.1),(179.9*pi/180,89.9),100,30,Spline);
draw(asbs,yellow +opacity(0.9),meshpen=black,render(compression=Low,merge=true));



/////////////////////
//CRITICAL STATE -> Red line
/////////////////////
real Kmn=(9-sin(phi)^2)/(1-sin(phi)^2);
real pc=pe/2;
real s=3*pc/sqrt(3);

real x(real theta) {
// MN Griffith f=0 kubische GL
// ax3 + bx2 + cx + d = 0.
real a= -sqrt(2)*Kmn*sin(3*theta); //a3
real b=(9-3*Kmn)*s; //a2
real c=0; //a1
real d=(2*Kmn-18)*s^3; //a0

real [] t=cubicroots(a,b,c,d);
t=t>= 0 ? t : null;
real t=min(t);
return (s/sqrt(3)+sqrt(2/3)*t*sin(theta-2*pi/3));}
real y(real theta) {
// MN Griffith f=0 kubische GL
// ax3 + bx2 + cx + d = 0.
real a= -sqrt(2)*Kmn*sin(3*theta); //a3
real b=(9-3*Kmn)*s; //a2
real c=0; //a1
real d=(2*Kmn-18)*s^3; //a0

real [] t=cubicroots(a,b,c,d);
t=t>= 0 ? t : null;
real t=min(t);
return (s/sqrt(3)+sqrt(2/3)*t*sin(theta));}
real z(real theta) {
// MN Griffith f=0 kubische GL
// ax3 + bx2 + cx + d = 0.
real a= -sqrt(2)*Kmn*sin(3*theta); //a3
real b=(9-3*Kmn)*s; //a2
real c=0; //a1
real d=(2*Kmn-18)*s^3; //a0

real [] t=cubicroots(a,b,c,d);
t=t>= 0 ? t : null;
real t=min(t);
return (s/sqrt(3)+sqrt(2/3)*t*sin(theta+2*pi/3));}

path3 Critline=graph(x,y,z,-179.9*pi/180,179.9*pi/180,operator ..);
draw(Critline,red+0.75pt);


// axes
/////////////////////////////////////////////////////////////////////////////
//$p' \sqrt{3} $
draw(Label("",1.2),(0,0,0) -- (p*1.1,p*1.1,p*1.1),blue);
real sig=80;
//draw(Label("$q \sqrt{2/3} $",1),(0,0,0) -- (-sig/2,-sig/2,sig),red,Arrow3);

//// ende
draw(Label("$\sigma_2'$",1),(-.2*p,0,0)--(1.5*p,0,0),black,Arrow3);
draw(Label("$\sigma_3'$",1),(0,-.2*p,0)--(0,1.5*p,0),black,Arrow3);
draw(Label("$\sigma_1'$",1),(0,0,-.2*p)--(0,0,1.5*p),black,Arrow3);