%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Simulation of distributed sychronization
% Author: K. Passino
% Version: Sept. 1, 2001
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

global N K omega


N=10;	% The number of oscillators

% To study phase transitions:
omega=pi*rand(1,N); % Set natural frequencies of the oscillators
%g0=(1/(sqrt((2*pi)^N))); % The g0 function for a Gaussian
g0=1/pi
Kc=(2/(pi*g0))

epsilon=Kc*0.1; % So can add or subtract 10% of the value to see the phase transition
K=Kc-epsilon;	% The strength of coupling

% For now just use:
omega=ones(1,N); % Keep them all the same
K=1;

% Define simulation parameters:

Tfinal=20; % Units are seconds (keep it even)
%Tspan=[0 Tfinal];
Tspan=0:0.01:Tfinal+0.01;

% Define initial conditions:

theta0=2*pi*rand(1,N)   % Pick a uniform distribution of inital phase angles for the oscillators

[t,z]=ode45('distributed_sych_func',Tspan,[theta0]);
%[t,z]=ode15s('distributed_sych_func',Tspan,[theta0]);

theta=z;   % Oscillator states

thetabar=mean(theta');

% Convert for plotting


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the plant signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1) % Not really useful:
clf
subplot(211)
polar(theta(:,1),1+0*theta(:,1),'k-');
hold on
for i=1:N
	polar(theta(:,i),1+0*theta(:,i),'r-');
	hold on
end
hold off
xlabel('Position')
ylabel('Velocity')
title('Oscillator trajectory')
subplot(212)
plot(sin(theta(:,1)),cos(theta(:,1)),'k-')
xlabel('Position')
ylabel('Velocity')
title('Oscillator trajectory')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2) 
clf % To get insights into dynamics+Lyap. functions
subplot(211)
for i=1:N-1
	plot(t,theta(:,i)-theta(:,i+1),'k-')
	hold on
end
ylabel('Phase error')
title('Phase errors (black)')
subplot(212)
for k=1:length(t)
	V(k,1)=0;
	for i=1:N-1
		V(k,1)=V(k,1)+(0.5)*((theta(k,i)-theta(k,i+1)))*((theta(k,i)-theta(k,i+1)));
	end
end	
plot(t,V,'b-')
xlabel('Time, sec.')
ylabel('V')
title('Lyap. candidate (blue)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To get insights into dynamics
figure(3) 
clf
subplot(211)
for i=1:N-1
	plot(t,sin(theta(:,i)-theta(:,i+1)),'r-')
	hold on
end
title('Sin and sin^2 of phase errors')
ylabel('sin(diff)')
subplot(212)
for i=1:N-1
	plot(t,sin(theta(:,i)-theta(:,i+1)).^2,'g-.')
	hold on
end
hold off
xlabel('Time, sec.')
ylabel('sin^2(diff)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To get insights into stability analysis, plot potential Lyapunov functions
figure(4) 
clf

Vtot=0*ones(length(t),1);

% The following goes to zero, but can increase (sin of difference)
for i=1:N-1
	V(:,i)=(sin(theta(:,i)-theta(:,i+1)));
	Vtot=Vtot+V(:,i);
end


for i=1:N-1
	plot(t,V(:,i),'r-')
	hold on
end
plot(t,Vtot,'b-')
ylabel('sum of sin(diff) (blue) and sin of each diff (red) (but only in order of numbers)')
xlabel('Time,sec.')
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5) 
clf
clear V
% Try another (sum of sin of differences)
Vsum=0*ones(length(t),1);

for k=1:length(t);	
for i=1:N

	summ=0;
	for j=1:N % Create the sum for each oscillator that links the oscillators
		difftheta=theta(k,j)-theta(k,i); 
		summ=summ+sin(difftheta);
	end
	
	V(k,i)=summ; 

end
Vsum(k,1)=sum(V(k,i));
end

for i=1:N
	plot(t,V(:,i),'r-')
	hold on
end

plot(t,Vsum(:,1),'b--')
hold off
xlabel('Time, sec.')
ylabel('Phase error')
title('Indiv. sums (red) and total (blue)')



pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Try to do it simultaneously... interleave the points, but only for two of them
for i=1:2:length(t)
	xx(i)=sin(theta(i,1));
	xx(i+1)=sin(theta(i,2));
	yy(i)=cos(theta(i,1));
	yy(i+1)=cos(theta(i,2));
end 


figure(7) % This does it - - with a line between the phase angles, but with some kind of error
clf
plot(sin(theta(:,1)),cos(theta(:,1)),'k.')
hold on
xlabel('Position')
ylabel('Velocity')
title('Oscillator trajectory')

comet(xx,yy,10) % Ok, this is odd, but it was the only way I could get it to keep the comet plot for printing.

% Note that if you do not really need the motion of the comet trajectory you could simply take the points and draw a line
% between them for each iteration.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%