%Modellering af sporstof
close all
clear
clc
%Konstanter
X=1000;             %Længde af åstrækning (m)
dX=5;               %Dellængder (m)
Q=0.34;             %Vandføring (m3)
A=1.87;              %Tværsnitsareal (m2)
v=Q/A               %Hastighed (m/s)
D=0.44;             %Dispersionshastighed (m2/s)
dt=10;              %Tidsskridt (s)
bereg_tid=2;      %timer
aat=bereg_tid*3600/dt;      %Antal tidsskridt
k_masse=0;
D_masse=0;

%--------------------Konstanter til magasinering

u=0.045;             %udveksling i m3
A_m=0.29      %Areal af magasin m2                  

%--------------------Startbetingelser

c(1,1:X/dX)=0;
c(1,1)=24.0*10^-5;    %kg/m^3 
Crc=v*dt/dX         %[-]
Crd=D*dt/dX^2       %[-]
Pe=Crc/Crd          %=v*dx/D
magasin=zeros(X/dX-1,1);

for t=1:aat;    
    c(t+1,1)=c(t,1)+dt*(-v*(c(t,1)+c(t,2))/(2*dX)+D*(c(t,2)-c(t,1))/dX^2);
    for x=2:X/dX-1;
    k=-v*(c(t,x+1)-c(t,x-1))/2/dX;
    m=D*(c(t,x-1)-2*c(t,x)+c(t,x+1))/dX/dX;
    c(t+1,x)=c(t,x)+dt*(k+m);
    magasin(x-1)=magasin(x-1)+c(t+1,x)*u/A_m/dX-magasin(x-1)*u/A_m/dX;
    c(t+1,x)=c(t+1,x)+magasin(x-1)*u/A/dX-u*c(t+1,x)/A/dX;
    end
% -------------Videoklip-laves----------------------------
% plot(c(t,:));
% set(gca,'xlim',[0 200],'ylim',[0 3e-5])
% set(gcf,'Color',[0.886 0.788 0.447 ])
% text(100,2.5e-5,['tid = (sek)',int2str(t*10)],'FontSize',18)
% text(175,0.05e-5,'1000 m','FontSize',12)
% text(10,2.8e-5,'3e-5 kg/m3','FontSize',12)
% F(t)=getframe(gca);
end
% mov = avifile('stoftransport_FD_model.avi')
% mov = addframe(mov,F);
% mov = close(mov);
% set(gcf,'Color',[0.8,0.88,0.6])
 x_akse=0:dt/60:bereg_tid*60;
 plot(x_akse,c(:,470/dX),'linewidth',2)
title('Rodaminkoncentration efter 470 m','fontsize',14)
 xlabel('minutter')
 ylabel('kg/m^3')
%målte data
%figure
hold on
graf=xlsread('rodamin.xls');
plot(graf(:,1)/60,graf(:,2)*10e-7,'r','linewidth',2)
xlabel('minutter')
ylabel('kg/m^3')

for i=1:X/dX
massebalance(i)=sum(c(:,i));
end
figure
set(gcf,'Color',[0.8,0.88,0.6])
plot(dX:dX:X,massebalance,'r','linewidth',2)
title('Massebevarelse','fontsize',14)
xlabel('Meter')
ylabel('Areal under kurven (kg*minutter/M^3)')
areal1=sum(c(1:4020/dt,470/dX))*dt
areal2=sum(graf(:,2))*3*dt
forskel=areal1-areal2