function [ PHISOL ] = calculgram(Nbtermes,NPT,T,dt,V,z0,dz0,alpha0,dalpha0,beta0)
m=1;
ro=1.;
J0=1;
k=1;
c=2;
a0=1.;
b0=1;
a=-.01;
S=1.;
L=.2;
%========
% Construction de M K et F
%========
M=[m, a*cos(alpha0);a*cos(alpha0),J0];
K=[k ,-ro*S*V^2*dcz0(alpha0)/2;0 c-ro*S*L*dcm0(alpha0)/2];
Rad=linsolve(M,K');
F0=[ro*S*V^2*cz0(alpha0)/2,ro*S*L*V^2*cm0(alpha0)/2];
B0=[ro*S*V^2*czg(beta0)/2,ro*S*L*V^2*cmg(beta0)/2];
C0=[ro*S*V^2*dczg(beta0)/2,ro*S*L*V^2*dcmg(beta0)/2];
C=linsolve(M,C0');
B=linsolve(M,B0');
F=linsolve(M,F0');
%==================
% donnees initiales
%==================
% construction de la matrice de Gram
GRAMM=zeros(4,4);
%=======
phi0=linspace(1,4,4);
deltaphi0=linspace(1,4,4);
Lsec=linspace(1,2,4)*0.;
     P1dot=linspace(1,2,NPT);
     P2dot=linspace(1,2,NPT);
     Q1dot=linspace(1,2,NPT);
     Q2dot=linspace(1,2,NPT);
%===========
for k1=1:4
  for i1=1:4
        if(i1 == k1)
                phi0(i1)=1;
        else
                phi0(i1)=0;
        end
  end
    [P1,P2]=soluedo(Rad,NPT,dt,phi0);
     P1dot(1)=phi0(2);
     P2dot(1)=phi0(4);
    
     for i=2:NPT
         P1dot(i)=(P1(i)-P1(i-1))/dt;
         P2dot(i)=(P2(i)-P2(i-1))/dt;
     end
     Lsec(k1)=dz0*phi0(1)+dalpha0*phi0(3)-z0*phi0(2)-alpha0*phi0(4);
     squatro=0;
     for jkl=1:NPT
         squatro=squatro+F(1)*P1(jkl)+F(2)*P2(jkl);
     end
     Lsec(k1)=Lsec(k1)+squatro*dt;
  %==============
  for k2=k1:4
       for i2=1:4
            if(i2 == k2)
               deltaphi0(i2)=1;
            else
               deltaphi0(i2)=0;
            end
       end
%======
% Remplissage des solution elementaire
%======
    [Q1,Q2]=soluedo(Rad,NPT,dt,deltaphi0);
     Q1dot(1)=deltaphi0(2);
     Q2dot(1)=deltaphi0(4);
     for i=2:NPT   
      Q1dot(i)=(Q1(i)-Q1(i-1))/dt;
      Q2dot(i)=(Q2(i)-Q2(i-1))/dt;
     end
%========   
% calcul des termes de GRAM
%========
  
    for n=1:Nbtermes
        coef=dt^2/(a0+b0*n^2*pi^2/T^2);
        sbis=0;
        sters=0;
        for jj=1:NPT
        sinusnjj=sin(n*pi*jj*dt/T);
        sbis=sbis+sinusnjj*(C(1)*P1dot(jj)+C(2)*P2dot(jj)-B(1)*P1(jj)-B(2)*P2(jj));
        sters=sters+sinusnjj*(C(1)*Q1dot(jj)+C(2)*Q2dot(jj)-B(1)*Q1(jj)-B(2)*Q2(jj));
        end
        
    end
    GRAMM(k1,k2)=GRAMM(k1,k2)+sbis*sters*coef;
    GRAMM(k2,k1)=GRAMM(k1,k2);
 end
end
GRAMM=(2/T)*GRAMM;
%======
%  calcul de PHI
%======
PHISOL=linsolve(GRAMM,Lsec');


end