clear all
close all

% Valg af partiklernes start position
k = menu('Skal Partiklerne i åen eller i indløbet?','I åen','I indløbet');
tic
z=0; %Viser hvor mange partikler der kommer ud af systemet

%Koordinatsystem defineres%
GridBredde = 150;
Gridhoejde = 55;

%Celler
CellLaengde  = 5;
AntalCellAdX = GridBredde / CellLaengde; %30
AntalCellAdY = Gridhoejde / CellLaengde; %11
StyreCell=zeros(AntalCellAdY,AntalCellAdX); %Kontrollere, hvor partiklerne er i det sidste tidsskridt

%Filter ind
FilterIndOffsetX  = 150;
FilterIndOffsetY  = 0;
FilterIndLaengdeX = 0;
FilterIndLaengdeY = 18;

%Filter ud
FilterUdOffsetX  = 0;
FilterUdOffsetY  = 37;
FilterUdLaengdeX = 0;
FilterUdLaengdeY = 18;

% Åen
Aax_start = 65;
Aax_slut = 85;
Aay_top = 55;
Aay_bund = 45;

% startbetingelser
beregning = 1125; %1125;              %Antal tidsskridt
%eff_por = 0.4;
eff_por=0.39*ones(11,30);
eff_por(1:3,13:18)=0.8;
eff_por(1:2,14:17)=0.9;
Ud = zeros(beregning,1);

if k == 2    
    partikler = 1040;              %Antal partikler
    Px_start = 150;
    Py_start = 0;
    Px = Px_start*ones(partikler,beregning);
    Py = Py_start*ones(partikler,beregning);
    for r=1:partikler
        Py(r,1)=FilterIndLaengdeY-(r*FilterIndLaengdeY/partikler);
    end
    % definition af hastighedsvektor
    Vy = 100*xlsread('vy.xls');
    Vx = 100*xlsread('vx.xls');
    Tid_tilsat=80;
else
    partikler = 4050;              %Antal partikler
    Px_start = (Aax_start+Aax_slut)/2;
    Py_start = (Aay_top+Aay_bund)/2;
    Px = Px_start*ones(partikler,beregning);
    Py = Py_start*ones(partikler,beregning);
    Tid_tilsat=150;
    position(:,1)=xlsread('Aa_part.xls','B4248:B1');
    position(:,2)=xlsread('Aa_part.xls','C4248:C1');
    for r=1:partikler
        l=int16(rand*4247)+1;
        Px(r,1)=position(l,1);
        Py(r,1)=position(l,2);
    end
%    for r=1:partikler
%        Px(r,1)=((Aax_slut-Aax_start)/partikler)*r+Aax_start;
%        Py(r,1)=((Aay_top-Aay_bund)/partikler)*r+Aay_bund;
%        ABE1(r,1)=Py(r,1);
%        ABE1(r,2)=Px(r,1);
%    end
    % definition af hastighedsvektor
    Vy = 100*xlsread('vy_vandloeb.xls');
    Vx = 100*xlsread('vx_vandloeb.xls');
end

%Definering af variabler%
dt = 60;
D=(2.13)*ones(11,30);
D(1:3,13:18)=5.76;
D(1:2,14:17)=0.0000001;

