%% Bakjes model Waternet
% TO 090910

XLSF='bakje_MP1_TO.xls';

if ismac, Meteo=xlsread(XLSF,'NE'    ,'','basic'); end
if ispc,  Meteo=xlsread(XLSF,'NE');                end
% on mac save xls file in excel5.0/95 format !!

Meteo(isnan(Meteo(:,1)),:)=[];
Meteo(isnan(Meteo(:,2)),:)=[];
Meteo(isnan(Meteo(:,3)),:)=[];
Meteo(:,1)=Meteo(:,1)+datenum(1899,12,31)-1; % Excel to Matlab dates

t =Meteo(:,1); NT=length(t);
P =Meteo(:,2);
EM=Meteo(:,3);

%% Get fixed data for all bakjes

if ismac, [Values,Labels]  =xlsread(XLSF,'Invoer','','basic'); end
if ispc,  [Values,Labels]  =xlsread(XLSF,'Invoer');            end

for i=1:length(Labels(:,1))
    if isempty(Labels{i,1}),
        Labels(i:end,:)=[];
        break;
    end
end
for i=1:length(Labels(1,:))
    if isempty(Labels{1,i}) || ~isempty(findstr('dummy',Labels{1,i}));
        Labels(:,i:end)=[];
        break;
    end
end
bakjeTypeName=lower(Labels(1,2:end));
Field=lower(Labels(2:end,1));

NBakjes=length(bakjeTypeName);


if all(isnan(Values(1,:))), Values(1,:)=[]; end  % verschil Mac en PC
Values=Values(1:length(Field),1:NBakjes);

bakje(NBakjes).name=bakjeTypeName{end};
for ib=1:NBakjes
    bakje(ib).name        =bakjeTypeName{ib};
    bakje(ib).opp         =Values(strmatch('perceelsopp',   Field),ib);
    bakje(ib).mv          =Values(strmatch('maaiv',         Field),ib);
    bakje(ib).drn         =Values(strmatch('drainage',      Field),ib);
    bakje(ib).sy          =Values(strmatch('yield',         Field),ib);
    bakje(ib).kwel        =Values(strmatch('kwel (mm/d',    Field),ib);
    bakje(ib).cropfac     =Values(strmatch('gewasverdampingsfactor',Field),ib);
    bakje(ib).grwflushfrac=Values(strmatch('fractie van berging die uitspoelt',Field),ib);
    bakje(ib).initialSpace=Values(strmatch('initiele ruimte in bak (0-1)',Field),ib);
    bakje(ib).cap =1000*...
        (bakje(ib).mv-bakje(ib).drn)*bakje(ib).sy;  % in mm
    
    % gerioleerd en afgekoppeld hebben nauwelijks bering, 1 mm aangehouden
    if ~isempty(findstr('gerioleerd', bakje(ib).name)), bakje(ib).cap=1.0; end
    if ~isempty(findstr('afgekoppeld',bakje(ib).name)), bakje(ib).cap=1.0; end
    
    bakje(ib).seep  =NaN(NT,1);
    bakje(ib).stored=NaN(NT,1);
    bakje(ib).space =NaN(NT,1);
    bakje(ib).runoff=NaN(NT,1);
    bakje(ib).shortage =NaN(NT,1);
    bakje(ib).Err   =NaN(NT,1);
   

