clc
clear all
close all
%%%%%%%%%%%%%%%%Parameters Settings
N = 512;
g = 5;
J0= 1/g;
a = 0.4;
tau = 3 ;
tau_v = 144;
dt = tau/10;
m = 150* tau/tau_v;
k = 5;% global inhibition factor
v = 4.36/3*1e-3;
T = 8*pi/v;
alpha=0.19;


%%%%%%%%%%%%%%%%%%Matrix Construction
J = zeros(N, 1);
Iext = zeros(N, 1);
U = zeros(N, 1);
V = zeros(N, 1);
r = 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));
		% 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);
%%%%%%%%%%%%%%
pos=pos';
loc=-pi*5/8;
centerx_U = zeros(1,length(T/dt));
centerx_I = zeros(1,length(T/dt));
r_t = zeros(N,length(T/dt));

j = 1;
t = 0;
while t < T
        loc = loc + v * dt;
    if loc >= pi
        loc = loc - 2 * pi;
    end
    dis = min(abs(pos - loc), 2 * pi - abs(pos - loc));
    Iext = alpha  * exp(-(dis.^2/(4*a^2)));
    Irec=ifft(Jfft.*fft(r));%;
    dU = dt * (-U - V + Iext+Irec)/tau;
	U = U + dU;
	dV = dt * (-V + m.*U) / tau_v;
	V = V + dV;
	U = max(U, 0);
	r = g.*U.^2./(1+k.*sum(U(:).^2));
    
    if t>T/4*3
        centerx_I(1,j)=loc;
        maxp = angle(exp(-1i.*pos)'*r);
        if loc > pi - 4*a && maxp < -pi+4*a
            maxp = maxp + 2 * pi;
        end

        if maxp > pi - 4*a && loc < -pi+4*a
            maxp = maxp - 2 * pi;
        end
        centerx_U(1,j)=maxp;
        r_t(:,j)=r(:);
        j=j+1;
    end
	t = t + dt;
    disp(t/T)
end
toc
time = dt:dt:(j-1)*dt;
[Time, Pos] = meshgrid(time,pos);
[CI, P] = meshgrid(centerx_I(1:length(time)),pos);
hold on

sigma = 5; %Only plot the most activated neurons
[r_max,pos_max] = max(r_t);
r_plot = zeros(N,length(time));
for i = 1:length(time)
    start_id = pos_max(i)-sigma;
    if start_id <= 0
        start_id = start_id+N;
    end
    end_id = pos_max(i)+sigma;
    if end_id > N
        end_id = end_id-N;
    end
    id = start_id:end_id;
    r_plot(id,i) = r_t(id,i);
end
relative_Pos = Pos-CI;
step_1 = 1;
step_2 = 5;
time_end = 3500;
time_s = 1800;
pos_end = 290;
pos_s = 110;
surf((Time(pos_s:step_1:pos_end,time_s:step_2:time_end))./1e3-0.55,...
    relative_Pos(pos_s:step_1:pos_end,time_s:step_2:time_end),...
    r_plot(pos_s:step_1:pos_end,time_s:step_2:time_end).*1e3)
shading flat
axis([0 0.4 -0.5 1])
view([0 0 1])
xlabel('Time(s)')
ylabel('Relative Loc.')
set(gca,'Xtick',0:0.1:0.5)
set(gca,'Ytick',-3:0.5:3)
colorbar
colormap hot
caxis([5 15])
box on
set(gcf,'unit','centimeters','position',[25,17,25,10])
set(gca,'linewidth',3,'fontsize',30,'fontname','Arial');