clc
clear all
close all
%%%%%%%%%%%%%%%%Parameters Settings
mbar = 110;%150
N = 512;
g = 5;
J0 = 1/g;
k = 5;
a =0.4;
tau = 3;
tau_v = 144;
dt = tau/10;
m = 110* tau/tau_v;
v = 4.36/3*1e-3;
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';
T=8*pi/v;
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;

% 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(:);

        j=j+1;
    end
	t = t + dt;
%     disp(t/T)
end
% toc

time=1:length(centerx_U);
L=centerx_U-centerx_I;
L_diff = L(time);
[~,Locs_diff] = findpeaks(L_diff);
period = mean(diff(time(Locs_diff)).*dt);
L_abs = abs(centerx_U(time));
LI_abs = abs(centerx_I(time));
theta_cycle = (time.*dt-time(Locs_diff(14)).*dt)./period;

%%%%%%%%fig_4B&C upper panel
figure
a_ht = abs(L_diff+centerx_I(time));
ht = centerx_U(time);
[~,Locs_ht] = findpeaks(-a_ht);
plot(theta_cycle,ht,'k','linewidth',3),hold on
plot(theta_cycle(Locs_ht),ht(Locs_ht),'r.','markersize',40)
plot([-11 11],[0 0],'k--','linewidth',2)

axis([-10.5 10 -2.5 2.5])
% axis([-2.5 2.5 -1 1]) %fig_4C
set(gcf,'unit','centimeters','position',[20,17,25,6])
box off
set(gca,'linewidth',3,'fontsize',15,'fontname','Arial');
xlabel('theta cycles')
ylabel('h(t)')

%%%%%%%%fig_4B lower panel & fig_4C middle panel
figure
[~,Locs_r] = findpeaks(r_t(N/2,time));
plot(theta_cycle,r_t(N/2,time).*1e3,'k','linewidth',3),hold on
plot(theta_cycle(Locs_r),r_t(N/2,time(Locs_r)).*1e3,'r.','markersize',40)
set(gca,'ytick',0:10:20);
axis([-10.5 10 0 1e-3*3.*5.*1e3])
% axis([-2.5 2.5 0 1e-3*3.*5.*1e3])
set(gcf,'unit','centimeters','position',[20,17,25,6])
box off
set(gca,'linewidth',3,'fontsize',15,'fontname','Arial');
% xlabel('v_{ext}')
xlabel('theta cycles')
ylabel('r(t)')

%%%%%%figure_4C lower panel
figure
T_cycle = (time.*dt-time(Locs_diff(1)).*dt)./(period);
phase_point = rem(T_cycle(Locs_r),1)*2*pi-pi*3/2;
phase_point(phase_point<-2.05) = phase_point(phase_point<-2.05)+2*pi;
phase_pt = phase_point(phase_point<2.4);
phase_pro = phase_point(phase_point>=2.4);
Tc = theta_cycle(Locs_r(phase_point<2.4));
Tc_pro = theta_cycle(Locs_r(phase_point>=2.4));
plot(Tc,phase_pt,'r.','markersize',40),hold on
plot(Tc_pro,phase_pro,'r*','markersize',20,'linewidth',3)
set(gca,'xtick',-15:20);
a=get(gca,'xticklabel');
% a = a(5:end-5);
b=cell(size(a));
b(mod(1:size(a,1),5)==1,:)=a(mod(1:length(a),5)==1,:);
set(gca,'xticklabel',b);
set(gca,'ytick',-3/2*pi:pi:2*pi);
set(gca,'YTickLabel',{'-3\pi/2','-\pi/2','\pi/2','3\pi/2'});
set(gcf,'unit','centimeters','position',[20,17,25,10])
box off
set(gca,'linewidth',3,'fontsize',30,'fontname','Arial');
% xlabel('v_{ext}')
xlabel('theta cycle')
ylabel('phase')
axis([-10.5 10 -pi 2*pi])
