% gjenbruk figure 1 fra finn_spor!
figure(1), hold on

tr = trackmatrix_new;
numpart=max(tr(:,4));
trny=[];
trnum = [];

%kontrollering av spor for å fjerne de partiklene som sitter fast
%loop gjennom alle partiklene for å kontrollere sporene.
figure(2), clf
hold on
m=1;
for n=1:numpart 
    I=find(tr(:,4)==n);
%     %de tre linjene nedenfor gir muligheten til å sjekke spor manuelt
%     figure(1), plot(tr(I,1),tr(I,2))
%      godta=input('Godta/Forkast sporet? g/f: ','s');
%     if godta=='g'
%         trny=[trny;tr(I,:)];
%         trnum(m)=n;
%         m=m+1;
%         figure(2),plot(tr(I,1),tr(I,2),'g') %plot godkjente spor i grønt
%     else
%         figure(2)
%         plot(tr(I,1),tr(I,2),'r') % og ikke godkjente i rødt
%     end
    % med de neste tre linjene tar man bare med partikler som har beveget
    % seg en viss lengde i x- og y-retning (alternativ til den manuelle
    % kontrollen)
    grense = 5; %sett inn et passende kriterium her, for å fjerne partikler som sitter fast
    ddxx = max(tr(I,1))-min(tr(I,1));
    ddyy = max(tr(I,2))-min(tr(I,2));
    if ddxx > grense && ddyy > grense
        figure(1), plot(tr(I,1),tr(I,2))
        figure(2), plot(tr(I,1),tr(I,2))
        trny=[trny;tr(I,:)];
        trnum(m)=n;
        m=m+1;
        figure(1),plot(tr(I,1),tr(I,2),'g') %plot godkjente spor i grønt
        figure(2),plot(tr(I,1),tr(I,2),'g') %plot godkjente spor i grønt
    else 
        figure(1),plot(tr(I,1),tr(I,2),'r.','MarkerSize',5) % og ikke godkjente i rødt
        figure(2),plot(tr(I,1),tr(I,2),'r.') %plot godkjente spor i grønt
    end
   
   %plot(trny(I,1))
end
figure(2),title 'Partikkelspor. Grønt spor: godkjent, rødt: forkastet'
hold off,axis equal
tr=trny;
numpart=length(trnum);

% kalkulere forflytningene
m=max(tr(:,3));
t=1:m;
sumdx=zeros(m,1); sumdx2=zeros(m,1);
sumdy=zeros(m,1); sumdy2=zeros(m,1); sumdxy2=zeros(m,1);
xstart = zeros(numpart,1);
ystart = zeros(numpart,1);
xstop = zeros(numpart,1);
ystop = zeros(numpart,1);
for n=1:numpart % loop gjennom alle partiklene
    I=find(tr(:,4)==trnum(n));
    % finne hvor langt partikkelen har flyttet seg fra startposisjonen ved 
    % hvert tidssteg
    dx=tr(I,1)-tr(I(1),1);
    dy=tr(I,2)-tr(I(1),2);
    sumdx=sumdx+dx.^2;
    sumdy=sumdy+dy.^2;
    sumdx2=sumdx2+dx.^4;
    sumdy2=sumdy2+dy.^4;
    sumdxy2=sumdxy2+(dx.^2+dy.^2).^2;
    xstart(n) = tr(I(1),1);
    ystart(n) = tr(I(1),2);
    xstop(n) = tr(I(end),1);
    ystop(n) = tr(I(end),2);
end
msdx=sumdx/numpart;
msdy=sumdy/numpart;
msdxy=(msdx+msdy);

msdx2 = sumdx2/numpart;
msdy2 = sumdy2/numpart;
msdxy2 = sumdxy2/numpart;

errdx = sqrt(msdx2-msdx.^2)/sqrt(numpart);
errdy = sqrt(msdy2-msdy.^2)/sqrt(numpart);
errdxy = sqrt(msdxy2-msdxy.^2)/sqrt(numpart);

% Plotte start og stopp punkter
figure(1), plot(xstart,ystart,'ro')
plot(xstop,ystop,'ko')
title('Partikkelspor med start (rød sirkel) og stopp (svart sirkel)')
%label('x [piksler]'), ylabel('y [piksler]')

% Plotte resultater 
figure(3), hold off
plotx=plot(msdx,'b');
xband = plotx.XData';
hold on;
% Plot error in y as function of x as a transparent blue band
f=fill([xband;flip(xband)],[msdx-errdx;flip(msdx+errdx)],'b','EdgeColor','None');
f.FaceAlpha=0.3;

plot(msdy,'g')
% Plot error in y as function of x as a transparent green band
f=fill([xband;flip(xband)],[msdy-errdy;flip(msdy+errdy)],'g','EdgeColor','None');
f.FaceAlpha=0.3;

plot(msdxy,'r')
% Plot error in y as function of x as a transparent red band
f=fill([xband;flip(xband)],[msdxy-errdxy;flip(msdxy+errdxy)],'r','EdgeColor','None');
f.FaceAlpha=0.3;

% finne stigningstall
% tvinge løsningen til å gå gjennom origo.
px = t(:)\msdx(:);
py = t(:)\msdy(:);
pxy = t(:)\msdxy(:);

xlabel('antall tidssteg')
ylabel('Midlere kvadratisk forflytning, piksler^2')
hold on, plot(t,px*t,'b'), plot(t,py*t,'g'), plot(t,pxy*t,'r'), hold off
title('Midlere forflytning av Brownske partikler')
axis tight

legend('<x^2>','\sigma_{<x^2>}','<y^2>','\sigma_{<y^2>}', '<x^2+y^2>','\sigma_{<x^2+y^2>}',...
    ['d<x^2>/dt=' num2str(px(1))],['d<y^2>/dt=' num2str(py(1))],['d<x^2+y^2>/dt=' num2str(pxy(1))],...
    'Location','NorthWest')

