program test_figures
!     Version: 00/00/03
! Test des subroutines figures3 et figures4 a partir des resultats 
! analytiques de Hubbard et al. Icarus (1974).

      implicit none
      integer n
      parameter (n=201)
      integer ndis,idisc(1),i
      real*8 mplanet,rplanet,omplanet,a1_2,j2_2,j4_2,j6_2,              &
     &     a1_3,j2_3,j4_3,j6_3,a1_4,j2_4,j4_4,j6_4,j8_4,pir,romoy,alpha
      real*8 m(n),r(n),ro(n),dror(n),deltaro(1),b(n),h(n)
      real*8 req_jup
      include './cepam.inc'
      include './Communs/ctephys.inc'

!----------------------------------------------------------------------
 2000 format(1p,6d12.5)
 2001 format(70('*'),/,a,' model from ',a,/,70('*'))
 2003 format(/,t20,'Req',t32,'J2',t44,'J4',t56,'J6',/,                  &
     &     'Reference',t19,1p,4d12.5,/,'Calculated',t19,4d12.5,/)
 2004 format(/,t20,'Req',t32,'J2',t44,'J4',t56,'J6',t68,'J8',/,         &
     &     'Reference',t19,1p,5d12.5,/,'Calculated',t19,5d12.5,/)
 2234 format(/,t20,'Req',t32,'J2',t44,'J4',t56,'J6',/,                  &
     &     'J3/J2-1  ',t19,1p,3d12.5,/,'J4/J3-1   ',t19,4d12.5,/)
      call constantes(' ')

      write(*,2001)'Polytropic n=1','Hubbard 1974'
      alpha=4.7
!      mplanet=1.9d30
!      rplanet=(mplanet*pi/4/alpha)**(1/3.d0)
      rplanet=6.98d9
      mplanet=4*alpha/pi*rplanet**3
      req_jup=71400d5
      omplanet=sqrt(0.0888*g*mplanet/req_jup**3)
      write(*,*)'Omega=',omplanet
      do i=1,n-1
       r(i)=(n-i)/(n-1.d0)
       pir=pi*r(i)
       ro(i)=alpha*sin(pir)/(pir)
       m(i)=1./pi*(sin(pir)-pi*r(i)*cos(pir))
       dror(i)=alpha*pi*(cos(pir)/pir-sin(pir)/pir/pir)
!       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
      enddo
      r(n)=0.d0
      m(n)=0.d0
      ro(n)=alpha
      dror(n)=0.d0
      ndis=1
      idisc(1)=1
      deltaro(1)=0

!      call figures2(mplanet,rplanet,omplanet,m,r,ro,dror,
!     &     n,ndis,idisc,deltaro,a1_2,j2_2,j4_2,b,h)
!      write(*,2000)a1_2,j2_2,j4_2

      call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
      j2_3=j2_3/(7.14d9/a1_3)**2
      j4_3=j4_3/(7.14d9/a1_3)**4
      j6_3=j6_3/(7.14d9/a1_3)**6
      write(*,2003)7.14d9,0.0139,-0.00052,0.000039,a1_3,j2_3,j4_3,j6_3

      call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
      write(*,2004)7.14d9,0.0139,-0.00052,0.000039,0.d0,                &
     &     a1_4,j2_4,j4_4,j6_4,j8_4

!----------------------------------------------------------------------
      write(*,2001)'Linear','Hubbard 1974'
      rplanet=6.97d9
      mplanet=5.04/3.*pi*rplanet**3
      do i=1,n
       m(i)=r(i)**3*(4-3*r(i))
       ro(i)=5.04*(1-r(i))
       dror(i)=-5.04
!       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
      enddo

      call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
      write(*,2003)7.14d9,0.0154,-0.00064,0.000039,a1_3,j2_3,j4_3,j6_3

      call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
      write(*,2004)7.14d9,0.0154,-0.00064,0.000039,0.d0,                &
     &     a1_4,j2_4,j4_4,j6_4,j8_4

!----------------------------------------------------------------------
!       write(*,"(70('*'))")
!       do i=1,n
!        m(i)=r(i)**3*(5-3*r(i)*r(i))/2.
!        ro(i)=3.34*(1-r(i)*r(i))
!        dror(i)=-3.34*2*r(i)
! c       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
!       enddo

!       write(*,2000)7.14d9,0.0181,-0.00083,0.000054
!       call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,
!      &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
!       write(*,2000)a1_3,j2_3,j4_3,j6_3

!       call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,
!      &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
!       write(*,2000)a1_4,j2_4,j4_4,j6_4,j8_4

