MODULE astronomy
  IMPLICIT NONE
  SAVE

  REAL :: aphelie, periheli, year_day, peri_day, obliquit, &
       timeperi, e_elips,p_elips
  REAL, PARAMETER :: unitastr=149.597927, & ! millions of km
       pi=2.*ASIN(1.)
CONTAINS
 
  SUBROUTINE solarlong(pday,psollong)

!=======================================================================
!      Calcul de la distance soleil-planete et de la declinaison
!   en fonction du jour de l'annee.
!
!   Methode:
!   --------
!      Calcul complet de l'ellipse
!
!   Input:
!   ------
!   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
!   lwrite        clef logique pour sorties de controle
!
!   Output:
!   -------
!   pdist_sol     distance entre le soleil et la planete
!                 ( en unite astronomique pour utiliser la constante 
!                  solaire terrestre 1370 Wm-2 )
!   pdecli        declinaison ( en radians )
!
!=======================================================================

    USE planet
    REAL, INTENT(IN) :: pday
    REAL, INTENT(OUT) :: psollong
    LOGICAL, PARAMETER ::  lwrite=.TRUE.
    REAL pdist_sol, pdecli

! Local:
! ------
    REAL zanom,xref,zx0,zdx,zteta,zz
    INTEGER iter

!--------------------------------------------------------
! calcul de l'angle polaire et de la distance au soleil :
! -------------------------------------------------------

!  calcul de l'zanomalie moyenne

    zz=(pday-peri_day)/year_day
    zanom=2.*pi*(zz-nint(zz))
    xref=abs(zanom)

!  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
!  methode de Newton
    
    zx0=xref+e_elips*sin(xref)
    DO iter=1,10
       zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
       zx0=zx0+zdx
    END DO
    zx0=zx0+zdx
    if(zanom.lt.0.) zx0=-zx0
    
    ! zteta est la longitude solaire
    
    zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
    
    psollong=zteta-timeperi
    
    IF(psollong.LT.0.) psollong=psollong+2.*pi
    IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
    !-----------------------------------------------------------------------
    !   sorties eventuelles:
    !   ---------------------
    
    IF (lwrite) THEN
       PRINT*,'jour de l"annee   :',pday
       PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
       PRINT*,'declinaison (en degres) :',pdecli*180./pi
    ENDIF
    
  END SUBROUTINE solarlong
  
END MODULE astronomy