%selve beregning%
for r=1:partikler
    Pxi=Px(r,1);
    Pyi=Py(r,1);
    for     t = 1:beregning
        if (partikler/t>partikler/Tid_tilsat) & (r>(partikler/Tid_tilsat)*t)
            Px(r,t)=Px(r,1);
            Py(r,t)=Px(r,1);
        else
            if      (Pyi == 55)
                    Pyi=Pyi-0.001;
            end
            if      (Pxi<5) & (Pyi==0)
                    cx = 1;
                    cy = int32(Pyi/CellLaengde);        
                    Vx1 = Vx(11 - cy,cx);
                    Vx2 = Vx(11 - cy,cx + 1);    
                    Vy1 = Vy(11 - cy,cx);
                    Vy2 = Vy(12 - cy,cx);
                    cx=0;
            elseif  (Pyi == 0)     
                    cx = int32(Pxi/CellLaengde-0.5);
                    cy = int32(Pyi/CellLaengde);        
                    Vx1 = Vx(11 - cy,cx);
                    Vx2 = Vx(11 - cy,cx + 1);    
                    Vy1 = Vy(11 - cy,cx);
                    Vy2 = Vy(12 - cy,cx);
            elseif  Pxi<5
                    cx = 1;
                    cy = int32(Pyi/CellLaengde-0.5);        
                    Vx1 = Vx(11 - cy,cx);
                    Vx2 = Vx(11 - cy,cx + 1);    
                    Vy1 = Vy(11 - cy,cx);
                    Vy2 = Vy(12 - cy,cx);
                    cx=0;
            elseif  (Pxi>=5) & (Pyi>=0)  
                    cx = int32(Pxi/CellLaengde-0.5);
                    cy = int32(Pyi/CellLaengde-0.5);
                    Vx1 = Vx(11 - cy,cx);
                    Vx2 = Vx(11 - cy,cx + 1);    
                    Vy1 = Vy(11 - cy,cx);
                    Vy2 = Vy(12 - cy,cx);
            end 
    %interpolering af Vx og Vy
            Vxx=[Vx1 Vx2];
            Vyy=[Vy1 Vy2];
            boksx(1,1)=double((cx)*CellLaengde);
            boksx(1,2)=double((cx+1)*CellLaengde);    
            Vxi=interp1(boksx,Vxx,Pxi,'linear');        
            boksy(1,1)=double((cy+1)*CellLaengde);
            boksy(1,2)=double(cy*CellLaengde);
            Vyi=interp1(boksy,Vyy,Pyi,'linear');    
            if      (Pxi<5) & (Pyi==0)
                    dx = 1;
                    dy = 11-int32(Pyi/CellLaengde);
            elseif  (Pxi<5) & (Pyi==55)
                    dx = 1;
                    dy = 1;
            elseif  (Pyi == 55)
                    dx = int32(Pxi/CellLaengde-0.5);
                    dy = 1;
            elseif  (Pyi == 0)
                    dx = int32(Pxi/CellLaengde-0.5);
                    dy = 11-int32(Pyi/CellLaengde);
            elseif  Pxi<5
                    dx = 1;
                    dy = 11-int32(Pyi/CellLaengde-0.5);
            elseif  (Pxi>=5) & (Pyi>=0)  
                    dx = int32(Pxi/CellLaengde-0.5);
                    dy = 11-int32(Pyi/CellLaengde-0.5);
            end        
            Px(r,t) = Pxi+Vxi*dt/eff_por(dy,dx)+2*(rand-0.5)*sqrt(sqrt((2*D(dy,dx)*Vxi*dt/eff_por(dy,dx))^2));
            Py(r,t) = Pyi+Vyi*dt/eff_por(dy,dx)+2*(rand-0.5)*sqrt(sqrt((2*D(dy,dx)*Vyi*dt/eff_por(dy,dx))^2));        
            if Px(r,t)<FilterUdLaengdeX & Py(r,t)>FilterUdOffsetY
                z=z+1; %Når der kommer endnu en partikel ud ligges der en til =)
                %z,r,t; %Viser hvormange partikler der er kommet ud, hvilken partikel og tiden
                Ud(t,1)=Ud(t,1)+1;
            break %Afslutter løkken for den ene partikel der er kommet ud, resten af beregningerne bliver nul
            end        
            if Px(r,t)>GridBredde
                Px(r,t)=(2*GridBredde)-Px(r,t); %Spejler partiklen tilbage i X retningen ved indløbet
            elseif Px(r,t)<0
                Px(r,t)=-1*Px(r,t); %Spejler partiklen tilbage i X retningen ved udløbet
            end        
            if Py(r,t)>Gridhoejde
                Py(r,t)=2*Gridhoejde-Py(r,t); %Spejler partiklen tilbage i Y retningen i toppen
            elseif Py(r,t)<0
                Py(r,t)=-1*Py(r,t); %Spejler partiklen tilbage i Y retningen i bunden
            end
    % Klargøring til næste beregningstrin
            Pxi=Px(r,t);
            Pyi=Py(r,t);
            %hold on
            %axis([0 150 0 55]);
            %plot(Px(r,t),Py(r,t),'.')
            %f(t)=getframe;
            %hold off
    % Når sidste tidsskridt er beregnet findes antallet af partiklernes i de forskellige kasser        
            if t==beregning
                cx = int32(Px(r,t)/CellLaengde-0.5);
                cy = int32(Py(r,t)/CellLaengde-0.5);
                if cx<=0
                    cx=1;
                end
                cx=double(cx);
                if cy<0
                    cy=0;
                end
                cy=double(cy);
                StyreCell(11-cy,cx)=StyreCell(11-cy,cx)+1;
            end
        end    
    end
end
%movie(f,1,3)
if k == 2
    M=xlsread('indloeb.xls','A1:A225');
    N=xlsread('indloeb.xls','B1:B225');
else
    M=xlsread('vandloeb.xls','A1:A225');
    N=xlsread('vandloeb.xls','B1:B225');
end
Ud=Ud/partikler;
for r=1:(beregning/5)
    L=r*5;
    Udi(r,1)=L;
    Udi(r,2)=Ud(L-4,1)+Ud(L-3,1)+Ud(L-2,1)+Ud(L-1,1)+Ud(L,1);
end
figure
hold on
axis([0 150 0 55]);
plot(Px(:,beregning),Py(:,beregning),'.')
hold off
Tid_for_beregning=toc
figure
c=cell(size(StyreCell));
for r=1:AntalCellAdY
    for t=1:AntalCellAdX
        c(r,t)={StyreCell(r,t)};
    end
end
cellplot(c)
%figure
%hold on
%plot(Ud)
%plot(M,N)
%hold off
%XYZ=smooth(Udi(:,1),Udi(:,2),10,'moving');
%ASP=sort(N);
%Udz=sort(Udi(:,2));
%figure
%plot(Udz,ASP)
figure
hold on
plot(M,N)
plot(Udi(:,1),Udi(:,2))
hold off