% tristable_log.m
% for alpha12 = alpha13, alpha21 = alpha23, alpha31 = alpha32
alpha1 = 10.^([1:3]);%transcription plus translation rate of pLac construct.  need to have same # of each alpha
alpha2 = 10.^([1:3]); %transcription plus translation rate of pTet construct
alpha3 = 10.^([1:3]); %transcription plus translation rate of pBAD construct
beta1 = 2; %Cooperativity ie when there are more repressors repression goes up by more than the increase in repressors
beta2 = 2;
beta3 = 2;
beta4 = 2;
beta5 = 2;
beta6 = 2;
J = zeros(3,3);%creates empty mat for Jacobian
digits(3);
count=0;
isstable = zeros(numel(alpha1), numel(alpha2), numel(alpha3));
%a one will replace each indece which stable (indeces do not
%necessarily correspond to alpha values)

% plotting parameters
mksize = 10;%size of marker dots in plot
x = 0;
for a1 = 1:length(alpha1)
    for a2 = 1:length(alpha2)
        count
        for a3 = 1:length(alpha3)
            count = count +1;
            syms u v w; %creates the variables
            f = vpa(alpha1(a1)/(1+v^beta1)+alpha1(a1)/(1+w^beta2)-u);
            g = vpa(alpha2(a2)/(1+u^beta3)+alpha2(a2)/(1+w^beta4)-v);
            h = vpa(alpha3(a3)/(1+u^beta5)+alpha3(a3)/(1+v^beta6)-w);
            [x,y,z] = solve(f,g,h,u,v,w); %solves eqns f, g and h at the roots for u, v and w using vpa standard # of digits
            sizex = size(x);
            eq = 0;
            for ss = 1:sizex(1)%takes first dimension of size x (always 21x1)
                uu = double(x(ss)); %must use double type otherwise isreal() always returns 0 
                vv = double(y(ss));
                ww = double(z(ss));
                if isreal(uu) && isreal(vv) && isreal(ww)
                    %compute jacobian
                    J(1,1) = -1;%du/dudt of alpha1/(v^beta1 +1)+alpha2/(z^beta2+1) -u
                    J(1,2) = -alpha1(a1)*beta1*vv^(beta1-1)/((1+vv^beta1)^2);
                    J(1,3) = -alpha1(a1)*beta2*ww^(beta2-1)/((1+ww^beta2)^2);
                    J(2,1) = -alpha2(a2)*beta3*uu^(beta3-1)/((1+uu^beta3)^2);
                    J(2,2) = -1;
                    J(2,3) = -alpha2(a2)*beta4*ww^(beta4-1)/((1+ww^beta4)^2);
                    J(3,1) = -alpha3(a3)*beta5*uu^(beta5-1)/((1+uu^beta5)^2);
                    J(3,2) = -alpha3(a3)*beta6*vv^(beta6-1)/((1+vv^beta6)^2);
                    J(3,3) = -1;
                    if isnan(J)==0 %if there are no NaN numbers in J
                        ev = eig(J);
                        if ev(1)<0 && ev(2)<0 && ev(3)<0 %if eigen value is negative then that set is stable, positive values are unstable
                            eq = eq + 1;
                        end 
                    end
                end    
            end
%             a1a1 = log10(alpha1(a1))
%             a1a2=log10(alpha1(a2))
%             a1a3=log10(alpha1(a3))
%             a2a1=log10(alpha2(a1))
%             a2a2=log10(alpha2(a2))
%             a2a3=log10(alpha2(a3))
%             a3a1=log10(alpha3(a1))
%             a3a2=log10(alpha3(a2))
%             a3a3=log10(alpha3(a3))

            if eq >= 3 %need three stable solutions for tristable switch
            isstable(a1,a2,a3) = 1;% puts 1 where switch is stable
            plot3(log10(alpha1(a1)),log10(alpha2(a2)),log10(alpha3(a3)),'b.','markersize',mksize +10)
            hold on
%             plot3(log10(alpha1(a1)),log10(alpha2(a3)),log10(alpha3(a2)),'b.','markersize',mksize +2)
%             hold on
%             plot3(log10(alpha1(a2)),log10(alpha2(a1)),log10(alpha3(a3)),'b.','markersize',mksize +2)
%             hold on
%             plot3(log10(alpha1(a2)),log10(alpha2(a3)),log10(alpha3(a1)),'b.','markersize',mksize +2)
%             hold on
%             plot3(log10(alpha1(a3)),log10(alpha2(a1)),log10(alpha3(a2)),'b.','markersize',mksize +2)
%             hold on
%             plot3(log10(alpha1(a3)),log10(alpha2(a2)),log10(alpha3(a1)),'b.','markersize',mksize +2)
%             hold on
            else
            plot3(log10(alpha1(a1)),log10(alpha2(a2)),log10(alpha3(a3)),'ro','markersize',mksize-4)
            hold on
%             plot3(log10(alpha1(a1)),log10(alpha2(a3)),log10(alpha3(a2)),'r.','markersize',mksize)
%             hold on
%             plot3(log10(alpha1(a2)),log10(alpha2(a1)),log10(alpha3(a3)),'r.','markersize',mksize)
%             hold on
%             plot3(log10(alpha1(a2)),log10(alpha2(a3)),log10(alpha3(a1)),'r.','markersize',mksize)
%             hold on
%             plot3(log10(alpha1(a3)),log10(alpha2(a1)),log10(alpha3(a2)),'r.','markersize',mksize)
%             hold on
%             plot3(log10(alpha1(a3)),log10(alpha2(a2)),log10(alpha3(a1)),'r.','markersize',mksize)
%             hold on

            end
            xlabel('alpha1 (log10)')%annotates plot
            ylabel('alpha2 (log10)')
            zlabel('alpha3 (log10)')
            title('Plots Region of Stable Alpha Values (production rate) in log10', 'fontsize', 14)
            legend('stable points', 'unstable points');
        end
    end
end
csvwrite('csv_p7to8.dat', isstable); % writes to csv file, 
isstable
count