#!/usr/bin/octave
global NH = 5300000 NK = 7900000; 
global EH = 1000    EK = 400;     
global GH = 8.6     GK = 4.9;     
global WH = 1100    WK = 1100;    
global LH = 11.5    LK = 12.5;    
global RH = 15      RK = 20;      
function FH = hundefutter(nh)
   global EH WK GK;
   fh = nh * EH;
   wk = WK * GK;
   FH = ceil(fh / wk);
end
function FK = katzenfutter(nk)
   global EK WH GH;
   fk = nk * EK;
   wh = WH * GH;
   FK = ceil(fk / wh);
end
H = [NH]; K = [NK]; 
for d = 1:1:365
   
   AK = K(end);
   AH = H(end);
   
   AK = AK - ceil(AK / LK / 365);
   AH = AH - ceil(AH / LH / 365);
   
   AK = AK + floor(AK / LK / 365 * RH);
   AH = AH + floor(AH / LH / 365 * RH);
   
   AK = AK - hundefutter(AH);
   AH = AH - katzenfutter(AK);
   
   if (AK <= 0) K = [K; NaN]; break; endif
   if (AH <= 0) H = [H; NaN]; break; endif
   
   K = [K; AK];
   H = [H; AH];
end
clf; hold on;
set(gca, ["fontsize"], 12);
ak = plot(K, "r;cats;o-");
ah = plot(H, "g;dogs;o-");
ylim([ylim()(1), ylim()(2) + 500000]);
set(gca, ["y" "ticklabel"], get(gca, ["y" "tick"]) / 1000000);
title("Germany: cats first");
xlabel("time [days]");
ylabel("population [mlns]");
print("d+c-germany1.png", "-dpng", "-S960,640");
print("d+c-germany1.svg", "-dsvg", "-S960,640");