ou_msd

clear all; clc;

ph1=0.08245;
ph2=0.3674;
beta1=0.008526;
beta2=0.3829;

gamma=0.441689378;
dt=0.001;
%countR=0;
%D=1;
D=0.078570873;

msd=zeros(1000);
X=zeros(1000);
Y=zeros(1000);
cou=zeros(1000);

 vx=0;
 vy=0;
 x=0;
 y=0;
 t=0;


for ksss=1:100



    for kjkj=1:1000

    for ssss=1:1000

           r1 = normrnd(0,1);
           r2 = normrnd(0,1);

           vx1=vx-gamma*vx*dt+sqrt(dt*D*2)*r1;
           vy1=vy-gamma*vy*dt+sqrt(dt*D*2)*r2;

           vx=vx1;
           vy=vy1;

           X(kjkj)=x+vx*dt; %подсчет координаты x
           Y(kjkj)=y+vy*dt; %подсчет координаты y

           %X(kjkj)=x;
           %Y(kjkj)=y;

           t=t+dt;
    end

    end

%подсчет MSD

for tt=1:1000
for tt0=1:1000

    dd=tt+tt0;

if (dd <1000)
msd(tt)=msd(tt)+(X(dd)-X(tt0))^2+(Y(dd)-Y(tt0))^2;
cou(tt)=cou(tt)+1;
end

end
end

    fid = fopen('OU_msd.dat','w');


    for ll=1:300
        tr=msd(ll)/cou(ll);
        analytic=4*D*(gamma*ll+exp(-gamma*ll)-1)/gamma^3;
        fprintf(fid, '%d;%.12f;%.12f \r\n',ll,tr,analytic);
        %fprintf(fid, '%d;%.12f \r\n',ll,analytic);
        end

    fclose(fid);

    ksss

end