%Grundvandsmodel for "Sandkassen" (alle enheder i meter, sekunder)

clear
clc
close all


%Variable inputparamatre---------------------------------------------------
r=11;                   %antal rækker 11,22,33,44,55
s=30;                   %antal søjler, 30, 60, ...
ha=0.456;                %Trykniveau i venstre kolonne, randbetingelse
hb=0.463;                %Trykniveau i højre kolonne, randbetingelse
hv=0.527;                %Trykniveau ved vandløb, randbetingelse
K1=0.000712%+1.424e-4;             %Konduktivitet for grovsand m/s
K2=(5.91e-5)+(5.910e-6);              %Konduktivitet for finsand m/s %7.23e-5 (0.01*d10^2)


vandloebsrand='j';      %'j' betyder, at åen er 'aktiv'. 
                        %'n' betyder, at modellen regner vandstanden i åen
filterkorrektion='j';   % skal der korrigeres for enkelttab ved filtre
stromlinier='n';        % ønsker du at plotte strømlinier? j,n
hastighedsvektorer='n'; % ønsker du at plotte hastighedsvektorer? j,n

 if filterkorrektion=='j'
 qa_exp=-0.0000045637      % forventet til beregning af c1 m3/s (ved ha)
                           % (hvis værdien er 0 medregnes c1 ikke)
 qb_exp=0.0000035409       % forventet Q til beregning af c2 m3/s (ved hb)
 qv_exp= 0.0000081046      % forventet Q til beregning af c3 m3/s (ved hv)
 c1=0.527;              % 1/s
 c2=0.527;              % 1/s
 c3=1.74;               % 1/s
 a1_exp=0.05*0.40;      % areal af filter i venstre side som gennemstrømmes
 end



%--Herfra nix pille--------------------------------------------------------
if filterkorrektion=='j'    %Her korrigeres for filtertab!
ha=ha-qa_exp/(a1_exp*c1)
disp('centimeter')
hb=hb+qb_exp/(0.182*0.4*c2)
disp('centimeter')
hv=hv-qv_exp/((0.20+0.06)*0.4*c3)
disp('centimeter')
filttertab_maks=1e3*max(abs([qa_exp/(a1_exp*c1);qb_exp/(0.182*0.4*c2);qv_exp/((0.20+0.06)*0.4*c3)]))
disp('milimeter')
end
B=1.5;                  %x
H=0.55;                 %y
dybde=0.4;              %z

B_vandloeb=0.3;         %Bredde af vandløb inkl. finsand
H_vandloeb=0.15;        %Højde af vandløb inkl. finsand
filter=0.182;           %højde af filter ved ha, hb
t_finsand=0.05;         %Tykkelse af finsand omkring vandløbet



rmax=r+2;
smax=s+2;
celler=r*s;             %Antal celler i modellen
         
%-Grid inddeling-----------------------------------------------------------
dy=H/r;          %m
dx=B/(s);
dz=dybde;

%koeficienter--------------------------------------------------------------
C_i_plus1_j    = dx*dz/dy;
C_i_j_plus1    = dy*dz/dx;
C_i_j_minus1   = dy*dz/dx;
C_i_minus1_j   = dx*dz/dy;

antal=int32(filter/dy);  %antal celler der udgør randene ved ha og hb

model= zeros(rmax,smax);            %Hjælpematrice som har 1-taller der hvor vi har en celle
model(2:rmax-1,2:smax-1)=1;                          
model_ij=model;                     %Denne er magen til model, bortset fra, at den aktiverer randcellerne så de indgår i Cij. Cij er IKKE ens for alle celler!!
model_ij(2:antal+1,1)=1;
model_ij(r+1-antal+1:r+1,s+2)=1;
if vandloebsrand=='j'
    model_ij(1,smax/2-1:smax/2+2)=1;
end

%Her opstilles K matricen og Keff regnes for x og y retningen--------------
model_K=model*K1;
model_K(2:int32(H_vandloeb/dy)+1,smax/2-int32(B_vandloeb/dx)/2+1:smax/2+int32(B_vandloeb/dx)/2)=K2;
model_K(2:int32((H_vandloeb-t_finsand)/dy)+1,smax/2-int32((B_vandloeb-2*t_finsand)/dx)/2+1:smax/2+int32((B_vandloeb-2*t_finsand)/dx)/2)=1;


model_Keff_x=zeros(r+2,s+1);
for i=1:(r+2)*(s+1)
    if (model_K(i)>0 & model_K(i+r+2)>0)
model_Keff_x(i)=dx/((dx/2)/model_K(i)+(dx/2)/model_K(i+r+2));
    end
end
model_Keff_x(2:antal+1,1)=2*K1;
model_Keff_x(r+1-antal+1:r+1,s+1)=2*K1;
model_Keff_x;

model_Keff_y=zeros(r+1,s+2);
for i=1:(r+1)
    for j=1:(s+2)
        if (model_K(i,j)>0 & model_K(i+1,j)>0)
model_Keff_y(i,j)=dy/((dy/2)/model_K(i,j)+(dy/2)/model_K(i+1,j));
        end
    end
end
if vandloebsrand=='j'
model_Keff_y(1,smax/2-1:smax/2+2)=2;
end
%--------------------------------------------------------------------------


j_plus1= model(2:rmax-1,3:smax)*C_i_j_plus1;
for i=1:r
    for j=1:s
    j_plus1(i,j)=j_plus1(i,j)*model_Keff_x(i+1,j+1); %Her ganges de rigtige Keff på 
    end
end

j_minus1= model(2:rmax-1,1:smax-2)*C_i_j_minus1;
for i=1:r
    for j=1:s
    j_minus1(i,j)=j_minus1(i,j)*model_Keff_x(i+1,j);
    end
end

