MODULE planete_h

implicit none

REAL,SAVE :: aphelie   ! Aphelion, in Mkm
REAL,SAVE :: periheli  ! Perihelion, in Mkm
REAL,SAVE :: year_day  ! Number of days in the year
!$OMP THREADPRIVATE(aphelie,periheli,year_day)
REAL,SAVE :: peri_day  ! Date of perihelion, in days
REAL,SAVE :: obliquit  ! Obliquity of the planet, in degrees
REAL,SAVE :: lmixmin
!$OMP THREADPRIVATE(peri_day,obliquit,lmixmin)
REAL,SAVE :: emin_turb
REAL,SAVE :: coefvis
REAL,SAVE :: coefir
!$OMP THREADPRIVATE(emin_turb,coefvis,coefir)
REAL,SAVE :: lsperi    ! Solar longitude of the perihelion, angle in rad
REAL,SAVE :: e_elips   ! Orbit eccentricity
REAL,SAVE :: p_elips   ! Ellipse parameter (semi-latus rectum)
!$OMP THREADPRIVATE(lsperi,e_elips,p_elips)
REAL,PARAMETER :: unitastr=149.597927 ! Astronomical unit AU, in Mkm

! NB: Ideally all module variables should be "protected"...

CONTAINS
!-----------------------------------------------------------------------
! Gives sol at perihelion for ls at perihelion (for precession cycles)
  SUBROUTINE lsp2solp(lsp,solp)

  use comcstfi_h, only: pi

  implicit none

  ! Arguments:
  real, intent(in) :: lsp ! ls at perihelion
  real, intent(out) :: solp ! sol at perihelion

  ! Locals:
  real :: zx0 ! eccentric anomaly at Ls=0
  real, parameter :: degrad = 180.d0/pi

  ! Compute orbit geometrical parameters
  e_elips = (aphelie - periheli)/(aphelie + periheli)
  zx0 = -2.0*atan(tan(0.5*lsp/degrad)*sqrt((1. - e_elips)/(1. + e_elips)))
  if (zx0 <= 0.) zx0 = zx0 + 2.*pi

  ! Compute sol at perihelion
  solp = year_day*(1. - (zx0 - e_elips*sin(zx0))/(2.*pi))

  END SUBROUTINE lsp2solp

!-----------------------------------------------------------------------
! Initialize orbital parameters
  SUBROUTINE iniorbit()

  use comcstfi_h, only: pi

  implicit none

  ! Locals:
  character(8), parameter :: myname = "iniorbit" 
  real    :: zxref, zanom, zz, zx0, zdx
  integer :: iter

  ! Info about orbital values
  write(*,*)myname,': Aphelion in Mkm                 =',aphelie
  write(*,*)myname,': Perihelion in Mkm               =',periheli
  write(*,*)myname,': Number of days in the year      =',year_day
  write(*,*)myname,': Date of perihelion in days      =',peri_day
  write(*,*)myname,': Obliquity in degrees            =',obliquit

  ! Compute orbit geometrical parameters
  e_elips = (aphelie - periheli)/(aphelie + periheli)
  p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
  write(*,*)myname,': Orbit eccentricity              =',e_elips
  write(*,*)myname,': Ellipse semi-latus rectum       =',p_elips
      
  ! Compute mean anomaly zanom
  zz = (year_day - peri_day)/year_day
  zanom = 2.*pi*(zz - nint(zz))
  zxref = abs(zanom)
  write(*,*)myname,': zanom                           =',zanom

  ! Solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
  ! using Newton method

  zx0 = zxref + e_elips*sin(zxref)
  DO iter = 1,100
    zdx = -(zx0 - e_elips*sin(zx0) - zxref)/(1. - e_elips*cos(zx0))
    if (abs(zdx).le.(1.e-12)) exit
    zx0 = zx0 + zdx
  ENDDO

  zx0 = zx0 + zdx
  if (zanom < 0.) zx0 = -zx0
  write(*,*)myname,': Eccentric anomaly zx0           =',zx0

  lsperi = -2.*atan(sqrt((1. + e_elips)/(1. - e_elips))*tan(zx0/2.))
  lsperi = modulo(lsperi,2.*pi)
  write(*,*)myname,': Perihelion solar long. Ls (deg) =', lsperi*180./pi

  END SUBROUTINE iniorbit

END MODULE planete_h