C4coeff.m

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