i_plus1= model(3:rmax,2:smax-1)*C_i_plus1_j;
for i=1:r
    for j=1:s
    i_plus1(i,j)=i_plus1(i,j)*model_Keff_y(i+1,j+1);
    end
end

i_minus1= model(1:rmax-2,2:smax-1)*C_i_minus1_j;
for i=1:r
    for j=1:s
    i_minus1(i,j)=i_minus1(i,j)*model_Keff_y(i,j+1);
    end
end

ij=model(2:rmax-1,2:smax-1);
for i=1:r
    for j=1:s
    ij(i,j)=-model_ij(i,j+1)*C_i_minus1_j*model_Keff_y(i,j+1)...
        -model_ij(i+2,j+1)*C_i_plus1_j*model_Keff_y(i+1,j+1)...
        -model_ij(i+1,j)* C_i_j_minus1*model_Keff_x(i+1,j)...
        -model_ij(i+1,j+2)*C_i_j_plus1*model_Keff_x(i+1,j+1);
    end
end

%-------Dannelse af koefficientmatrice------------------------------------
 
 A=zeros(celler);

 for m=1:celler
    for n=1:celler
        %m=række
        %n=søjle
        if n==m
        A(m,n)=ij(m);
        end

        if (n==m+1 |n==m-(celler-1))
        A(m,n)= i_plus1(m);
        end

        if (n==m+r | n==m-(celler-r))
        A(m,n)= j_plus1(m);
        end

        if (n==m+(celler-r) | n==m-r)
        A(m,n)= j_minus1(m);    
        end

        if (n==m+(celler-1) | n==m-1) 
        A(m,n)= i_minus1(m);% duer ikke ved 1 række
        end
    end
 end

%------D matrice-----------------------------------------------------------
D=zeros(celler,1);

D(1:antal)=-K1*dy*dz/(0.5*dx)*ha; %Skal altid være K1???
D(celler-antal+1:celler)=-K1*dy*dz/(0.5*dx)*hb; %Her er trykniveauerne ganget på men ikke i modellen ellers... 
if vandloebsrand=='j'
    D((s/2-2)*r+1)=-dx*dz/(0.5*dy)*hv;
    D((s/2-1)*r+1)=-dx*dz/(0.5*dy)*hv;
    D((s/2)*r+1)=-dx*dz/(0.5*dy)*hv;
    D((s/2+1)*r+1)=-dx*dz/(0.5*dy)*hv;
    end
h=A\D;


%Hastigheder og massebalance-----------------------------------------------
h_m=zeros(r,s+2); % sætter h op i en matrice
for i=1:r*s
    h_m(i+r)=h(i);
end
h_m(1:antal,1)=ha;
h_m(r-antal+1:r,s+2)=hb;
%h_m

Vx=zeros(r,s+1);
for i=1:r
    for j=1:(s+1)
    Vx(i,j)=-model_Keff_x(i+1,j)*(h_m(i,j+1)-h_m(i,j))/dx;
    end
end

Vy=zeros(r-1,s);
for i=1:(r-1)
    for j=1:s
    Vy(i,j)=-model_Keff_y(i+1,j+1)*(h_m(i+1,j+1)-h_m(i,j+1))/dy;
    end
end

Qb=sum(Vx(r-antal+1:r,s+1)*dy*dz);%m3/s

Qa=sum(Vx(1:antal,1)*dy*dz); %m3/s
if vandloebsrand=='j'
    Qv=-Qa+Qb
else
masseafvigelse=Qa-Qb    
end

    
h1=(h_m(antal,2)+h_m(antal-1,2))/2;
h2=(h_m(r,2)+h_m(r-1,2))/2;
h3=(h_m(int32(H_vandloeb/dy)+1,smax/2)+...
    h_m(int32(H_vandloeb/dy)+1,smax/2+1))/2;
h4=(h_m(r,smax/2)+...
    h_m(r,smax/2+1))/2;
h5=(h_m(antal,s+1)+h_m(antal-1,s+1))/2;
h6=(h_m(r,s+1)+h_m(r-1,s+1))/2;
hvand=(h_m(1,smax/2));
h_beregnet=[K1*1e3;K2*1e4;Qa*1e5;Qb*1e5;ha;hb;h1;h2;h3;h4;h5;h6;hvand]
%----------------Hastighedsvektorer plottes--------------------------------
for i=1:r
    for j=1:s
Vx_plot(i,j)=(Vx(i,j)+Vx(i,j+1))/2;
    end
end
Vx_plot(antal+1:r,1)=0;
Vx_plot(2:r-antal,s)=0;
Vy_plot=zeros(r,s);
for i=1:r-2
    for j=1:s
        Vy_plot(i+1,j)=(Vy(i,j)+Vy(i+1,j))/2;
    end
end
Vy_plot(1,:)=0;
Vy_plot(r,:)=0;
x_koor=zeros(r,s);
for i=1:r
    for j=1:s
        x_koor(i,j)=(j-1)*dx+dx/2;
    end
end
y_koor=zeros(r,s);
for i=1:r
    for j=1:s
        y_koor(i,j)=(0.55-dy/2)-(i-1)*dy;
    end
end
    if hastighedsvektorer=='j'
        scale=2;
        quiver(x_koor,y_koor,Vx_plot,-Vy_plot,scale,'r')
    end
        %------strømlinieplot------------------------------------------------------
    if stromlinier =='j'
        if vandloebsrand=='j'
        streamline(x_koor,y_koor,Vx_plot,-Vy_plot,0.65:0.01:0.85,0.45*ones(21,1))
        end
        if vandloebsrand=='n'
        streamline(x_koor,y_koor,Vx_plot,-Vy_plot,1.475*ones(21,1),0:0.01:0.2)
        end
    end

Q_beregnet=[Qa;Qb;Qv]