% try flight of a bumblebee

% idea constant speed with diffusing direction theta
N = 1000; % number of flyers
% initial position
c = zeros(N,2);
theta = 2*pi*rand(N,1);  % random initial direction
dt = 0.001;  %time step
D = 1; % specify the diffusion coefficient for theta
sc = sqrt(2*D*dt);
v = 1;  % specify the ballistic velocity
n_time_steps = 5000;  %specify the number of time steps
mn = [];
T = [];


for j = 1:n_time_steps % number of time steps
dtheta = sc*randn(N,1);
%take a forward time step
theta = theta + dtheta;
cnew = c + v* dt*[cos(theta),sin(theta)];
c = cnew;
rsq = c(:,1).^2+c(:,2).^2;
mn(j) = sum(rsq)/N;  % mean squared displacement
T(j) = j*dt;
end

figure(1)
plot(c(:,1),c(:,2),'*')
%figure(2)
%plot(csv(:,1).^2+csv(:,2).^2)
% do statistics on radial distibution  
rsq = c(:,1).^2+c(:,2).^2;
rstat = sort(rsq);

[p,s] = polyfit(rstat(1:N-1),log(1-[1:N-1]'/N),1)
y=polyval(p,rstat);
figure(2)
plot(rstat(1:N-1),log(1-[1:N-1]/N),rstat,y,'--')

figure(3)
plot(T,mn,'linewidth',2)
xlabel('Time','fontsize',20)
ylabel('Means Squared Displacement','fontsize',20)