%%%%%%%%%%Intrinsic speed of the network bump
clc
close all
clear all

%%%Parameter settings
tau = 3;
tau_v = 144;
M = 0:0.05:5;
M = M.*tau./tau_v;
N=512;
density = N / (2*pi);
k= 0.7;
a =0.4;
g = 5;
J0 = 1/g;
dt = tau / 10;
%%%%%%%%%%%%%%%%%%Matrix Construction
J = zeros(N, 1);
% Map all neurons to [-pi, pi)
x=linspace(-pi, pi, N+1);
pos = x(1: N);
for i = 1: N
		dx = min(pos(i)-pos(1), pi-pos(i)); %Periodical boundary
		% choose the weight matrix to be a 2D Gaussian Distribution
		J(i) = J0/(sqrt(2*pi)*a) * exp(-(dx^2)/(2*a^2));
end
Jfft = fft(J);
exppos = exp(1i*pos);
pos=pos';
v = 4.36 / 3 *1e-3;
v_int = zeros(1,length(M));
v_int_s = zeros(1,length(M));
T = 2*1e4;
centerx_U=zeros(1,length(T/dt));
centerx_I=zeros(1,length(T/dt));
for i = 1:length(M)
m = M(i);

j=1;
t = 0;
U = zeros(N, 1);
V = zeros(N, 1);
r = zeros(N, 1);
loc=-(pi/2+pi/8);

while t<T
    m=M(i);
    if t<10
        loc = loc + v * dt;
        if loc >= pi
            loc = loc - 2 * pi;
        end
        %loc
        loc_cal=loc;
        dis = min(abs(pos - loc_cal), 2 * pi - abs(pos - loc_cal));
        Iext = 0.5  * exp(-(dis.^2/(4*a^2)));
        else
            Iext=zeros(N,1);
    end
    dU = dt * (-U + ifft(Jfft.*fft(r(:))) - V + Iext)/tau;
	dV = dt * (-V + m*U) / tau_v;
	U = U + dU;
	V = V + dV;
	U = max(U, 0);
	r(:) = g.*U.^2./(1+k*sum(U(:).^2));
    centerx_U(1,j)=angle(exppos*r);
    j=j+1;
	t = t + dt;
end
    L = centerx_U;
    if m <= tau/tau_v || max(L)-min(L) < 0.1
        v_int_s(i) = 0;
        v_int(i) = 0;
    else
        v_int(i)=2*a/tau_v*sqrt(m*tau_v/tau-sqrt(m*tau_v/tau));
        [~,pos_max]=max(L);
        pos_max=pos_max(end);
        T_interval = pi/v_int(i);
        v_int_s(i)=(L(pos_max)-L(pos_max-floor(T_interval/dt)))/(T_interval);
    end
    disp(i/length(M));
end
%%
figure
% v_int = smoothdata(v_int_s);
plot(M,v_int_s.*1e3./3,'b','linewidth',1),hold on;

axis([0 max(M) -0.5 3])
set(gca,'YTick',0:1:10);
set(gca,'XTick',0:0.05:1);
xlabel('feedback inh. stre. m','fontname','Arial')
ylabel('v_{int}(1/s)','fontname','Arial')
set(gcf,'unit','centimeters','position',[25,17,20,15])
grid off
box on
set(gca,'linewidth',1,'fontsize',25,'fontname','Arial');