clc
clear all

syms b1(rho_1,alpha,beta,theta)
syms b2(rho_1,alpha,beta,theta)

b1(rho_1,alpha,beta,theta) = (3.*rho_1./2).*beta.*(1-theta.*(1+alpha));

b2(rho_1,alpha,beta,theta) = ((3.*rho_1)./(4.*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));



rho_1 = 1.0;
theta= 0.25;
alpha = 2;
beta = 5;

x = 0:0.01:1;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
figure
plot(x,y)
hold on 
alpha = -10;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
alpha = -1;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
alpha = 10;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
legend('alpha=2', 'alpha=-10', 'alpha=-1', 'alpha=10')
title(['b2 for beta =', num2str(beta)])
xlabel('theta')
ylabel('b2')


alpha = 2;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
figure
plot(x,y)
hold on 
alpha = -10;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
alpha = -1;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
alpha = 10;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
legend('alpha=2', 'alpha=-10', 'alpha=-1', 'alpha=10')
title(['b1 for beta =', num2str(beta)])
xlabel('theta')
ylabel('b1')





%
%


alpha = 2;
beta = 0.5;

x = 0:0.01:1;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
figure
plot(x,y)
hold on 
beta = 1;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 2;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 5;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 10;
y = arrayfun(@(t)double(b2(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
legend('beta=0.5', 'beta=1', 'beta=2', 'beta=5', 'beta=10')
title(['b2 for alpha =', num2str(alpha)])
xlabel('theta')
ylabel('b2')


beta = 0.5;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
figure
plot(x,y)
hold on 
beta = 1;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 2;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 5;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
hold on 
beta = 10;
y = arrayfun(@(t)double(b1(rho_1,alpha,beta,t)),x,'UniformOutput', false);
y = cell2mat(y);
plot(x,y)
legend('beta=0.5', 'beta=1', 'beta=2', 'beta=5', 'beta=10')
title(['b1 for alpha =', num2str(alpha)])
xlabel('theta')
ylabel('b1')

%
%
%
%


figure
fsurf(@(alpha,beta) b1(rho_1, alpha ,beta, theta ), [-20  20 0  20]  ) 
title(['Fixed Theta Plot: theta = ', num2str(theta) ])
xlabel('alpha')
ylabel('beta')

figure
% fsurf(@(alpha,beta) b2(rho_1, alpha ,beta, theta) ) 
% singularity if beta = -(1-theta)/theta  ...beta is only positive
%sing = -(1-theta)/theta  
fsurf(@(alpha,beta) b2(rho_1, alpha ,beta, theta),[-20  20 0 20] ) 
title(['Fixed Theta Plot: theta = ', num2str(theta) ])
xlabel('alpha')
ylabel('beta')




% figure
% % fsurf(@(alpha,theta) b1(rho_1, alpha ,beta, theta) , [-20  20 0 1] ) 
% fsurf(@(theta,alpha) b1(rho_1, alpha ,beta, theta) , [0 1 -20 20] ) 
% title(['Fixed Beta Plot: beta = ', num2str(beta) ])
% xlabel('theta')
% ylabel('alpha')

figure
fsurf(@(alpha,theta) b1(rho_1, alpha ,beta, theta) , [-20  20 0 1] ) 
title(['Fixed Beta Plot: beta = ', num2str(beta) ])
xlabel('alpha')
ylabel('theta')




figure
fsurf(@(alpha,theta) b2(rho_1, alpha ,beta, theta),[-20  20 0 1] )
title(['Fixed Beta Plot: beta = ', num2str(beta) ])
xlabel('alpha')
ylabel('theta')


figure
fsurf(@(beta,theta) b1(rho_1, alpha ,beta, theta) , [0.1  20 0 1] ) 
title(['Fixed Beta Plot: alpha = ', num2str(alpha) ])
xlabel('beta')
ylabel('theta')




figure
fsurf(@(beta,theta) b2(rho_1, alpha ,beta, theta),[0.1  20 0 1] )
title(['Fixed Beta Plot: beta = ', num2str(alpha) ])
xlabel('beta')
ylabel('theta')

