syms alpha beta theta rho_1 mu_1 
assume(beta,'positive')
assume(theta,'positive')
assume(mu_1,'positive')
assume(rho_1,'positive')

% 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,rho_1) = (3*rho_1/2).*beta.*(1-theta.*(1+alpha));
b2(alpha,beta,theta,rho_1) = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha));




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

% 
% 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/4);
% q3_star = subs(q3_star,theta,thetav);
% b1 = subs(b1,theta,thetav);
% b2 = subs(b2,theta,thetav);
% q1 = subs(q1,theta,thetav);
% q2 = subs(q2,theta,thetav);


% theta = thetav
eq1 = (3*rho_1/2).*beta.*(1-theta.*(1+alpha)) == (q1(beta,theta,mu_1)*q2(beta,theta,mu_1))^(1/2)

eq2 = (3*rho_1/(4*((1-theta)+theta.*beta))).*(1-theta.*(1+beta.*alpha)) == mu_1.*(beta./(theta+(1-theta).*beta)) 

eqns = [eq1, eq2] 

% eqns = subs(eqns,{theta,alpha,rho_1},{1/20,2,1} ) % works! 
% eqns = subs(eqns,{rho_1,alpha,theta},{1,2,1/4} )
eqns = subs(eqns,{theta},{1/4} )                        % fix only theta...
% eqns = subs(eqns,{theta,mu_1,rho_1},{1/4,10,1} )


% eqns = subs(eqns,{mu_1,beta},{387/1250,7741/10000} )
% simplify(eqns)

% eqns = subs(eqns,{theta,mu_1,rho_1},{1/4,10,5} )
% eqns = subs(eqns,beta,2.0)

% S = solve(eqns,mu_1,beta,'MaxDegree' , 5)
% S = solve(eqns,theta,beta,'MaxDegree' , 5)
% S = solve(eqns,theta,mu_1,'MaxDegree' , 5)


%interesting:
% eqns = subs(eqns,{theta,alpha,rho_1},{1/20,2,1} )
% S = solve(eqns,mu_1,beta,'MaxDegree' , 5)


% S.mu_1
% S.beta
% double(S.mu_1)
% double(S.beta)
% 
% 
% eqns = subs(eqns,{mu_1,beta},{0.3096,0.7741} )
% simplify(eqns)


S = solve(eqns, alpha,beta, rho_1,mu_1, 'MaxDegree' , 5)
% S = solve(eqns, mu_1,beta, 'MaxDegree' , 5)

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

S.alpha
S.beta
double(S.alpha)
double(S.beta)

eqns = subs(eqns,{alpha,beta},{95.2970,0.1577} )
% simplify(eqns)


% 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);




% 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('--------------------')