function ret = p21(nk, nm, np)
ret = (1 - (1 - 0.94)^nk) * (1 - (1 - (1 - (1 - 0.93)^np) * (1 - (1 - 0.93)^nm))^3);
endfunction

function ret = c21(nk, nm, np)
ret = 3 * np * 15 + 3 * nm * 6 + nk * 2;
endfunction

function ret = m21(nk, nm, np)
ret = 1 / (0.2 - (0.7) / (23 * nk)) + 1 / ((1/23) - 0.7 / (23 * 3 * np)) + 1 / ((1/13) - 0.7 / (23 * 3 * nm));
endfunction

function ret = lam(nk, nm, np)
ret = min([nk / 5, 3 * np / 23, 3 * nm / 13]);
endfunction

function ret = table(nk, nm, np)
ret = [p21(nk, nm, np), m21(nk, nm, np), c21(nk, nm, np), lam(nk, nm, np)];
endfunction

for nk = 1:93
    for nm = 1:11
        for np = 1:5
            a = table(nk, nm, np);
            if (1 - a(1) < 1e-9 && a(2) <= 43.4 && a(3) <= 250)
                [nk, nm, np]
                [1 - a(1), a(2), a(3), a(4)]
            endif
        endfor
    endfor
endfor

for nk = 1:93
    for nm = 1:11
        for np = 1:5
            a = table(nk, nm, np);
            if (1 - a(1) < 1e-9 && a(2) <= 43.4 && a(3) <= 250)
                [nk, nm, np]
                [1 - (1 - a(1)) / 1e-9, 1 - (a(2) - 41) / (43.4 - 41), 1 - (a(3) - 240) / 10, (a(4) - 0.35) / (0.6 - 0.35)]
            endif
        endfor
    endfor
endfor

