1function C4x = C4coeff(n)
2%C4COEFF  Evaluate coefficients for C_4
3%
4%   C4x = C4COEFF(n) evaluates the coefficients of epsilon^l in expansion
5%   of the area (Eq. (65) expressed in terms of n and epsi).  n is a
6%   scalar.  C4x is a 1 x 21 array.
7
8  persistent coeff nC4 nC4x
9  if isempty(coeff)
10    nC4 = 6;
11    nC4x = (nC4 * (nC4 + 1)) / 2;
12    coeff = [ ...
13        97, 15015, ...
14        1088, 156, 45045, ...
15        -224, -4784, 1573, 45045, ...
16        -10656, 14144, -4576, -858, 45045, ...
17        64, 624, -4576, 6864, -3003, 15015, ...
18        100, 208, 572, 3432, -12012, 30030, 45045, ...
19        1, 9009, ...
20        -2944, 468, 135135, ...
21        5792, 1040, -1287, 135135, ...
22        5952, -11648, 9152, -2574, 135135, ...
23        -64, -624, 4576, -6864, 3003, 135135, ...
24        8, 10725, ...
25        1856, -936, 225225, ...
26        -8448, 4992, -1144, 225225, ...
27        -1440, 4160, -4576, 1716, 225225, ...
28        -136, 63063, ...
29        1024, -208, 105105, ...
30        3584, -3328, 1144, 315315, ...
31        -128, 135135, ...
32        -2560, 832, 405405, ...
33        128, 99099, ...
34            ];
35  end
36  C4x = zeros(1, nC4x);
37  o = 1;
38  k = 1;
39  for l = 0 : nC4 - 1
40    for j = nC4 - 1 : -1 : l
41      m = nC4 - j - 1;
42      C4x(k) = polyval(coeff(o : o + m), n) / coeff(o + m + 1);
43      k = k + 1;
44      o = o + m + 2;
45    end
46  end
47end
48