clc
clear all
close all
n1=200;
n2=200;
alpha=linspace(0.14,0.35,n1); %Range of input strength
m=linspace(80,200,n2)./48;%Range of feedback inh. stre.


tic
omega = zeros(n1,n2);
for i=1:n1
    parfor j=1:n2
        [omega(i,j),~] = omega_simulate(alpha(i),m(j));
    end
    disp(i/n1);
end
toc
save omega_final.mat

%%%%%%%%%%%%%%%%%%fig_3B
clc
clear all
close all
load omega_final.mat
[M,Alpha]=meshgrid(m,alpha);
surf(M,Alpha,omega),
colorbar
view([0 0 1])
set(gcf,'unit','centimeters','position',[25,17,12,8])
axis([min(m) max(m) min(alpha) max(alpha)])
xlabel('feedback inh. stre. m','fontname','Arial')
ylabel('input stre. \alpha','fontname','Arial')
colorbar
box on
set(gca,'linewidth',3,'fontsize',15,'fontname','Arial');

%%%%%%%%%%%%%%%%%%fig_3C: left panel

a = 0.4;
J0 = 1;
tau = 3;
tau_v = 144;
k = 5;% global inhibition factor
idx = 93; %index of the probe feedback inh. stre. m
m_probe = m(idx);

%%%Calculate the theoretical soluation of oscillaiton frequency vs. m
A_u = zeros(1,n1);
w = zeros(1,n1);
for i=1:n1
    A_u(i)=(J0+2*sqrt(pi)*a*k.*alpha(i))./(2*sqrt(pi)*a*k.*(1+m_probe));
    w(i)=sqrt(alpha(i)./(tau*tau_v*A_u(i)))/2/pi*1e3;
    if omega(i,idx)==0
        w(i)=0;
    end
end
figure
start_p = 30;
end_p = 150;
plot(alpha(start_p:end_p),w(start_p:end_p),'b','linewidth',3), hold on %analytical results
plot(alpha(start_p:10:end_p),omega(start_p:10:end_p,idx),'r.','markersize',30) %simulation results
axis([min(alpha(start_p:end_p))  max(alpha(start_p:end_p)) 9 14])
set(gca,'XTick',0.2:0.05:0.5);
set(gca,'YTick',8:2:16);
legend('Analytical results','Simulation results')
legend box off
legend location best
xlabel('input stre. \alpha','fontname','Arial')
ylabel('Osc. freq. \omega','fontname','Arial')
box off
set(gca,'linewidth',3,'fontsize',15,'fontname','Arial');
set(gcf,'unit','centimeters','position',[25,13,10,12])

%%%%%%%%%%%%%%%%%%fig_3C: right panel
figure
idx = 75; %index of the probe alpha
alpha_probe = alpha(idx);
%%%Calculate the theoretical soluation of oscillaiton frequency vs. alpha 
A_u = zeros(1,n2);
w = zeros(1,n2);
for i=1:n2
    A_u(i)=(J0+2*sqrt(pi)*a*k.*alpha_probe)./(2*sqrt(pi)*a*k.*(1+m(i)));
    w(i)=sqrt(alpha_probe./(tau*tau_v*A_u(i)))/2/pi*1e3;
    if omega(idx,i)==0
        w(i)=0;
    end
end

start_p = 35;
end_p = 150;
plot(m(start_p:end_p),w(start_p:end_p),'b','linewidth',3), hold on %analytical results
plot(m(start_p:end_p),omega(idx,start_p:end_p),'r.','markersize',30) %simulation results
axis([min(m(start_p:end_p))  max(m(start_p:end_p)) 9 14])
set(gca,'YTick',8:2:16);
legend('Analytical results','Simulation results')
legend box off
legend location best
xlabel('feedback inh. stre. m','fontname','Arial')
ylabel('Osc. freq. \omega','fontname','Arial')
box off
set(gca,'linewidth',3,'fontsize',15,'fontname','Arial');
set(gcf,'unit','centimeters','position',[25,13,10,12])