syms f(alpha,beta,theta)
syms b1(alpha,beta,theta)
syms b2(alpha,beta,theta)
syms q1(beta,theta)
syms q2(beta,theta)

mu_1 = sym(1);
rho_1 = sym(1);

%TEST (new)
b1(alpha,beta,theta) = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
b2(alpha,beta,theta) = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));

q1(beta,theta) = mu_1.*(beta./(theta+(1-theta).*beta)) ;
q2(beta,theta) = mu_1.*((1-theta)+theta.*beta); 



fprintf('Functions to be minimized')
% f(beta,theta,alpha) = (mu_1.*(beta./(theta+(1-theta).*beta))).*((3*rho_1/2).*beta.*(1-theta.*(1+alpha))).^2 - (mu_1.*((1-theta)+theta.*beta)).*(3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha)).^2;
f(beta,theta,alpha) = q1(beta,theta).*b1(alpha,beta,theta).^2 - (q2(beta,theta).*b2(alpha,beta,theta).^2) ;

% Fix Theta ... 
thetav = sym(1/2);
f = subs(f,theta,thetav);
b1 = subs(b1,theta,thetav);
b2 = subs(b2,theta,thetav);
q1 = subs(q1,theta,thetav);
q2 = subs(q2,theta,thetav);



% f = subs(f,theta,0.25);
% b1 = subs(b1,theta,0.25);
% b2 = subs(b2,theta,0.25);
% q1 = subs(q1,theta,0.25);
% q2 = subs(q2,theta,0.25);



% f = subs(f,alpha,3);


eq = f == 0;


S = solve(eq,alpha, beta,'MaxDegree' , 5)
% S = solve(eq,alpha, beta,theta,'MaxDegree' , 5)

S.alpha
S.beta


% CHECK
b1 = subs(b1,alpha,S.alpha(1));
b1 = subs(b1,beta,S.beta(1));
b2 = subs(b2,alpha,S.alpha(1));
b2 = subs(b2,beta,S.beta(1));

q1 = subs(q1,beta,S.beta(1));
q2 = subs(q2,beta,S.beta(1));

b1 = double(b1)
b2 = double(b2)
q1 = double(q1)
q2 = double(q2)


b1.*((q1).^2)- ( b2.*(q2.^2)) 



f = subs(f,alpha,S.alpha(1));
f = subs(f,beta,S.beta(1))
simplify(f)


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



fprintf('--------------------')