pkg load statistics;

function y= get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= ones(N,1);
end

N = 200;
x0=0; y0=0;
x = zeros(2,N);

for i=1:1:N
    l = get_stride(1);
    th = get_angle(1);
    deltax = l*cos(th); deltay = l*sin(th);
    x(1,i) = deltax;
    x(2,i) = deltay;
end    
xinc = x(1,:); yinc = x(2,:);
xtraj = cumsum(xinc); ytraj = cumsum(yinc);
plot(xtraj, ytraj,"linewidth",4);

function y= get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= ones(N,1);
end

N = 200;
Np = 10;

for j=1:1:Np
    x0=0; y0=0;
    x = zeros(2,N);

    for i=1:1:N
        l = get_stride(1);
        th = get_angle(1);
        deltax = l*cos(th); deltay = l*sin(th);
        x(1,i) = deltax;
        x(2,i) = deltay;
    end
    xinc = x(1,:); yinc = x(2,:);
    xtraj = cumsum(xinc); ytraj = cumsum(yinc);
    plot(xtraj, ytraj,"linewidth",4);
    hold on;
end    
hold off;

function y= get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= ones(N,1);
end


function y=get_final_loc(N)
    x0=0; y0=0;
    x = zeros(2,N);

    for i=1:1:N
        l = get_stride(1);
        th = get_angle(1);
        deltax = l*cos(th); deltay = l*sin(th);
        x(1,i) = deltax;
        x(2,i) = deltay;
    end
    xinc = x(1,:); yinc = x(2,:);
    xtraj = sum(xinc); ytraj = sum(yinc); rorigin = (xtraj.^2 + ytraj.^2).^0.5;
    y= rorigin;
end

N = 200;
Np = 1000;
r = zeros(Np,1);

for j=1:1:Np
    r(j) = get_final_loc(N);
end

[counts, bins]=hist(r, 40);
bar(bins, counts/trapz(bins,counts));

s1= get_angle(10);
s2=get_stride(10);
dx = s2.*cos(s1); dy = s2.*sin(s1);
lx = sum(dx); ly = sum(dy);
disp(lx); disp(ly);

-2.3314
-1.5183

function y=get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= abs(cauchy_rnd(0, 1,[N,1]));
end

N = 200;
x0=0; y0=0;
x = zeros(2,N);

for i=1:1:N
    l = get_stride(1);
    th = get_angle(1);
    deltax = l*cos(th); deltay = l*sin(th);
    x(1,i) = deltax;
    x(2,i) = deltay;
end
xinc = x(1,:); yinc = x(2,:);
xtraj = cumsum(xinc); ytraj = cumsum(yinc);
plot(xtraj, ytraj);

function y=get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= abs(cauchy_rnd(0, 1,[N,1]));
end

N = 200;
Np = 10;

for j=1:1:Np
    x0=0; y0=0;
    x = zeros(2,N);
    for i=1:1:N
         l = get_stride(1);
        th = get_angle(1);
        deltax = l*cos(th); deltay = l*sin(th);
        x(1,i) = deltax;
        x(2,i) = deltay;
    end

    xinc = x(1,:); yinc = x(2,:);
    xtraj = cumsum(xinc); ytraj = cumsum(yinc);
    plot(xtraj, ytraj);
    hold on;
end

daspect([1 1 1]);

r = linspace(0, 10); f = 1./(pi^2*r.*(1+r.^2));
loglog(r, f);

warning: axis: omitting non-positive data in log plot
warning: called from
    __plt__>__plt2vv__ at line 495 column 10
    __plt__>__plt2__ at line 242 column 14
    __plt__ at line 107 column 18
    loglog at line 60 column 10

function y=get_angle(N) # N is the time steps
    y= random("uniform",0, 2*pi, [N,1]);
end

function y= get_stride(N)
    y= abs(cauchy_rnd(0, 1,[N,1]));
end

function y=get_final_loc(N)
    l = get_stride(N);
    th = get_angle(N);
    dx = l.*cos(th); dy = l.*sin(th);
    xtraj = sum(dx); ytraj = sum(dy); rorigin = (xtraj.^2 + ytraj.^2).^0.5;
    y= rorigin;
end


N = 200;
Np = 10000;
r = zeros(Np,1);

for j=1:1:Np
    r(j) = get_final_loc(N);
end

binsp = linspace(0, 1000);
[counts, bins]=hist(r, binsp);
bar(bins, counts/trapz(bins,counts));