!*******************************************************
      
      program test_poly
!     Version: 16-Sep-2010
      
      implicit none
      include './cepam.inc'
      integer pn_d,pn_f
      parameter (pn_d=pn+pn_noy+pndis,pn_f=pn+pn_noy)
      integer i,ierr
      integer*4 leng
      character*1 oui
! Pour test polytropes
      real*8 qq,mm,j2hub,j4hub,j6hub
!------------------------------ Communs ----------------------------
      real*8 req,j2,j4,j6,j8,nu0
      common/gravmmts/req,j2,j4,j6,j8,nu0
      include './Communs/ctephys.inc'
      include './Communs/planetes.inc'
      include './Communs/modelp.inc'
      integer nd_f,nx_f,m_f,nr_f,knot_f,id_f(0:pndis+1)
      real*8 x_f(pn_f),xt_f(pn_f+pmc_S*(pndis+2)),xr_f(pn_f),           &
     &     s_f(pn_d*pmc_S*2),pot_f(2,pn_d),rplanet
      common/figures_i/nd_f,nx_f,m_f,nr_f,knot_f,id_f
      common/figures_r/x_f,xt_f,xr_f,pot_f,s_f,rplanet
!----------------------------------------------------------------------

      write(*,101)
 101  format(30('*'),' C E P A M ',30('*'),//,                          &
     &     'Programme test des subroutines de rotation',/,              &
     &     ' (utiliser "etat_poly")',/,                                 &
     &     'Indiquer le nom du modele de depart (ex: jup)')
      read(*,'(a)')modele

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     INITIALISATION DE CEPAM                                          c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      call cepam(1,ierr)

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     On ouvre un fichier <modele>.poly                                c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      open(unit=76,file='test_poly.dat',                                &
     &     form='formatted')
      write(76,*)'C TEST/CEPAM & THEORIE DES FIGURES: '
      write(76,*)'C Masse utilisee: [g] ',mplanet*mjup
      write(76,755)

 755   format('c',t5,'q',t17,'m',t25,'Omega [s-1]',t38,'Req [cm]',t49,  &
     &      'Rmoy [cm]',t61,'J2,J4,J6,J8,J2H,J4H,J6H...')
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     1ere boucle: theorie a l'ordre 3                                 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      f_rot='figures3'
      write(76,*)'C THEORIE: ',f_rot
      write(76,*)15
      do i=1,15
       mm=(i/33.d0)**2
       omega=sqrt(mm/(rplanet**3/g/mplanet/mjup))/omjup
!       omega=i/10.d0
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     CALL CEPAM                                                       c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       call cepam(2,ierr)

       qq=(omega*omjup)**2*req**3/g/mplanet/mjup
!       mm=(omega*omjup)**2*rplanet**3/g/mplanet/mjup
       j2hub=0.173273*qq-0.197027*qq**2+0.15*qq**3
       j4hub=-0.081092*qq**2+0.15*qq**3
       j6hub=0.056329*qq**3
       
       write(76,756)qq,mm,omega*omjup,req,rplanet,j2,j4,j6,j8,j2hub,    &
     &      j4hub,j6hub
       write(*,755)
       write(*,756)qq,mm,omega*omjup,req,rplanet,j2,j4,j6,j8,j2hub,     &
     &      j4hub,j6hub

 756   format(1p,6d12.5,7d14.6)
      enddo

!      close(unit=76)
!      write(*,*)'J''ai cree le fichier test_poly.dat'
!      stop

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     2eme boucle: theorie a l'ordre 4                                 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      f_rot='figures4b'
      write(76,*)'C THEORIE: ',f_rot
      write(76,*)15
      do i=1,15
       mm=(i/33.d0)**2
       omega=sqrt(mm/(rplanet**3/g/mplanet/mjup))/omjup
!       omega=i/10.d0
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     CALL CEPAM                                                       c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       call cepam(2,ierr)

       qq=(omega*omjup)**2*req**3/g/mplanet/mjup
!       mm=(omega*omjup)**2*rplanet**3/g/mplanet/mjup
       j2hub=0.173273*qq-0.197027*qq**2+0.15*qq**3
       j4hub=-0.081092*qq**2+0.15*qq**3
       j6hub=0.056329*qq**3
       
       write(76,756)qq,mm,omega*omjup,req,rplanet,j2,j4,j6,j8,j2hub,    &
     &      j4hub,j6hub
       write(*,755)
       write(*,756)qq,mm,omega*omjup,req,rplanet,j2,j4,j6,j8,j2hub,     &
     &      j4hub,j6hub

      enddo

      close(unit=76)
      write(*,*)'J''ai cree le fichier test_poly.dat'

      end