# -------- Initialize ----------
isochrones = rand (73160, 29);
j(1, 1, 1, :) = isochrones(:, 15);
k(1, 1, 1, :) = isochrones(:, 17);
ukidss = rand (2, 77);
oj(1, 1, :, 1) = ukidss(:, 40);
ok(1, 1, :, 1) = ukidss(:, 76);
sj(1, 1, :, 1) = ukidss(:, 41);
sk(1, 1, :, 1) = ukidss(:, 77);
e = 0.1/2:0.1:19.75;
KJ = 0.2833; KK = 0.11424;
d = (0:0.1:12)';
m = convn (d, [1; 1]/2, 'valid');
md =  5*log10 (m*1000) - 5;

# -------- The processing ----------
# The maximum number of elements in 32-bit addressing is of order 268402689
# hence split the arrays by splitting the isochrones array by N (~15)
N = 1;
t = floor (linspace (1, rows (isochrones), N + 1)); t(end)++;
for i = 1:columns (t) - 1
  jk_C{i, 1} = j(1, 1, 1, t(i):t(i + 1) - 1);
  jk_C{i, 2} = k(1, 1, 1, t(i):t(i + 1) - 1);
endfor;
# Replicate the other variables
[md_C{1:N, 1:2}] = deal (md);
[e_C{1:N, 1:2}] = deal (e);
[K_C{1:N, 1}] = deal (KJ);
[K_C{1:N, 2}] = deal (KK);
[sjk_C{1:N, 1}] = deal (sj);
[sjk_C{1:N, 2}] = deal (sk);
[ojk_C{1:N, 1}] = deal (oj);
[ojk_C{1:N, 2}] = deal (ok);

# Perform the parallel computation without the need for 64-bit addressing
x_C = parcellfun (N*2, @exp_arg, jk_C, md_C, e_C, K_C, sjk_C, ojk_C,
                  "UniformOutput", false);
