clc
clear all


mu_1 = 1;
rho_1 = 1;


% HYPERBOLIC 
b1 = @(a,b,t) (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
b2 = @(a,b,t) mu_1.*(rho_1/8).*(1-t.*(1+b.*a));

h = @(a,b,t) b1(a,b,t).*b2(a,b,t);


% fix alpha
% a=1;


% fix theta , value in (0,1)
theta= 0.55;




% ELLIPTIC 

q1 = @(b,t) mu_1.*(b./(t+(1-t).*b)); %harmonic mean
q2 = @(b,t) mu_1.*((1-t)+t.*b);    % mu_bar
q3 = @(b,t) q1(b,t);

b1 = @(a,b,t) (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
b2 = @(a,b,t) mu_1.*(rho_1/8).*(1-t.*(1+b.*a));

a1 = @(a,b,t) (q2(b,t).*q1(b,t).*b1(a,b,t) - q3(b,t).*q2(b,t).*b2(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
a2 = @(a,b,t) (q1(b,t).*q2(b,t).*b2(a,b,t) - q3(b,t).*q1(b,t).*b1(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 

e = @(a,b,t) a1(a,b,t).*a2(a,b,t);


x = -20:0.2:20;
y = 0:0.1:20;
[X,Y] = meshgrid(x,y);
T = theta*ones(size(X));


V = e(X,Y,T);

V2 = h(X,Y,T);

% COLOR-Test
% C = double((V>=0));
% 
% surf(X,Y,V,C, 'FaceAlpha',0.5,'EdgeColor','none') 
% xlabel('alpha');
% ylabel('beta');
% axis([-20 20 0 20 -100 100])
% % colorbar
% hold on
% 
% C2 = double((V2>=0));
% surf(X,Y,V2,C2,'FaceAlpha',0.5,'EdgeColor','none') 
% xlabel('alpha');
% ylabel('beta');
% axis([-20 20 0 20 -100 100])
% colorbar
% mycolors = [1 0 0 ; 0 0 1];
% colormap(mycolors);
%  % view(90,0)



% %plot values above zero 

surf(X,Y,(V>=0).*V,'FaceAlpha',0.5,'EdgeColor','none','FaceColor','r') 
xlabel('alpha');
ylabel('beta');
axis([-20 20 0 20 -100 100])
colorbar
hold on

surf(X,Y,(V<0).*V,'FaceAlpha',0.5,'EdgeColor','none','FaceColor','b') 
xlabel('alpha');
ylabel('beta');
axis([-20 20 0 20 -100 100])
colorbar
% view(90,0)

hold on
surf(X,Y,(V2>=0).*V2,'FaceAlpha',0.5,'EdgeColor','none','FaceColor','black') 
xlabel('alpha');
ylabel('beta');
axis([-20 20 0 20 -100 100])
colorbar
hold on
surf(X,Y,(V2<0).*V2,'FaceAlpha',0.5,'EdgeColor','none','FaceColor','g') 
xlabel('alpha');
ylabel('beta');
axis([-20 20 0 20 -100 100])
colorbar