      program Mercury ! BB June 01 2017.
c Calculation of angular momentum and other parameters.
c Gravitational units G=c=1.
      implicit real*8 (a-h,o-z)
      logical ltest
      parameter(iuo=6)
      parameter(c=2.998d+08) ! 2.998d+08 [m/s] = 1 -> c.
      parameter(G=6.674d-11) ! 6.674d-11 [m^3 kg^{-1} s^{-2}] = 1 -> G.
      parameter(y=365.d00*24.d00*3600.d0)  ! Year in [s].
      parameter(ym=116.d00*24.d00*3600.d0) ! Mercury year in [s].
      parameter(r=70006.D06) ! Mercury distance to sun [m].
      parameter(sm=1.989d30) ! Sun mass [kg].
      ltest=.false.
      if(ltest) stop "Mercury: ltest."
      ltest=.true.
      write(iuo,'(/," sm [kg], r [m] =",2g12.4)') sm,r
c
      write(iuo,'(/," In units of seconds:")') 
      s=1.d00
      xm=s/c            ! meter    in seconds.
      xkg=G*xm**3/s**2  ! kilogram in seconds.
      write(iuo,'(/,"  m [s], kg [s] =",2g12.4)') xm,xkg
c
      smsec=sm*xkg
      rsec=r*xm
      write(iuo,'(" sm [s],  r [s] =",2g12.4)') smsec,rsec
      omega=sqrt(smsec/rsec**3)
      write(iuo,'(" omega [1/s]    =",1g12.4)') omega
      beta=rsec*omega
      phiyear=omega*ym
      write(iuo,'(" phiyear [rad]  =",1g12.4)') phiyear
      tpi=2*acos(-1.d00)
      phiyear=360*phiyear/tpi
      write(iuo,'(" phiyear [deg]  =",1g12.4,/)') phiyear
c
      write(iuo,'(" velocity beta  =",1g12.4)') beta
      h=rsec*beta
      write(iuo,'(" angular mom h  =",1g12.4)') h
      write(iuo,'(" Contributions: ")')
      xn1=2*smsec/rsec
      write(iuo,'(" Newton 1       =",1g12.4)') xn1
      xn2=(h/rsec)**2
      write(iuo,'(" Newton 2       =",1g12.4)') xn2
      xn=xn1-xn2
      write(iuo,'(" Newton         =",1g12.4)') xn
      s=2*smsec*h**2/rsec**3 ! Schwarzschild contribution.
      write(iuo,'(" GR term        =",1g12.4)') s
      write(iuo,'(" rsec/smsec     =",1g12.4)') rsec/smsec
      write(iuo,'(" Ratios:        =",3g12.4,/)') s/xn1,s/xn2,s/xn
c      
      write(iuo,'(" In units of meter:",/)') 
      xm=1.d00
      sec=xm*c
      xkg=G*xm**3/sec**2
      smm=sm*xkg
      write(iuo,'(" sm [m], r [m]  =",2g12.4)') smm,r
      h=r*beta
      write(iuo,'(" angular mom h  =",1g12.4)') h
      write(iuo,'(" Contributions: ")')
      xn1=2*smm/r
      write(iuo,'(" Newton 1       =",1g13.5)') xn1
      xn2=(h/r)**2
      write(iuo,'(" Newton 2       =",1g13.5)') xn2
      xn=xn1-xn2
      write(iuo,'(" Newton         =",1g13.5)') xn
      s=2*smm*h**2/r**3 ! Schwarzschild contribution.
      write(iuo,'(" GR term        =",1g13.5)') s
      write(iuo,'(" r/smm          =",1g12.4)') r/smm
      write(iuo,'(" Ratios:        =",3g13.5,/)') s/xn1,s/xn2,s/xn
c
      stop "Mercury: all done."
      end
