clc
close all
clear all
N = 512;
%%%%%%%%%%%%%%%%Parameters Settings
g = 5;
J0 = 1/g;
a = 0.4;
tau = 3;
tau_v = 144;
dt = tau/10;
mbar = 151.2;
m = (mbar).*tau/tau_v;
density = N / (2*pi);
k = 5; %global inhibition factor
alpha = 0.19; %Input strength: Oscillatory tracking
v= 4.36/3*1e-3; %Moving speed of external input
T = 1.3*pi/v; %Simulation time
loc = -pi*3/4; %Initial position of the input

%%%%%%%%%%%%%%%%%%Matrix Construction
J = zeros(N, 1);
Iext = 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';
centerx_U=zeros(1,length(T/dt));
centerx_I=zeros(1,length(T/dt));
U = zeros(N, 1);
V = zeros(N, 1);
r = zeros(N, 1);
r_t=[];
time=[];
exppos = exp(1i*pos);
%%Video initiation
myVideo = VideoWriter('burst.mp4','MPEG-4'); 
myVideo.FrameRate = 40; 
open(myVideo); 
% pic_num=1;
t = 0;
while t < T
    loc = loc + v * dt;%Input position
    if loc >= pi
        loc = loc - 2 * pi;
    end
    %Calculate input
    dis = min(abs(pos - loc), 2 * pi - abs(pos - loc));
    Iext = alpha.*exp(-(dis.^2/(4*a^2)));%External input
    Irec=  ifft(Jfft.*fft(r)); %Recurrent input 
    %Update network state
    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));
     %Video generation
   count=floor(t/dt) ;
   r_t=[r_t,r(N/2)];
   time=[time,t];
    if count/5 == floor(count/5)
        subplot(1,2,1)
    plot(pos,Iext.*1e2/2,'k','linewidth',2),hold on
    plot(pos,r.*1e3,'b','linewidth',2)
    plot(pos(N/2),r(N/2).*1e3,'r.','markersize',20),hold off
    title('Population activity')
    axis([-pi pi 0 15])
    xlabel('x')
    ylabel('firing rate')
    set(gcf,'unit','centimeters','position',[20,10,20,13])
    set(gca,'linewidth',3,'fontsize',15,'fontname','Cambria Math');
%     legend('I^{ext}(x,t)','r(x,t)','Probe neuron')
        subplot(1,2,2)
    plot(time.*1e-3,r_t.*1e3,'b','linewidth',2),hold on
    plot(time(end).*1e-3,r_t(end).*1e3,'r.','markersize',20),hold off
    title('Probe neuron activity')
    axis([0 T.*0.9.*1e-3 0 15])
    xlabel('time')
    ylabel('firing rate')
    set(gcf,'unit','centimeters','position',[10,10,30,10])
    set(gca,'linewidth',3,'fontsize',15,'fontname','Cambria Math');
    drawnow
    frame = getframe(gcf);
    im = frame2im(frame); 
    writeVideo(myVideo,im); 
    end
	t = t + dt;
end
close(myVideo)