function C4x = C4coeff(n)
%C4COEFF  Evaluate coefficients for C_4
%
%   C4x = C4COEFF(n) evaluates the coefficients of epsilon^l in expansion
%   of the area (Eq. (65) expressed in terms of n and epsi).  n is a
%   scalar.  C4x is a 1 x 21 array.

persistent coeff nC4 nC4x
if isempty(coeff)
nC4 = 6;
nC4x = (nC4 * (nC4 + 1)) / 2;
coeff = [ ...
97, 15015, ...
1088, 156, 45045, ...
-224, -4784, 1573, 45045, ...
-10656, 14144, -4576, -858, 45045, ...
64, 624, -4576, 6864, -3003, 15015, ...
100, 208, 572, 3432, -12012, 30030, 45045, ...
1, 9009, ...
-2944, 468, 135135, ...
5792, 1040, -1287, 135135, ...
5952, -11648, 9152, -2574, 135135, ...
-64, -624, 4576, -6864, 3003, 135135, ...
8, 10725, ...
1856, -936, 225225, ...
-8448, 4992, -1144, 225225, ...
-1440, 4160, -4576, 1716, 225225, ...
-136, 63063, ...
1024, -208, 105105, ...
3584, -3328, 1144, 315315, ...
-128, 135135, ...
-2560, 832, 405405, ...
128, 99099, ...
];
end
C4x = zeros(1, nC4x);
o = 1;
k = 1;
for l = 0 : nC4 - 1
for j = nC4 - 1 : -1 : l
m = nC4 - j - 1;
C4x(k) = polyval(coeff(o : o + m), n) / coeff(o + m + 1);
k = k + 1;
o = o + m + 2;
end
end
end