format compact 

function err = crit_electrons(ref_electrons, measured_electrons, reference_isos)
    isos_estimated  = 100 ./ (measured_electrons / ref_electrons);
    err = norm(log2(isos_estimated) - log2(reference_isos));
end

function ref_electrons = fit_electrons(measured_electrons, isos_dxo)
    ref_electrons = fminsearch(@(e) crit_electrons(e, measured_electrons, isos_dxo), 50000);
    ref_electrons = round(ref_electrons);
end

isos                = [ 100    200    400    800    1600   3200   6400 ];
electrons_photo_5d3 = [ 67684 34146  16835  8735   4288   2197   1085  ];
electrons_photo_6d  = [ 79546 42220  20643  10482  5215   2568   1271  ]
isos_dxo_5d3        = [ 80    160    323    641    1280   2518   5179  ];
isos_dxo_6d         = [ 80    153    311    616    1210   2400   4991  ];

iso_difference      = log2(mean(isos_dxo_5d3 ./ isos_dxo_6d));
ref_electrons_5d3   = fit_electrons(electrons_photo_5d3, isos * 2.^(iso_difference/2));
ref_electrons_6d    = fit_electrons(electrons_photo_6d, isos * 2.^(-iso_difference/2));

# quantum efficiency 5D3 vs 6D
qe_ratio = log2(ref_electrons_6d / ref_electrons_5d3 * (5496*3670) / (5796*3870))


# 5D3
isos =      [ 100    200    400    800    1600   3200   6400  12800 ];
dr_720 =    [ 11.02  10.96  10.93  10.85  10.67  10.36  9.89  8.89  ];
dr_1080 =   [ 11.10  11.05  11.02  10.93  10.73  10.39  9.88  8.88  ];
dr_crop =   [ 11.22  11.13  11.00  10.82  10.32   9.68  8.92  7.70  ];
dr_photo =  [ 10.92  10.83  10.69  10.46  10.02   9.40  8.49        ];
dr_max = max(dr_720, dr_1080);
dr_max = max(dr_max, dr_crop);

isos_720    = round(100 ./ ([1053242 643219 322555 163387 83083 41740 20180 10231] / ref_electrons_5d3 / 15))
isos_1080   = round(100 ./ ([770966 394170 199880 101404 50773 25974 12471 6278] / ref_electrons_5d3 / 9))
isos_crop   = round(100 ./ ([88129 44552 22018 11124 5648 2810 1376 717] / ref_electrons_5d3))
isos_photo  = round(100 ./ ([67684 34146 16835 8735 4288 2197 1085] / ref_electrons_5d3))

semilogx(isos_720, dr_720, '-or;720p;', 'linewidth', 2, 'markerfacecolor', 'r'); hold on;
semilogx(isos_1080, dr_1080, '-ob;1080p;', 'linewidth', 2, 'markerfacecolor', 'b');
semilogx(isos_crop, dr_crop, '-og;crop;', 'linewidth', 2, 'markerfacecolor', 'g');
semilogx(isos_photo, dr_photo, '-oy;photo;', 'linewidth', 2, 'markerfacecolor', 'y');

semilogx([50 25600], log2([ref_electrons_5d3*2 ref_electrons_5d3/256]), 'g--;1:1 shot noise limit;')
semilogx([50 25600], log2([ref_electrons_5d3*9*2 ref_electrons_5d3*9/256]), 'b--;3x3 binning shot noise limit;')
semilogx([50 25600], log2([ref_electrons_5d3*15*2 ref_electrons_5d3*15/256]), 'r--;5x3 binning shot noise limit;')

set(gca,'xtick', isos)
set(gca,'xticklabel', isos)
grid on
xlabel('ISO');
ylabel('Dynamic range (EV)');
title('5D3 - movie mode');
legend("location", "southwest")
axis([50 25600 7 12.5]);

offsets = [-5 20 80 150 400 1000 2500 5000];
for i = 1:length(isos)
    label = sprintf("\nISO %d", isos(i));
    if (i > 1), label = sprintf("\n%d", isos(i)); end
    text(isos_1080(i)+offsets(i), dr_max(i)+0.5, label, 'horizontalalignment','center');
end
print -dpng -r70 -FHelvetica:14 5d3-movie.png

isos_1080_5d3 = isos_1080;
dr_1080_5d3 = dr_1080;
isos_photo_5d3 = isos_photo;
dr_photo_5d3 = dr_photo;

# 6D

clf
isos =      [ 100    200    400    800    1600   3200   6400  12800 ];
dr_720 =    [ 11.87  11.85  11.81  11.67  11.40  10.96  10.48 9.35  ];
dr_1080 =   [ 11.84  11.82  11.75  11.62  11.38  10.95  10.47 9.44  ];
dr_crop =   [ 11.78  11.72  11.59  11.29  10.85  10.25  9.65  8.51  ];
dr_photo =  [ 11.48  11.49  11.28  11.00  10.55   9.87  9.11  8.04  ];

dr_max = max(dr_720, dr_1080);
dr_max = max(dr_max, dr_crop);

