% Repressilator Model Solver.
% Call Equation funciton of Repressilator rep.m
% 6 variables
% Created by Kai, 2011/05/25
% modified by Bruce,  2011/06/08
% dual repressilator, opamp simulation
clc;
clear all;


t = [0 1000];% Simulation time

y0= [0.001 0 0 .001 0 0 ];% Initial Values 1-3 are rep1 --  4-6 are rep2

% Call ODE solver
[T0,Y0] = ode15s(@opampX2odePhaseSpace,t,y0);

figure(1); clf;
plot(T0,Y0(:,1),T0,Y0(:,2),T0,Y0(:,3));
legend('v1','v2','v3');
hold on
plot(T0,Y0(:,4),'.',T0,Y0(:,5),'.',T0,Y0(:,6),'.');
title('opamp ode15s')

figure(2); clf
% ignore the first third when plotting phase space
% to see all, set start=1
start = fix(length(Y0(:,1))/3) ;
%start=1;
plot3(Y0(start:end,1),Y0(start:end,2),Y0(start:end,3));
hold on
plot3(Y0(start:end,4),Y0(start:end,5),Y0(start:end,6),'r');
xlabel('v1')
ylabel('v2')
zlabel('v3')
box on
legend('rep1','rep2');
%set(gca, 'Projection', 'perspective')
set(gcf, 'Renderer', 'opengl')
axis vis3d

% return
% % parallel axes projection -- not too useful in this case
% figure(3); clf;
% box on
% %view(22,63)
% axis vis3d
% sh=surf(Y0);
% box on
% %view(22,63)
% axis vis3d
% % xs=1:3;
% % ys=0:dt:tmax;
% % zs=pos';
% % sh=surf(xs,ys,zs)
% set(sh, 'edgecolor','none')
%  colormap(hot)
% % set(gca, 'zlim', [-20 50])
% % set(gca, 'xtick', [1 2 3])
% % set(gca, 'xticklabel', {'x1','x2','x3'})
% % xlabel('variable');ylabel('time');zlabel('amp');
% % rotate3d on