MODULE solarlong_mod

      IMPLICIT NONE

      CONTAINS
      
      SUBROUTINE solarlong(pday,psollong)
      use planete_h, only: lsperi, peri_day, year_day, e_elips
      use comcstfi_h, only: pi
      IMPLICIT NONE

c=======================================================================
c   Compute solar longitude psollong (in radians) for a given Mars date
c   pday (where pday==0 at northern Spring Equinox and pday in sols 
c   and fractions thereof)
c=======================================================================

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

      REAL,INTENT(IN) :: pday
      REAL,INTENT(OUT) :: psollong
      
c Local variables:
c ----------------

      REAL :: xref ! mean anomaly
      REAL :: zx0 ! eccentric anomaly
      REAL :: zteta ! true anomaly
      REAL :: zanom,zdx,zz
      INTEGER :: iter


c  compute mean anomaly

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

c  solve equation  zx0 - e * sin (zx0) = xref
c  using Newton method

      zx0=xref+e_elips*sin(xref)
      DO iter=1,10 ! 10 is overkill, typically converges in 1-3 iterations
         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
         if(abs(zdx).le.(1.e-7)) exit
         zx0=zx0+zdx
      ENDDO

      zx0=zx0+zdx
      if(zanom.lt.0.) zx0=-zx0

c compute true anomaly zteta, now that eccentric anomaly zx0 is known

      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))

c compute solar longitude from zteta
      psollong=zteta+lsperi

c handle limit cases where psollong lands outside [0:2*pi]
      IF(psollong.LT.0.) psollong=psollong+2.*pi
      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi

      END SUBROUTINE solarlong

      END MODULE solarlong_mod