function P = polnom(ltau1,ltau2,zeta,lc)
%
% Return the coeficients of the polynomial P

% (c) S. Balac - 1999

global r0
global Z    % coordinates of z in the local frame
global beta 
global g

lgeta = beta/(g*r0.^4);
gcl = 1;

A = Z(3) + gcl*Z(1);
B = Z(1)*ltau1 + Z(2)*ltau2 + Z(3)*lc - Z(1)*gcl*zeta;
C = 1+gcl.^2;
D = 2*gcl*(ltau1-gcl*zeta) + 2*lc;
E = (ltau1-gcl*zeta).^2+ ltau2.^2 +lc.^2;

P(13) =  - lgeta.^2.*E.^2 - 9.*lgeta.^2.*B.^4 + 6.*lgeta.^2.*B.^2.*E ...
      + E.^5.*zeta.^2;
P(12) = 5.*D.*E.^4.*zeta.^2 - 2.*E.^5.*zeta - 36.*lgeta.^2.*A.*B.^3 ...
      + 12.*lgeta.^2.*A.*B.*E + 6.*lgeta.^2.*B.^2.*D - 2.*lgeta.^2.*D.*E;
P(11) = 12.*lgeta.^2.*A.*B.*D - 54.*lgeta.^2.*A.^2.*B.^2 + E.^5 ...
      + 5.*C.*E.^4.*zeta.^2 ...
      + 10.*D.^2.*E.^3.*zeta.^2 + 6.*lgeta.^2.*B.^2.*C - 2.*lgeta.^2.*C.*E ...
      - lgeta.^2.*D.^2 + 6.*lgeta.^2.*A.^2.*E - 10.*D.*E.^4.*zeta;
P(10) = 20.*C.*D.*E.^3.*zeta.^2 - 20.*D.^2.*E.^3.*zeta - 10.*C.*E.^4.*zeta ...
      - 36.*lgeta.^2.*A.^3.*B + 12.*lgeta.^2.*A.*B.*C - 2.*lgeta.^2.*C.*D ...
      + 6.*lgeta.^2.*A.^2.*D + 5.*D.*E.^4 + 10.*D.^3.*E.^2.*zeta.^2;
P(9) = 6.*lgeta.^2.*A.^2.*C - 20.*D.^3.*E.^2.*zeta - lgeta.^2.*C.^2 ...
     + 10.*D.^2.*E.^3 + 10.*C.^2.*E.^3.*zeta.^2 ...
     + 30.*C.*D.^2.*E.^2.*zeta.^2 + 5.*D.^4.*E.*zeta.^2 ...
     + 5.*C.*E.^4 - 40.*C.*D.*E.^3.*zeta - 9.*lgeta.^2.*A.^4;
P(8) = 20.*C.*D.^3.*E.*zeta.^2 + 10.*D.^3.*E.^2 - 10.*D.^4.*E.*zeta ...
     + 30.*C.^2.*D.*E.^2.*zeta.^2 + D.^5.*zeta.^2 - 20.*C.^2.*E.^3.*zeta ...
     - 60.*C.*D.^2.*E.^2.*zeta + 20.*C.*D.*E.^3;
P(7) = 5.*D.^4.*E - 2.*D.^5.*zeta + 30.*C.*D.^2.*E.^2 + 10.*C.^2.*E.^3 ...
     - 40.*C.*D.^3.*E.*zeta + 30.*C.^2.*D.^2.*E.*zeta.^2 ...
     - 60.*C.^2.*D.*E.^2.*zeta + 5.*C.*D.^4.*zeta.^2  ...
     + 10.*C.^3.*E.^2.*zeta.^2;
P(6) = - 60.*C.^2.*D.^2.*E.*zeta + D.^5 + 20.*C.^3.*D.*E.*zeta.^2 ...
     - 10.*C.*D.^4.*zeta + 20.*C.*D.^3.*E + 10.*C.^2.*D.^3.*zeta.^2 ...
     + 30.*C.^2.*D.*E.^2 - 20.*C.^3.*E.^2.*zeta;
P(5) = 10.*C.^3.*E.^2 - 20.*C.^2.*D.^3.*zeta + 30.*C.^2.*D.^2.*E ...
     + 10.*C.^3.*D.^2.*zeta.^2 + 5.*C.*D.^4 - 40.*C.^3.*D.*E.*zeta ...
     + 5.*C.^4.*E.*zeta.^2;
P(4) =  - 20.*C.^3.*D.^2.*zeta + 10.*C.^2.*D.^3 + 20.*C.^3.*D.*E ...
     - 10.*C.^4.*E.*zeta + 5.*C.^4.*D.*zeta.^2;
P(3) = - 10.*C.^4.*D.*zeta + 10.*C.^3.*D.^2 + 5.*C.^4.*E + C.^5.*zeta.^2;
P(2) = - 2.*C.^5.*zeta + 5.*C.^4.*D;
P(1) = C.^5;
