coupleN.m

% function omegout = coupleN(node,gext)
% input variables:
%          node is a network structure
%          gext is a scalar for external coupling
% output variable omegout is the oscilation freq. of each oscillator
% Oscillators are Poincare oscillators


function [omegout,y] = coupleN(node,gext)

[dum,N] = size(node);
mnomega = 1.0;

rndphase = 2*pi*rand(1,N);

for nodeloop = 1:N
    omega(nodeloop) = node(nodeloop).element;
end

%keyboard

y0 = rndphase;

% Settle-down loop
tspan = [0 1];  %800
[t,y] = ode45(@f5,tspan,y0);
[sy,dum] = size(y);

%keyboard

y1 = y(sy,:);
clear y

% Calculate the trajectory
tspan = [0 2000];
options = odeset('RelTol',1e-4*std(omega));
[t,y] = ode45(@f5,tspan,y1,options);
[sy,~] = size(y);

% Fit the frequency
for omloop = 1:N
    
    if rem(sy,4) == 0
        [mtmp(1),btmp(1)] = linfit(t(1:sy/2),y(1:sy/2,omloop),0);
        [mtmp(2),btmp(2)] = linfit(t(sy/2+1:sy),y(sy/2+1:sy,omloop),0);
        [mtmp(3),btmp(3)] = linfit(t(sy/4+1:3*sy/4),y(sy/4+1:3*sy/4,omloop),0);
        [mtmp(4),btmp(4)] = linfit(t,y(:,omloop),0);
    else
        sytmp = 4*floor(sy/4);
        [mtmp(1),btmp(1)] = linfit(t(1:sytmp/2),y(1:sytmp/2,omloop),0);
        [mtmp(2),btmp(2)] = linfit(t(sytmp/2+1:sytmp),y(sytmp/2+1:sytmp,omloop),0);
        [mtmp(3),btmp(3)] = linfit(t(sytmp/4+1:3*sytmp/4),y(sytmp/4+1:3*sytmp/4,omloop),0);
        [mtmp(4),btmp(4)] = linfit(t(1:sytmp),y(1:sytmp,omloop),0);
    end
    
    %[m(omloop),b(omloop)] = linfit(t,y(:,omloop),0);  
    
    m(omloop) = mean(mtmp);
    %m(omloop) = median(mtmp);
    
    w(omloop) = mnomega + m(omloop);
%     disp(strcat('Freq',num2str(omloop),' = ',num2str(m(omloop))))
%     figure(99)
%     plot(t,y(:,omloop))
%     pause(0.1)
end

omegout = m;

    function yd = f5(t,y)

        for omloop = 1:N
            temp = omega(omloop);
            linksz = node(omloop).numlink;
            for cloop = 1:linksz
                cindex = node(omloop).link(cloop);
                g = node(omloop).coupling(cloop);
                %temp = temp + (g/linksz)*sin(y(cindex)-y(omloop)) - gext*sin(y(omloop));   % Normalize by number of links in summation
                temp = temp + g*sin(y(cindex)-y(omloop)) - gext*sin(y(omloop));
            end
            yp(omloop) = temp;
        end

        for omloop = 1:N
            yd(omloop,1) = yp(omloop);
        end

    end     % end f5


end % end vandpol