#!/usr/bin/octave

% Grundlagen:
%      Hunde        Katzen
global NH = 1000000 NK = 2091694; % Anzahl zu Beginn - Eigenwert
global EH = 1000    EK = 400;     % Ernährungsbedarf - kcal/Tag
global GH = 8.6     GK = 4.9;     % Gewicht - kg
global WH = 1100    WK = 1100;    % Nährwert - kcal/kg
global LH = 11.5    LK = 12.5;    % Lebenserwartung - Jahre
global RH = 15      RK = 20;      % Reproduktionsrate - Nachkommen/Leben

% weitere Annahmen:
% - Hunger wird immer gestillt.
% - Geschlachtet wird einmal am Tag:
%   zuerst die Hunde, dann die Katzen.
%   Der Rest des Tages wird weggeworfen.
% - Natürlich gestorben wird vor dem Werfen:
%   Diese Kadaver bilden kein Futter.

% Bedarf an Hundefutter in Katzen
function FH = hundefutter(nh)
   global EH WK GK;
   fh = nh * EH;
   wk = WK * GK;
   FH = ceil(fh / wk);
end

% Bedarf an Katzenfutter in Hunden
function FK = katzenfutter(nk)
   global EK WH GH;
   fk = nk * EK;
   wh = WH * GH;
   FK = ceil(fk / wh);
end

% naive Berechnung, numerisch
H = [NH]; K = [NK]; % täglich aktuelle Anzahl Hunde, Katzen
E = [NH * GH * WH + NK * GK * WK]; % Gesamte Energie im System
for d = 1:1:365
   % ein neuer Tag
   AK = K(end);
   AH = H(end);
   % natürlicher Tod, pessimistische Schätzung
   AK = AK - ceil(AK / LK / 365);
   AH = AH - ceil(AH / LH / 365);
   % Reproduktion, pessimistische Schätzung
   AK = AK + floor(AK / LK / 365 * RH);
   AH = AH + floor(AH / LH / 365 * RH);
   % Schlachtung
   AH = AH - katzenfutter(AK);
   AK = AK - hundefutter(AH);
   % Ende bevor alle Hunde oder Katzen tot sind
   if (AK <= 0) K = [K; NaN]; break; endif
   if (AH <= 0) H = [H; NaN]; break; endif
   % es ging sich aus
   E = [E; AH * GH * WH + AK * GK * WK];
   K = [K; AK];
   H = [H; AH];
end

% Zeichnen
clf; hold on;
set(gca, ["fontsize"], 12);
ak = plot(K, "r;cats;o-");
ah = plot(H, "g;dogs;o-");
gj = plot(E .* 0.0000041868, "b;energy;o-")
set(gca, ["y" "ticklabel"], get(gca, ["y" "tick"]) / 1000000);
title("Eigenwert: dogs first");
xlabel("time [days]");
ylabel("population [mlns]");
print("d+c-eigenwert2.png", "-dpng", "-S960,640");
print("d+c-eigenwert2.svg", "-dsvg", "-S960,640");