C=============================================================
C Subroutine to return distance, longitude, and declination of 
C the sun in planetocentric coordinates from orbital elements
C
C INPUTS: 
C edays = time since perihelion (earth days)
C a = semimajor axis (AU)
C ecc = eccentricity
C omega = Ls of perihelion, relative to equinox (radians)
C eps = obliquity (radians)
C
C OUTPUTS:
C Ls = areocentric longitude (radians)
C dec = planetocentric solar declination (radians)
C r = heliocentric distance (AU)
C=============================================================

      SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r)
      implicit none
      real*8 edays,a,ecc,omega,eps  ! input
      real*8 Ls,dec,r  ! output
      real*8 pi,d2r
      parameter (pi=3.1415926535897932,d2r=pi/180.d0)
      integer j
      real*8 M,E,nu,T,Eold

c     T = orbital period (days)
      T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400.

c     M = mean anomaly (radians)
      M = -omega + 2.*pi*edays/T   ! bad for roundoff when edays is large
      !M = 2.*pi*edays/T  ! M=0 at perihelion

c     E = eccentric anomaly 
c     solve M = E - ecc*sin(E) by newton method
      E = M
      do j=1,10
         Eold = E
         E = E - (E - ecc*sin(E) - M)/(1.-ecc*cos(E))
         if (abs(E-Eold)<1.e-8) exit
      enddo

c     nu = true anomaly
      !nu = acos(cos(E)-ecc/(1.-ecc*cos(E)))
      !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E))
      !nu = atan(sqrt(1-ecc^2)*sin(E)/(1-cos(E)))
      nu = 2.*atan(sqrt((1.+ecc)/(1.-ecc))*tan(E/2.))

      !r = a*(1.-ecc**2)/(1.+ecc*cos(nu))
      r = a*(1-ecc*cos(E))
      Ls = mod(nu + omega,2.*pi)   
      dec = asin(sin(eps)*sin(Ls))

      END