%two opamp repressilators coupled

clc;
clear all;

dt = 0.02 ;
t = 0:dt:1000 ;% Simulation time

%set up output array for 3 voltages
v1 = zeros(length(t),3) ;
v1(1,1:3) = [ 0.1 0 0];
v2 = zeros(length(t),3) ;
v2(1,1:3) = [ 0.1 0 0];

% model parameters
tau1 = 2; %Time constant in Sec of RC
tau2 = 3; %Time constant in Sec of RC
%k = tau/dt ;% constant from differential equation
g = -12*[1 1 1] ; % gain of opamp inverter

% coupling strength
% For gain -12 and tau1=2 and tau2=3
%   1:1 at beta=0.10 
%   period 6 at beta=0.0995
%   period 1 at beta=0.0996
%   period 2 at beta .01 
%   periodic (40 cycles?) at beta .04 
%   periodic with  period (4) at beta=0.09
%   approx periodic with  period (8) at beta=0.08?
%   longer period with beta=0.2 
% There is a STRONG interaction between gain, beta and tau2
%beta = .09957315; % period 1
%beta = .09957314; period 8

beta = 0.105

% For gain -12 and tau1=2 and tau2=3
% beta = 0.1: phase space show period 1
% beta = 0.105: phase space show period=1 for v1, period=2 for v2
% beta = 0.11: phase space show period=1 for v1, continum for v2
% beta = 0.12: phase space show period=2 for v1, period=8? for v2
% beta = 0.15: phase space shows chaos?
% beta = 0.24: phase space show period=2 for v1, period=4 for v2
% beta = 0.3: phase space show period=2 for v1, period=2 for v2
% beta = 0.4: phase space show period=1 for v1, period=1 for v2

% negative betas?
%beta = -.05 ;

%solve by Euler iteration
for i=2:length(t)
    % output of opamps
    va1 = g .* [v1(i-1,3)+beta*v2(i-1,3), v1(i-1,1), v1(i-1,2)];
    va1(va1>10) = 10;
    va1(va1<-1) = -1;
    % change of voltage on capacitor
    v1(i,1:3) = v1(i-1,1:3)+ dt/tau1*(va1 - v1(i-1,1:3));
   
     va2 = g .* [v2(i-1,3)+beta*v1(i-1,3), v2(i-1,1), v2(i-1,2)];
    va2(va2>10) = 10;
    va2(va2<-1) = -1;
    % change of voltage on capacitor
    v2(i,1:3) = v2(i-1,1:3)+ dt/tau2*(va2 - v2(i-1,1:3));
end

figure(1); clf;
plot(t,v1(:,1),t,v1(:,2),t,v1(:,3));
hold on
plot(t,v2(:,1),'.',t,v2(:,2),'.',t,v2(:,3),'.');
legend('v1','v2','v3');
title('opamp Euler')

figure(2); clf
% only use the last 1/2 of the run for phase space
start = fix(length(t)/2);
plot3(v1(start:end,1),v1(start:end,2),v1(start:end,3));
hold on
plot3(v2(start:end,1),v2(start:end,2),v2(start:end,3),'r');

xlabel('v1')
ylabel('v2')
zlabel('v3')
box on
%set(gca, 'Projection', 'perspective')
%set(gcf, 'Renderer', 'opengl')
axis vis3d