isos_720    = round(100 ./ ([242833 110182 66885 36324 18721 9748 4874 2380] / ref_electrons_6d / 3))
isos_1080   = round(100 ./ ([274301 142616 75858 39219 19868 10092 4931 2408] / ref_electrons_6d / 3))
isos_crop   = round(100 ./ ([97890 48335 25171 12958 6508 3354 1663 831] / ref_electrons_6d))
isos_photo  = round(100 ./ ([79546 42220 20643 10482 5215 2568 1271 556] / ref_electrons_6d))

semilogx(isos_720, dr_720, '-or;720p;', 'linewidth', 2, 'markerfacecolor', 'r'); hold on;
semilogx(isos_1080, dr_1080, '-ob;1080p;', 'linewidth', 2, 'markerfacecolor', 'b');
semilogx(isos_crop, dr_crop, '-og;crop;', 'linewidth', 2, 'markerfacecolor', 'g');
semilogx(isos_photo, dr_photo, '-oy;photo;', 'linewidth', 2, 'markerfacecolor', 'y');

semilogx([50 25600], log2([ref_electrons_6d*2 ref_electrons_6d/256]), 'g--;1:1 shot noise limit;')
semilogx([50 25600], log2([ref_electrons_6d*3*2 ref_electrons_6d*3/256]), 'b--;1x3 binning shot noise limit;')

set(gca,'xtick', isos)
set(gca,'xticklabel', isos)
grid on
xlabel('ISO');
ylabel('Dynamic range (EV)');
title('6D - movie mode');
legend("location", "southwest")
axis([50 25600 7 12.5]);

offsets = [-5 20 80 150 400 1000 2500 5000];
for i = 1:length(isos)
    label = sprintf("\nISO %d", isos(i));
    if (i > 1), label = sprintf("\n%d", isos(i)); end
    text(isos_1080(i)+offsets(i), dr_max(i)+0.5, label, 'horizontalalignment','center');
end
print -dpng -r70 -FHelvetica:14 6d-movie.png

isos_1080_6d = isos_1080;
dr_1080_6d = dr_1080;
isos_photo_6d = isos_photo;
dr_photo_6d = dr_photo;

# 5D3 vs 6d 1080p normalized
clf

norm_5d3 = log2(sqrt((1920*1080))) - log2(sqrt(1920*1080));
norm_6d = log2(sqrt(1792*1008)) - log2(sqrt(1920*1080));
semilogx(isos_1080_5d3, dr_1080_5d3 + norm_5d3, '-ob;5D3 1080p;', 'linewidth', 2, 'markerfacecolor', 'b'); hold on;
semilogx(isos_1080_6d, dr_1080_6d + norm_6d, '-or;6D 1080p, normalized from 1792x1008;', 'linewidth', 2, 'markerfacecolor', 'r');

semilogx([50 25600], log2([ref_electrons_6d*3*2 ref_electrons_6d*3/256]) + norm_6d, 'r--;6D shot noise limit (1x3 binning);')
semilogx([50 25600], log2([ref_electrons_5d3*9*2 ref_electrons_5d3*9/256]) + norm_5d3, 'b--;5D3 shot noise limit (3x3 binning);')

set(gca,'xtick', isos)
set(gca,'xticklabel', isos)
grid on
xlabel('ISO');
ylabel('Dynamic range (EV)');
title('5D3 vs 6D - 1080p movie mode, normalized to 1920x1080');
legend("location", "southwest")
axis([50 25600 7 12.5]);

dr_max = dr_1080_6d - 0.1

for i = 1:length(isos)
    label = sprintf("\nISO %d", isos(i));
    if (i > 1), label = sprintf("\n%d", isos(i)); end
    text(isos_1080(i)+offsets(i), dr_max(i)+0.5, label, 'horizontalalignment','center');
end

print -dpng -r70 -FHelvetica:14 5d3-vs-6d-movie.png

# 5D3 vs 6d photo normalized
clf
norm_5d3 = log2(sqrt((5796*3870))) - log2(sqrt(8e6));
norm_6d = log2(sqrt(5496*3670)) - log2(sqrt(8e6));
semilogx(isos_photo_5d3, dr_photo_5d3 + norm_5d3, '-ob;5D3 photo;', 'linewidth', 2, 'markerfacecolor', 'b'); hold on;
semilogx(isos_photo_6d, dr_photo_6d + norm_6d, '-or;6D photo;', 'linewidth', 2, 'markerfacecolor', 'r');

semilogx([50 25600], log2([ref_electrons_6d*2 ref_electrons_6d/256]) + norm_6d, 'r--;6D shot noise limit;')
semilogx([50 25600], log2([ref_electrons_5d3*2 ref_electrons_5d3/256]) + norm_5d3, 'b--;5D3 shot noise limit;')

set(gca,'xtick', isos)
set(gca,'xticklabel', isos)
grid on
xlabel('ISO');
ylabel('Dynamic range (EV)');
title('5D3 vs 6D - photo mode, normalized to 8MP');
legend("location", "southwest")
axis([50 25600 7 12.5]);

dr_photo_6d + norm_6d
dr_photo_5d3 + norm_5d3
dr_max = dr_photo_6d + norm_6d - 0.1

for i = 1:length(isos)
    label = sprintf("\nISO %d", isos(i));
    if (i > 1), label = sprintf("\n%d", isos(i)); end
    text(isos_1080(i)+offsets(i), dr_max(i)+0.5, label, 'horizontalalignment','center');
end

print -dpng -r70 -FHelvetica:14 5d3-vs-6d-photo-8mp.png