!----------------------------------------------------------------------
      write(*,2001)'Jupiter Linear','Zharkov & Trubytsyn 1978 p.275'
      mplanet=1.9d30
      rplanet=6.98d9
      omplanet=sqrt(0.0830*g*mplanet/rplanet**3)
      write(*,*)omplanet
      romoy=3*mplanet/4./pi/rplanet**3 !romoyen
      do i=1,n
       m(i)=r(i)**3*(4-3*r(i))
       ro(i)=4*(1-r(i))*romoy
       dror(i)=-4*romoy
!       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
      enddo

      write(*,2000)

      call figures2(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_2,j2_2,j4_2,b,h)
!      write(*,2000)a1_2,j2_2,j4_2

      call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
      write(*,2003)7.140d9,1.4798d-2,-5.877d-4,3.528d-5,                &
     &     a1_3,j2_3,j4_3,j6_3

      call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
      write(*,2004)7.140d9,1.4798d-2,-5.929d-4,3.457d-5,-2.77d-6,       &
     &     a1_4,j2_4,j4_4,j6_4,j8_4

      write(*,2234)a1_3/a1_2-1,j2_3/j2_2-1,j4_3/j4_2-1,                 &
     &     a1_4/a1_3-1,j2_4/j2_3-1,j4_4/j4_3-1,j6_4/j6_3-1

!----------------------------------------------------------------------
      write(*,2001)'Saturn Linear','Zharkov & Trubytsyn 1978 p.275'
      mplanet=0.569d30
      rplanet=5.78d9
      omplanet=sqrt(0.143*g*mplanet/rplanet**3)
      write(*,*)omplanet
      romoy=3*mplanet/4./pi/rplanet**3 !romoyen
      do i=1,n
       m(i)=r(i)**3*(4-3*r(i))
       ro(i)=4*(1-r(i))*romoy
       dror(i)=-4*romoy
!       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
      enddo

      call figures2(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_2,j2_2,j4_2,b,h)
!     write(*,2000)a1_3,j2_3,j4_3

      call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
      write(*,2003)6.020d9,2.512d-2,-1.665d-3,1.726d-4,                 &
     &     a1_3,j2_3,j4_3,j6_3

      call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
      write(*,2004)6.020d9,2.510d-2,-1.713d-3,1.658d-4,-2.43d-5,        &
     &     a1_4,j2_4,j4_4,j6_4,j8_4

      write(*,2234)a1_3/a1_2-1,j2_3/j2_2-1,j4_3/j4_2-1,                 &
     &     a1_4/a1_3-1,j2_4/j2_3-1,j4_4/j4_3-1,j6_4/j6_3-1

!----------------------------------------------------------------------
      write(*,2001)'Saturn Polytropic n=1',                             &
     &     'Zharkov & Trubytsyn 1978 p.291'
      mplanet=0.569d30
      omplanet=1.677d-4
      rplanet=(0.141*g*mplanet/omplanet/omplanet)**(1/3.d0)
      req_jup=59800.d5
      romoy=3*mplanet/4/pi/rplanet**3
      write(*,*)'rplanet=',rplanet,'romoy=',romoy,'= 0.7142?'
      alpha=pi*pi*romoy/3.
      write(*,*)'alpha=',alpha
      do i=1,n-1
       pir=pi*r(i)
       ro(i)=alpha*sin(pir)/(pir)
       m(i)=1./pi*(sin(pir)-pi*r(i)*cos(pir))
       dror(i)=alpha*pi*(cos(pir)/pir-sin(pir)/pir/pir)
!       write(*,'(1p,4d12.5)')r(i),m(i),ro(i),dror(i)
      enddo
      r(n)=0.d0
      m(n)=0.d0
      ro(n)=4.7
      dror(n)=0.d0
      ndis=1
      idisc(1)=1
      deltaro(1)=0

!      call figures2(mplanet,rplanet,omplanet,m,r,ro,dror,
!     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,b,h)
!      write(*,2000)a1_3,j2_3,j4_3

      call figures3(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_3,j2_3,j4_3,j6_3,b,h)
      write(*,2003)5.98d9,2.2324d-2,-1.4755d-3,1.399d-4,                &
     &     a1_3,j2_3,j4_3,j6_3

      call figures4(mplanet,rplanet,omplanet,m,r,ro,dror,               &
     &     n,ndis,idisc,deltaro,a1_4,j2_4,j4_4,j6_4,j8_4,b,h)
      write(*,2004)5.98d9,2.2324d-2,-1.4755d-3,1.399d-4,-1.676d-5,      &
     &     a1_4,j2_4,j4_4,j6_4,j8_4

      end