function P = polnom(tau1a,tau2a,zetaa,gl,gp,gc,elx,epx,alpha,beta,L)
%
% Return the coeficients of the polynom P
% 
% (c) S. Balac - MAPLY - INSA Lyon - mars 2001
%
% INPUT:
% tau1a = dimensionless tau1 parameter
% tau2a = dimensionless tau2 parameter
% zetaa = dimensionless zeta parameter
% gl,gc,gp = gradients intensity (read-out, slice selection and phase)
% elx = dot product of nl and ex vectors
% epx = dot product of np and ex vectors
% alpha, beta = constants in the magnetic field expression*
% 
% OUTPUT:
%  the polynom P = P(1)*X^5 + ... + P(4)*X + P(5).

ind = 1;
P(0+ind) =-zetaa^5*gc^5*L^3-2*zetaa*gc*L^3*tau1a^2*gl^4*tau2a^2-zetaa*gc*L^3*tau1a^4*gl^4  ...
       -zetaa*gc*L^3*tau2a^4*gl^4+4*zetaa^2*gc^2*L^3*tau1a^3*gl^3  ...
       -6*zetaa^3*gc^3*L^3*tau1a^2*gl^2+4*zetaa^4*gc^4*L^3*tau1a*gl  ...
       -2*zetaa^3*gc^3*L^3*tau2a^2*gl^2+4*zetaa^2*gc^2*L^3*tau1a*gl^3*tau2a^2  ...
       +alpha*L^2*gc^4*zetaa^4+alpha*L^2*tau2a^4*gl^4+2*beta*gl^4*elx^2*tau1a^2  ...
      +2*beta*gl^4*epx^2*tau2a^2-beta*gl^2*gc^2*zetaa^2+2*beta*gl^3*tau1a*gc*zetaa  ...
       -4*beta*gl^3*elx^2*tau1a*gc*zetaa+4*beta*gl^4*elx*tau1a*epx*tau2a  ...
      +2*beta*gl^2*elx^2*gc^2*zetaa^2-4*beta*gl^3*elx*gc*zetaa*epx*tau2a  ...
      -beta*gl^4*tau1a^2  ...
       -beta*gl^4*tau2a^2+alpha*L^2*tau1a^4*gl^4+2*alpha*L^2*tau1a^2*gl^4*tau2a^2  ...
       -4*alpha*L^2*tau1a^3*gl^3*gc*zetaa+6*alpha*L^2*tau1a^2*gl^2*gc^2*zetaa^2  ...
       -4*alpha*L^2*tau1a*gl*gc^3*zetaa^3+2*alpha*L^2*gc^2*zetaa^2*tau2a^2*gl^2  ...
       -4*alpha*L^2*tau1a*gl^3*gc*zetaa*tau2a^2;

P(1+ind) = -2*beta*gl^3*tau1a*gc+2*gc*L^3*tau1a^2*gl^4*tau2a^2+gc*L^3*tau1a^4*gl^4 ...
      +gc*L^3*tau2a^4*gl^4-8*gc^2*L^3*tau1a^3*gl^3*zetaa+18*gc^3*L^3*tau1a^2*gl^2*zetaa^2 ...
      -16*gc^4*L^3*tau1a*gl*zetaa^3+6*gc^3*L^3*zetaa^2*tau2a^2*gl^2 ...
      -8*gc^2*L^3*tau1a*gl^3*zetaa*tau2a^2+5*gc^5*L^3*zetaa^4+2*beta*gl^2*gc^2*zetaa ...
      +4*beta*gl^3*elx^2*tau1a*gc-4*beta*gl^2*elx^2*gc^2*zetaa+4*beta*gl^3*elx*gc*epx*tau2a ...
      -4*alpha*L^2*gc^4*zetaa^3+12*alpha*L^2*tau1a*gl*gc^3*zetaa^2+4*alpha*L^2*tau1a^3*gl^3*gc ...
      -12*alpha*L^2*tau1a^2*gl^2*gc^2*zetaa+4*alpha*L^2*tau1a*gl^3*gc*tau2a^2 ...
      -4*alpha*L^2*gc^2*zetaa*tau2a^2*gl^2;

P(2+ind) = 24*gc^4*L^3*tau1a*gl*zetaa^2+4*gc^2*L^3*tau1a^3*gl^3-18*gc^3*L^3*tau1a^2*gl^2*zetaa ...
      +4*gc^2*L^3*tau1a*gl^3*tau2a^2-6*gc^3*L^3*zetaa*tau2a^2*gl^2-beta*gl^2*gc^2 ...
      -10*gc^5*L^3*zetaa^3+2*beta*gl^2*elx^2*gc^2+6*alpha*L^2*gc^4*zetaa^2 ...
      +6*alpha*L^2*tau1a^2*gl^2*gc^2+2*alpha*L^2*gc^2*tau2a^2*gl^2 ...
      -12*alpha*L^2*tau1a*gl*gc^3*zetaa;

P(3+ind) = 4*alpha*L^2*tau1a*gl*gc^3+6*gc^3*L^3*tau1a^2*gl^2+2*gc^3*L^3*tau2a^2*gl^2 ...
      -16*gc^4*L^3*tau1a*gl*zetaa+10*gc^5*L^3*zetaa^2-4*alpha*L^2*gc^4*zetaa;

P(4+ind) = -5*gc^5*L^3*zetaa+4*gc^4*L^3*tau1a*gl+alpha*L^2*gc^4;

P(5+ind) = L^3*gc^5;

P = fliplr(P); % the matlab convention for writting polynomial is
               % P(1)*X^N + ... + P(N)*X + P(N+1)
