#!/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:
%   Katzen und Hunde je eins, bis alle satt sind.
%   Der (kleine) Rest des Tages wird weggeworfen.
% - Natürlich gestorben wird vor dem Werfen:
%   Diese Kadaver bilden Futter!

% Sättigungsfaktor pro geschlachtetem Individuum
global FK = WK * GK / EH;
global FH = WH * GH / EK;

% nachhaltiges Schlachten
function [nh, nk] = schlachtung(nh, nk, th, tk)
   global FK FH;
   % hungrige Tiere, die auf Futter warten
   hk = nk; hh = nh;
   % Nachhaltigkeit: Aas einrechnen
   nh = nh + th; nk = nk + tk;
   % Solange gehungert wird und es Futter gibt
   while((hh > 0 || hk > 0) && nh >= th && nk >= tk)
      % eine Katze schlachten
      if (hh > 0)
         nk = nk - 1;
         hk = hk - 1; % Nachhaltigkeit: weniger Mäuler
         hh = hh - FK; % Nachhaltigkeit: Rest behalten
      endif
      % einen Hund schlachten
      if (hk > 0)
         nh = nh - 1;
         hh = hh - 1; % Nachhaltigkeit: weniger Mäuler
         hk = hk - FH; % Nachhaltigkeit: Rest behalten
      endif
   end
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
   TK = ceil(AK / LK / 365); AK = AK - TK;
   TH = ceil(AH / LH / 365); AH = AH - TH;
   % Reproduktion, pessimistische Schätzung
   AK = AK + floor(AK / LK / 365 * RH);
   AH = AH + floor(AH / LH / 365 * RH);
   % Schlachtung, auch Aas
   [AH, AK] = schlachtung(AH, AK, TH, TK);
   % 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];
   printf("%3d, %e\n", d, E(end)); % Fortschrittsanzeige
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: sustainable slaughter");
xlabel("time [days]");
ylabel("population [mlns]");
print("d+c-eigenwert3a.png", "-dpng", "-S960,640");
print("d+c-eigenwert3a.svg", "-dsvg", "-S960,640");