%%  Compute, everything in mm

    for it=1:NT
        if it==1,  % initiele waarden
            spacePrev =bakje(ib).cap*bakje(ib).initialSpace;
            storedPrev=bakje(ib).cap-spacePrev;
        else       % waarden onderweg voor it>1
            spacePrev =bakje(ib).space(it-1);
            storedPrev=bakje(ib).stored(it-1);
        end

        P =Meteo(it,2);  % Neerslag
        EM=Meteo(it,3);  % Makkink verdamping

        bakje(ib).seep(it)=storedPrev*bakje(ib).grwflushfrac;  % downward seepage loss

        %% Bijzonder bakjes (onderliggende) apart behandelen
        if isempty(strfind('onder drain',bakje(ib).name))
            In= P-EM*bakje(ib).cropfac+bakje(ib).kwel-bakje(ib).seep(it);
        else
            In= bakje(ib-1).seep(it)+...
                bakje(ib-1).runoff(it)-bakje(ib-1).shortage(it)+...
                bakje(ib).kwel-bakje(ib).seep(it);
        end
        
        
        %% Kern van het rekenhart
        if In>0  % toevoer is positief
            if In>spacePrev,                % berging stroomt over
                bakje(ib).space(it)=0;                 % ruimte voor berging is op
                bakje(ib).stored(it)=bakje(ib).cap;    % berging zit vol
                bakje(ib).runoff(it)=In-spacePrev;     % runoff treedt op
                bakje(ib).shortage(it)=0;              % er is geen tekort
            else                            % berging stroomt niet over
                bakje(ib).space(it) =spacePrev -In;    % ruimte voor bering neemt af
                bakje(ib).stored(it)=storedPrev+In;    % bering neemt toe
                bakje(ib).runoff(it)=0;                % er is geen runoff
                bakje(ib).shortage(it) =0;             % er is geen tekort
            end
        else  % In<=0  toevoer is negatief
            if storedPrev+In<0              % berging raakt geheel op
                bakje(ib).stored(it)=0;                % berging is nu leeg
                bakje(ib).space(it)=bakje(ib).cap;     % ruimte voor berging is maximaal
                bakje(ib).runoff(it)=0;                % er is een runoff
                bakje(ib).shortage(it)=-In-storedPrev;  % maar wel een tekort
            else                            % berging raakt niet leeg
                bakje(ib).stored(it)=storedPrev+In;    % berging neemt af
                bakje(ib).space(it) =spacePrev -In;    % ruimte voor berging neemt toe
                bakje(ib).runoff(it)=0;                % er is geen surface runoff
                bakje(ib).shortage(it)=0;              % er is geen tekort
            end
        end
        %% Sluitterm (geldt voor boven en onderbakjes evenzeer
        bakje(ib).Err(it)=In-bakje(ib).runoff(it)+bakje(ib).shortage(it)-(bakje(ib).stored(it)-storedPrev); % balance should be zero

        %  Voor gecombineerde bakjes kunnen de termen worden opgesteld
        if ~isempty(strfind('onder drain',bakje(ib).name))
            bakje(ib).Err(it)=P-E*bakje(ib-1).cropfactor...
                              +bakje(ib).kwel-bakje(ib).seep(i)...   % dit is de In voor het onder drain bakje
                              -bakje(ib-1).runoff(it)-bakje(ib).runoff(it)...       % runoff bovenbakje
                             -(bakje(ib-1).stored(it)-bakje(ib-1).stored(it-1))...
                             -(bakje(ib  ).stored(it)-bakje(ib  ).stored(it-1)); % balance should be zero
        end

    end

end


%% Visualisatie

close all

%% Peil boven en onder de drain

figure
plot(t,bakje(1).mv-bakje(1).space/1000./bakje(1).sy,'g'); hold on;
plot(t,bakje(2).mv-bakje(2).space/1000./bakje(2).sy,'b'); hold on;
legend('peil bakje 1 (boven drain)','peil bakje 2 (onder drain)');
title('peil bakje(1) en bakje(2)'); ylabel('m NAP'); datetick

%% Inzicht in de verschillende bakjes

for ib=1:length(bakje)
    figure
    subplot(3,1,1);
    plot(t,bakje(ib).stored); hold on
    title(['bakje: ',bakje(ib).name,'; Item: stored']); ylabel('mm'); datetick;
    subplot(3,1,2);
    plot(t,bakje(ib).shortage,'r'); hold on
    title(['bakje: ',bakje(ib).name,'; Item: shortage']); ylabel('mm/d'); datetick;
    subplot(3,1,3);
    plot(t,bakje(ib).runoff,'b'); hold on
    title(['bakje: ',bakje(ib).name,'; Item: runoff']); ylabel('mm/d'); datetick;
end