function [omega,type]=omega_simulate(alpha, m)%, filename

%%%%Parameter settings
type = 0;
omega = 0;

N = 512; %Neuron number
a =0.4; %Connection width
g = 5; %Gain factor
J0 = 1/g;%Connection strength
tau = 3; %time constant of U, 3ms
tau_v = 144; %time constant of V, 144ms
dt = tau / 10; %simulation time step
k = 5; %global inhibition strength
J = zeros(N, 1); %Connection matrix
x=linspace(-pi, pi, N+1); 
pos = x(1: N);%Preferred location of place cells
for i = 1: N
		dx = min(pos(i)-pos(1), pi-pos(i)); %Periodic condition
		% 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';
v = 4.36 / 3 *1e-3; %Moving speed of external input
T = 27000; %Simulation time

centerx_U=zeros(1,length(T/dt)); %Mass center position of U
centerx_I=zeros(1,length(T/dt)); %Mass center position of U

%%%%Initialization
j=1;
t = 0;
U = zeros(N, 1);
V = zeros(N, 1);
r = zeros(N, 1);
loc=-(pi/2+pi/8);

while t < T
    
    %%%Update input
        loc = loc + v * dt;
    if loc >= pi
        loc = loc - 2 * pi; 
    end
    dis = min(abs(pos - loc), 2 * pi - abs(pos - loc));%Periodic condition
    Iext = alpha  * exp(-(dis.^2/(4*a^2))); %External input
    centerx_I(1,j)=loc;
   
    %Update network state
    dU = dt * (-U + ifft(Jfft.*fft(r(:))) - V + Iext) / tau;
	dV = dt * (-V + m*U) / tau_v;
	U = U + dU;
	V = V + dV;
	U = max(U, 0);
	r(:) = g.*U.^2./(1+k*sum(U(:).^2));
    
    %Calculate bump center position
    maxp = angle(r'*exp(1i*pos)); %Center mass of U in polar coordinates
    %Eliminating the boundary effects
    if loc > pi - 2*a && maxp < -pi+2*a
        maxp = maxp + 2 * pi; 
    end
    
    if maxp > pi - 2*a && loc < -pi+2*a
        maxp = maxp - 2 * pi;
    end
    centerx_U(1,j)=maxp;
    j=j+1;
	t = t + dt;
end

checkpoint = 4500/dt; %A time period for the network state to be stable
L=centerx_U-centerx_I; % Relative positions of the U bump with respect to input location
data = L(checkpoint:end)-mean(L(checkpoint:end)); % Normalized oscillation curve





if max(data)-min(data) > 1 %The bump gets too far away from input
    type = 4; %Traveling wave state
    
end    
% l_search = 180;
% if type == 0 %The network is not at traveling wave state
%     for i = checkpoint:(T - 200)/dt
%         if max(L(i-l_search/dt:i+l_search/dt))-min(L(i-l_search/dt:i+l_search/dt))>0.1
%             type = 3;
%             omega = abs(fshift(max_o));
%         end
%     end
% end
if type == 0 %The network is not at traveling wave state
    if max(data)-min(data)>0.1 %The offset between U and Iext is in reasonable range
        type = 3; %Oscillatory tracking state
        %%%Calculate the frequency by doing Fourier transform
        fL = fft(data);
        fL_shift = fftshift(fL);    
        n=length(data);
        power = abs(fL_shift).^2/n; 
        fs = 1/dt*1e3; %The unit of frequency is Hertz
        fshift = (-n/2:n/2-1)*(fs/n);
        [~,max_o] = max(power);
        omega = abs(fshift(max_o));
    end
end
if type == 0 %The network is at smooth tracking state
        if sum(L(checkpoint:length(L))) < 0 
            type = 1; %Delayed tracking state
        else 
            type = 2; %Anticipative tracking state
        end
end


end

