z= menu('hill yield criteria','from yield strength','from r-values');
if z== 1
sig0 = input('Enter the value of yield stress along rolling direction:');
angle1 = (pi/180)*(input('Enter the value of angle 1:'));
sig1 = input('Enter the value of yield stress along angle 1:');
angle2 = (pi/180)*(input('Enter the value of angle 2:'));
sig2 = input('Enter the value of yield stress along angle 2:');
sigbiaxial = input('Enter the value of yield stress in biaxial condition:');
c= [1 (sig0/sig1)^2 (sig0/sig2)^2 (sig0/sigbiaxial)^2];
A = [0,1,1,0;(sin(angle1))^4,(cos(angle1))^4,((cos(angle1))^2 - (sin(angle1))^2)^2,2*(sin(angle1))^2 * (cos(angle1))^2;(sin(angle2))^4,(cos(angle2))^4,((cos(angle2))^2 - (sin(angle2))^2)^2,2*(sin(angle2))^2 * (cos(angle2))^2;1,1,0,0];
B = [1;c(2);c(3);c(4)];
X = linsolve(A,B);
q = msgbox({['F = ',num2str(X(1))]
['G = ',num2str(X(2))]
['H = ',num2str(X(3))]
['N = ',num2str(X(4))]});
a=(X(2)+X(3))/(sig0)^2;b=(X(3)+X(1))/((sig0)^2);c=-(X(3))/((sig0)^2);d=(2*X(4))/((sig0)^2);
alp = (b/(a*b-c^2))^0.5;
beta = (a/(a*b-c^2))^0.5;
x= linspace(-alp,alp,100);
subplot(1,2,1)
plot(x,(-2*c*x)./(2*b) + ((((2*c*x).^2-4*b*(a*x.^2-1)).^0.5)./(2*b)),'r',x,(-2*c*x)./(2*b) - ((((2*c*x).^2-4*b*(a*x.^2-1)).^0.5)./(2*b)),'r','LineWidth',2)
grid on
title('Hill48 Yield Criteria')
axis([-1.2*alp 1.2*alp -1.2*beta 1.2*beta])
subplot(1,2,2)
alp1 = (round((1.2*((b/(a*b-c^2))^0.5))/10))*10;
beta1 = (round((1.2*((a/(a*b-c^2))^0.5))/10))*10;
l=gcd(alp1,beta1);
gama1=10*(l);
fun=@(x,y,z)(a*X.^2+b*Y.^2+d*(Z).^2+c*(X.*Y)-1) ;
[X,Y,Z]=meshgrid(linspace(-alp1,alp1,l+1),linspace(-beta1,beta1,l+1),linspace(-gama1,gama1,l+1))
val=fun(X,Y,Z);
fv=isosurface(X,Y,Z,val,0);
p = patch(fv);
isonormals(X,Y,Z,val,p)
set(p,'FaceColor' , 'red');
set(p,'EdgeColor' , 'none');
daspect([1,1,1])
view(3); axis tight
camlight
lighting phong