clear;

DMAX = 1000;

VMAX = 10000;

L = 100;

r = 1;

lambda = 0.1;

phi = (-pi/2) + lambda;

d_phi = (pi - (2 * lambda)) / DMAX;

dv = lambda / VMAX;

for i = 1:DMAX;

d_out = ((r * sin(phi))^2 + (L - r * cos(phi))^2)^0.5;

d_back = ((r * sin(phi + lambda))^2 + (L - r * cos(phi + lambda))^2)^0.5;

t_roundtrip = (d_out + d_back) / 2.9979e8;

velocity = (lambda * L) / t_roundtrip;

vc(i) = velocity / 2.9979e8;

phi = phi + d_phi;

end

for i = 1:VMAX;

v(i) = 0;

x(i) = (i-1)*dv;

for j = 1:DMAX;

if vc(j) > x(i);

if vc(j) <= (x(i)+dv)

v(i) = v(i) + (1/DMAX);

end

end

end

end

plot(x,v);

axis([0.0495 0.051 0 0.1]);