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 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
      
  SUBROUTINE iniorbit(paphelie,pperiheli,pyear_day,pperi_day,pobliq)
  use comcstfi_h, only: pi
  IMPLICIT NONE
  ! initialize orbital parameters
  
  REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq

  CHARACTER(LEN=8),PARAMETER :: myname="iniorbit" 
  REAL :: zxref,zanom,zz,zx0,zdx
  INTEGER :: iter

  ! copy over input values
  aphelie =paphelie
  periheli=pperiheli
  year_day=pyear_day
  obliquit=pobliq
  peri_day=pperi_day
  
  write(*,*)myname,': Perihelion in Mkm  ',periheli
  write(*,*)myname,': Aphelion  in Mkm  ',aphelie
  write(*,*)myname,': Obliquity in degrees  :',obliquit
      
  ! compute orbit geometrical parameters
  e_elips=(aphelie-periheli)/(periheli+aphelie)
  p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
      
  write(*,*)myname,': e_elips',e_elips
  write(*,*)myname,': p_elips',p_elips
      
  !  compute mean anomaly zanom

  zz=(year_day-pperi_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.lt.0.) zx0=-zx0
  write(*,*)myname,': 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