%%%%%%%%%%%%%%%%Smooth tracking
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 = 6 * tau/tau_v;
Alpha=0.1;
k = 0.7;% global inhibition factor
v = 4.36 / 3 *1e-3;
T = 8 * pi / v;
%%%%%%%%%%%%%%%%%%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));
j=1;
t = 0;
r_t = zeros(N,length(T/dt));
U_t = zeros(N,length(T/dt));
I_t = zeros(N,length(T/dt));

% m.*U

% tic
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));%.*(ratio)
    %if t>10
    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(:);
        U_t(:,j) = U(:);
        I_t(:,j) = Iext;
        Ar(j) = max(r);
        j=j+1;
    end
	t = t + dt;
%     disp(t/T)
end
% toc

time=1:90:length(centerx_U)/2-500;
L=centerx_U-centerx_I;

[Pos,Time] = meshgrid(pos,time);
c=zeros(size(Time));
factor = max(max(U_t))/1.2;
start_time = floor(length(centerx_U)/10+50);
for i=1:length(c(:,1,1))
    for j=1:length(c(1,:,1))
        c(i,j,1)=1-U_t(j,time(i)+start_time)/1.5/factor;
        c(i,j,2)=1-U_t(j,time(i)+start_time)/3/factor;
        c(i,j,3)=1-U_t(j,time(i)+start_time)/3/factor;
    end
end
surf(Pos,Time.*dt*1e-3,U_t(:,time+start_time)',c),hold on
c2=zeros(size(Time));
factor = max(max(I_t))/1.5;
for i=1:length(c2(:,1,1))
    for j=1:length(c2(1,:,1))
        c2(i,j,1)=0.95-I_t(j,time(i)+start_time)/3/factor;
        c2(i,j,2)=0.95-I_t(j,time(i)+start_time)/1.5/factor;
        c2(i,j,3)=0.95-I_t(j,time(i)+start_time)/1.5/factor;
    end
end
scale = max(max(U_t))/max(max(I_t));
surf(Pos,Time.*dt*1e-3,I_t(:,time+start_time)'.*scale,c2),hold on
shading flat
alpha(.7);
view([-0.05 -1 3])
set(gca,'XTick',-pi:pi:pi);
set(gca,'XTicklabel',{'-\pi','0','\pi'})
set(gca,'ZTick',0:1);
axis([-pi pi 0 max(time)*dt*1e-3])
plot3(centerx_U(time+start_time),time*dt*1e-3,max(U_t(:,time+start_time))+0.01,'b','linewidth',1)
plot3(centerx_U(time+start_time),time*dt*1e-3,max(U_t(:,time+start_time))+0.01,'b.','markersize',22)
plot3(centerx_I(time+start_time),time*dt*1e-3,max(I_t(:,time+start_time).*scale)+0.01,'color',[0.6,0,0],'linewidth',1)
plot3(centerx_I(time+start_time),time*dt*1e-3,max(I_t(:,time+start_time).*scale)+0.01,'.','color',[0.6,0,0],'markersize',22)
grid off
xlabel('Neural position x')
ylabel('Time')
zlabel('U(x,t)')
set(gcf,'unit','centimeters','position',[20,17,20,15])
box off
set(gca,'linewidth',2,'fontsize',25,'fontname','Arial');
