Hoi4 Dispersed Vs Concentrated

clc
clear all
close all

% messures of growth of production efficency for production caps of
% 50%,60%,70% and 80%

pe50 = [
    1010
    1026
    1050
    1074
    1097
    1120
    1143
    1165
    1186
    1207
    1228
    1248
    1268
    1288
    1307
    1327
    1346
    1364
    1382
    1401
    1418
    1505
    1602
    1708
    1808
    1890
    1929
    1942
    1968
    1993
    2018
    2056
    2092
];

dpe50 = [
    247
    241
    238
    232
    226
    223
    217
    214
    211
    204
    201
    198
    195
    192
    189
    186
    183
    183
    180
    177
    174
    165
    156
    143
    137
    131
    128
    128
    125
    125
    122
    119
    119
];

pe60 = [
    1010
    1072
    1138
    1170
    1201
    1230
    1260
    1288
    1317
    1344
    1371
    1397
    1448
    1497
    1545
    1569
    1614
    1659
    1702
    1744
    1785
    1825
    1865
    1903
    1941
    1978
    2014
    2025
    2085
];

dpe60 = [
    357
    336
    314
    305
    299
    290
    284
    278
    272
    266
    262
    256
    247
    238
    232
    229
    223
    217
    211
    204
    201
    195
    192
    189
    183
    180
    177
    174
    171
];

pe70 = [
  1010
  1050
  1096
  1141
  1184
  1226
  1305
  1379
  1449
  1548
  1580
  1611
  1642
  1672
  1701
  1786
  1814
  1867
  1894
  1920
  1970
  1995
  2020
  2044
  2092
];

dpe70 = [
    488
    467
    446
    427
    412
    400
    375
    354
    336
    314
    308
    302
    296
    293
    287
    272
    269
    262
    256
    253
    247
    244
    241
    235
    232
];


pe80 = [
    1001
    1065
    1125
    1182
    1236
    1288
    1338
    1432
    1520
    1563
    1603
    1643
    1682
    1721
    1794
    1830
    1899
    1933
    1999
    2031
];

dpe80 = [
    638
    601
    568
    540
    516
    494
    476
    446
    421
    409
    397
    388
    378
    369
    354
    348
    336
    330
    317
    314
];

% i copied just the digits from in game percentages so fix it here
pe50=pe50./10000;
dpe50=dpe50./100000;

pe60=pe60./10000;
dpe60=dpe60./100000;

pe70=pe70./10000;
dpe70=dpe70./100000;

pe80=pe80./10000;
dpe80=dpe80./100000;

% pe cap for each messured growth, note that those are a bit off from
% the the expected round values, but the game reports just like it
peCap = [0.5005 0.6006 0.7007 0.8008];

%
% do wierd matlab stuff ...
%

% find out that the rate of change is the inverse of the function value
% times production efficency cap squared
dpe = @(pe,cap) 1/1000 * cap.^2./ pe ;

%
% check rate of change against game data
%
figure
plot([pe50 ;pe60 ;pe70; pe80],[dpe50 ;dpe60 ;dpe70 ;dpe80] ,'ko')
hold on
plot(pe50,dpe(pe50,peCap(1)),'r')
plot(pe60,dpe(pe60,peCap(2)),'g')
plot(pe70,dpe(pe70,peCap(3)),'b')
plot(pe80,dpe(pe80,peCap(4)),'c')
legend('In Game Data','cap 50%','cap 60%','cap 70%','cap 80%')
xlabel('P_e')
ylabel('dP_e')
hold off

%
% reproduce paradox graph result
%

% integrate to get production efficency at t days given base efficency p0
pe=@(t,cap,p0) min(cap , sqrt(2) * sqrt(cap^2*t/1000+1/2*p0^2));

days=0:500;
cap=0.75;
base=0.10;
% guessed old growth function
oldEfficencyGrowthFunction=@(t) min(0.002407*t+base,cap);
figure
plot(days,oldEfficencyGrowthFunction(days),'-r','LineWidth',2)
hold on
plot(days,pe(days,cap,base),'b','LineWidth',2);
legend('old linear','new log like approximation')
hold off

%
% figure out if concentrated is better than dispersed
%

% free trade bonus and tech level
tech=3;
trade=0.15;

baseConcentrated = 1+trade+tech*0.15;
baseDispersed= 1+trade+tech*0.10;
cap = 0.5005 + tech * 0.1001;

p0Concentrated = 0.10;
p0Dispersed = 0.10+tech*0.05;

days = 0:1000;
figure
plot(days, pe(days,cap,p0Concentrated) , ...
     days, pe(days,cap,p0Dispersed))

xlabel('days')
ylabel('P_e')
legend('concentrated','disparsed')

figure
plot(days,Fout(days,cap,baseConcentrated,p0Concentrated), ...
     days,Fout(days,cap,baseDispersed,p0Dispersed));

xlabel('days')
ylabel('F_{out}')
legend('concentrated','disparsed')

figure
days=10:1000;
largeTimeAsymptote=(baseDispersed-baseConcentrated)/baseConcentrated;
plot(days, ( F_out(days,cap,baseDispersed,p0Dispersed)...
              -F_out(days,cap,baseConcentrated,p0Concentrated))....
              ./F_out(days,cap,baseConcentrated,p0Concentrated),'b',...
              [0 1000], [0 0] ,'k',...
              [0 1000], [largeTimeAsymptote largeTimeAsymptote] ,'--k');
legend('relative output gained from disparsed')
xlabel('days')