SUBROUTINE orbite(pls,pdist_sol,pdecli)
      use planet_h          
      IMPLICIT NONE

c=======================================================================
c
c   Objet:
c   ------
c
c   Distance from sun and declimation as a function of the solar
c   longitude Ls
c
c   Interface:
c   ----------
c
c
c
c   Arguments:
c   ----------
c
c   Input:
c   ------
c   pls          Ls
c
c   Output:
c   -------
c   pdist_sol     Distance Sun-Planet in UA
c   pdecli        declinaison ( en radians )
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

!#include "planete.h"
#include "comcstfi.h"

c arguments:
c ----------

      REAL pday,pdist_sol,pdecli,pls,i

c      REAL zanneedebut,zan,zalpha0,zalpha1,zbeta0,zbeta1,zint 
c      REAL zanneereference 

c-----------------------------------------------------------------------
c Initialisations

c      zalpha0 = 76.312*3.14159/180.
c      zalpha1 = 0.0381
c      zalpha1 = 0.025
c      zbeta0 = 18.503*3.15149/180.
c      zbeta1 = 0.00986
c      zanneedebut = 1985. 
c      zanneereference = 1820.
c      daysec=551837.

c      p_elips : parametre p de l'ellispe

c Distance Sun-Planet
c equation polaire ellipse

      pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi))

c Solar declination

c  transformation de pday en annees terrestres (zday=0 
c  correspond au premier janvier 0 A.D.)
c  31557595 est la duree approx d'une annee terrestre en secondes
c  pour le rest de la formule magique cf these benchkoura p110
c  et pour le 569 : c'est un fit sur la courbe dans l'article
c  de Yelle 1995 
      
c      zan  = zanneedebut+(pday*daysec/31557595.)-zanneereference
c      zint = 0.436 * cos(zalpha0 - zalpha1 * zan)
c      zint = zint +0.3495*sin(zalpha0-zbeta0+(zbeta1-zalpha1)*zan)
c      zint = zint +0.0251*sin(zalpha0+zbeta0-(zalpha1+zbeta1)*zan)	
c      pdecli = asin (zint)

c ********************* version before 01/01/2000 *******

      pdecli= asin (sin(pls)*sin(obliquit*pi/180.))

c ********************* version after 01/01/2000 *******
c      i=obliquit*pi/180.
c      pdecli=asin(sin(pls)*sin(i)/sqrt(sin(pls)**2+
c     & cos(pls)**2*cos(i)**2))
c ******************************************************



      RETURN
      END