riset.f90

Go to the documentation of this file.

!> \file riset.f90  Compute rise, transit and set times, or beginning and end of twilight for libTheSky
 
 
!  Copyright (c) 2002-2025  Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
!   
!  This file is part of the libTheSky package, 
!  see: https://libthesky.sf.net/
!  
!  This is free software: you can redistribute it and/or modify it under the terms of the
!  European Union Public Licence 1.2 (EUPL 1.2).
!  
!  This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
!  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!  See the EU Public Licence for more details.
!  
!  You should have received a copy of the European Union Public Licence along with this code.
!  If not, see <https://www.eupl.eu/1.2/en/>.
 
 
 
!***********************************************************************************************************************************
!> \brief Procedures to determine rise, transit and set times
 
module thesky_riset
  implicit none
  save
  
contains
  
  !*********************************************************************************************************************************
  !> \brief  Rise, transit and set times routine for Sun, Moon, planets and asteroids.
  !!         This version recomputes positions (low accuracy first, then full accuracy) during the convergence process.
  !!         See riset_ipol() for a version using interpolation, to be used for planets only(!)
  !! 
  !! \param jd         Julian day number.  The rise/transit and set data are computed for the (calendar) DAY of this JD (midnight-midnight),
  !!                   unless for_night=.true. in which case this happens for the NIGHT starting at JD.
  !! \param pl         Planet/object number  (-1 - 9: ES, Moon, Mer-Plu; >10000: asteroids)
  !! 
  !! \param rt         Rise time (hours) (output)
  !! \param tt         Transit time (hours) (output)
  !! \param st         Set time (hours) (output)
  !! 
  !! \param rh         Rising wind direction (rad) (output)
  !! \param ta         Transit altitude (rad) (output)
  !! \param sh         Setting wind direction (rad) (output)
  !! 
  !! \param rsAlt      Altitude to return rise/set data for (degrees; 0. is actual rise/set).  rsAlt>90: compute transit only
  !! \param for_night  Compute rise/tranit/set data for the NIGHT starting at the indicated JD, rather than the calendar day.
  !!                   Times between 0h and 12h are for the NEXT calendar day! (optional; default: false)
  !! 
  !! \param ltime      Passed to planet_position(). If .true., include light time, doubling the CPU time while gaining a bit of 
  !!                     accuracy (optional; default: false)
  !! \param cWarn      Warn upon convergence failure (optional; default: true)
  !! 
  !! \param converge   Number of iterations needed to converge (optional) (output)
  !! 
  !! \note
  !! - for rsAlt = 0.d0, rise and set times are computed
  !! - for (rsAlt.ne.0), the routine calculates when alt=rsAlt is reached
  !! - use rsAlt=-6,-12,-18 for the sun for twilight calculations (rsAlt is expressed in degrees).
  !! - Returns times, rise/set azimuth and transit altitude
  !! - Moon transit error ~0.15s (?)
  !! 
  !! \todo
  !! - This version sometimes finds the answer tmRad(i)>2pi, which is the answer for the next day...
  !!   - (only?) solution: use a solver on JD for az,alt
  !! 
  !! \note Speed wrt riset_ipol(), transit only, when using low-accuracy approximations:
  !! - Moon: 270x faster (used to be slowest by factor 3-5 wrt planets, factor 6.5 wrt Sun)
  !! - Sun:  340x faster - full accuracy (VSOP): ~2-3x SLOWER!
  !! - Mer:   20% slower
  !! - Ven:    4% slower
  !! - Mar:   10% slower
  !! - Jup:   25% slower
  !! - Nep:   30% slower
  !! 
  !! \see
  !! - riset_ipol()
  !! - Meeus, Astronomical algorithms, Ch.15, but with geographic longitude east of Greenwich defined as > 0
  
  
  subroutine riset(jd,pl,  rt,tt,st, rh,ta,sh,  rsAlt, for_night, ltime, cWarn, converge)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2, d2r,am2r, r2h, enpname, earthr,au
    use sufr_angles, only: rev, rev2
    use sufr_date_and_time, only: cal2jd,jd2cal
    use sufr_time2string, only: hm
    use sufr_numerics, only: deq0, deq
    
    use thesky_local, only: tz, lat0,lon0,deltat
    use thesky_planets, only: planet_position, planet_position_la
    use thesky_planetdata, only: planpos
    use thesky_coordinates, only: refract
    
    implicit none
    integer, intent(in) :: pl
    real(double), intent(in) :: jd, rsAlt
    real(double), intent(out) :: rt,tt,st,rh,ta,sh
    logical, intent(in), optional :: for_night, ltime, cWarn
    integer, intent(out), optional :: converge(3)
    
    integer :: evi,iter,yr,mnt,tc,evMax, lconverge(3)
    real(double) :: dy,day0,  jd0,jd1,tmRad(3),tmhr(3),  ra,dec
    real(double) :: rsa,cosH0,h0,agst0,th0,dTmRad,accur,  ha,alt,azAlt(3)
    character :: event(3)*(13)
    logical :: use_vsop, lfor_night, lltime, lcWarn
    
    
    ! Take care of optional input variables:
    lfor_night = .false.  ! Compute rise/transit/set times for a given NIGHT, rather than a calendar day
    if(present(for_night)) lfor_night = for_night
    
    lltime = .false.                    ! Call planet_position() IGNORING light time by default (faster, lower accuracy)
    if(present(ltime)) lltime = ltime
    
    lcwarn = .true.                     ! Warn upon convergence failure
    if(present(cwarn)) lcwarn = cwarn
    
    ! Use the old interpolation routine for all but Moon and Sun:
    if(pl.ne.0.and.pl.ne.3) then  ! Not the Moon or the Sun, but a planet, comet or asteroid
       call riset_ipol(jd,pl, rt,tt,st, rh,ta,sh, rsalt, for_night=lfor_night, ltime=lltime, cwarn=lcwarn)
       return
    end if
    
    
    alt=0.d0; ha=0.d0; h0=0.d0; azalt=0.d0; tmrad=0.d0
    tc = 1  ! 0: geocentric (unused), 1: topocentric;  see also different rsa for the Moon below
    event = ['Transit time ','Rise time    ','Set time     ']
    
    if(abs(rsalt).gt.1.d-9) then
       rsa = rsalt*d2r  ! Use a user-specified altitude
    else
       rsa = -0.5667d0*d2r                      ! Standard rise/set altitude for planets
       if(pl.eq.3) rsa = rsa - 16.d0*am2r       ! Correct for radius of Sun, -16 arcminutes  ! rsa = -0.8333d0*d2r - default for Sun
       if(pl.eq.0) rsa = 0.125d0*d2r            ! Approximate standard rise/set altitude for Moon, including parallax
    end if
    
    evmax = 3                     ! 'Maximum' event to compute: transit only: evMax=1, +rise/set: evMax=3
    if(rsalt.gt.90.d0) evmax = 1  ! Compute transit time and altitude only
    
    call jd2cal(jd, yr,mnt,dy)
    day0 = dble(floor(dy))-tz/24.d0     ! Midnight local time at the start of the desired DAY, needed for agst0
    if(lfor_night) day0 = day0 + 0.5d0  ! Noon local time at the start of the desired NIGHT, needed for agst0
    jd0  = cal2jd(yr,mnt,day0)          ! Midnight local time at the start of the desired day, needed for agst0
    
    call planet_position_la(jd0, pl, 3, 60)  ! Compute low-accuracy positions - calc=2 computes ra,dec, calc=3 computes AGST
    
    ra    = planpos(5+tc*20)
    dec   = planpos(6+tc*20)
    agst0 = planpos(45)       ! AGST for midnight
    
    if(pl.eq.0) then  ! Moon
       if(tc.eq.0) then  ! Geocentric (unused)
          rsa = asin(earthr/(planpos(4)*au))*0.7275d0-0.5667d0*d2r  ! Exact altitude for Moon, ~0.1-0.2deg, for geocentric coord.
       else  ! tc.eq.1 - topocentric:
          rsa = -0.8333d0*d2r  ! For Moon, in combination with topocentric coordinates
       end if
    end if
    
    
    if(evmax.eq.3) then  ! Computing transit, rise and set
       cosh0 = (sin(rsa)-sin(lat0)*sin(dec)) / (cos(lat0)*cos(dec))  ! Cosine of the hour angle of rise/set; Meeus, Eq.15.1
       if(abs(cosh0).gt.1.d0) then  ! Body never rises/sets
          evmax = 1                 ! Compute transit time and altitude only
       else
          h0 = rev(2*acos(cosh0))/2.d0  ! Hour angle of rise/set
       end if
    end if
    
    
    tmrad(1) = rev(ra - lon0 - agst0)  ! Transit time in radians; tmRad(1)=m0, tmRad(2)=m1, tmRad(3)=m2 in Meeus, but lon0 > 0 for E
    if(evmax.eq.3) then
       tmrad(2) = rev(tmrad(1) - h0)   ! Rise time in radians
       tmrad(3) = rev(tmrad(1) + h0)   ! Set time in radians
    end if
    
    
    do evi=1,evmax         ! Transit, rise, set
       iter = 0
       accur = 1.d-3       ! Accuracy (time in rad).  Initially 1d-3 (~14s), later 1d-5 (~0.14s). Don't make this smaller than 1d-16
       use_vsop = .false.  ! Initially
       
       dtmrad = huge(dtmrad)
       do while(abs(dtmrad).ge.accur .or. .not.use_vsop)
          th0 = agst0 + 1.002737909350795d0*tmrad(evi)   ! Solar day in sidereal days in 2000; Expl.Suppl.tt.Astr.Almanac 3rdEd 
          jd1 = jd0 + tmrad(evi)/pi2 + deltat/86400.d0   ! ra, s -> days                       Eq.3.17 (removed '...37...' typo)
          
          if(abs(dtmrad).le.accur) then
             use_vsop = .true.
             accur = 1.d-5  ! 1d-5 rad ~ 0.14s.  Changing this to 1.d-4 (~1s) speeds the code yearly_moon_table code up by ~30%.  Don't make this smaller than 1d-16.
          end if
          
          if(use_vsop) then
             ! Accuracy of 1 min = 0.25 deg = 4e-3 rad.  1.d-6,1.d-2 saves ~43% CPU time for the Sun, <1 round-off errors/year
             call planet_position(jd1,pl, lbaccur=1.d-6,raccur=1.d-2, ltime=lltime)  ! Uses truncated VSOP87
          else
             call planet_position_la(jd1,pl, 2,60)  ! Computes low-accuracy positions - calc=2 computes ra,dec - use 60 terms
          end if
          
          ra  = planpos(5+tc*20)  ! Right ascension
          dec = planpos(6+tc*20)  ! Declination
          
          ha  = rev2(th0 + lon0 - ra)                                  ! Hour angle;  Meeus p.103
          ! ha  = planpos(8+tc*20)                                       ! Hour angle -/- DeltaT
          alt = asin(sin(lat0)*sin(dec) + cos(lat0)*cos(dec)*cos(ha))  ! Altitude;  Meeus, Eq.13.6
          
          ! Correction to transit/rise/set times:
          if(evi.eq.1) then  ! Transit
             dtmrad = -ha
          else               ! Rise/set
             dtmrad = (alt-rsa)/(cos(dec)*cos(lat0)*sin(ha))
          end if
          tmrad(evi) = tmrad(evi) + dtmrad
          
          iter = iter+1
          if(iter.gt.30) exit  ! do-while loop
       end do  ! do while(abs(dTmRad).ge.accur .or. .not.use_vsop)
       
       
       if(iter.gt.30) then  ! Convergence failed
          if(lcwarn .and. (pl.ne.3.or.nint(rsalt).ne.-18)) &
               write(0,'(A,F10.3,A)') '  * WARNING:  riset():  Riset failed to '// &
               'converge: '//trim(enpname(pl))//'  '//trim(event(evi)),rsalt,'d *'
          
          tmrad(evi) = 0.d0
          azalt(evi) = 0.d0
       else                    ! Result converged, store it
          if(evi.eq.1) then
             azalt(evi) = alt  ! Transit altitude
          else
             azalt(evi) = atan2(sin(ha),(cos(ha)*sin(lat0)-tan(dec)*cos(lat0)))  ! Rise,set hour angle -> azimuth
          end if
       end if
       
       ! The Moon may rise/transit/set AFTER MIDNIGHT (t>24h; tmRad > 2pi):
       if(pl.eq.0 .and. tmrad(evi).gt.pi2) then  ! Moon;  in case after m_i = m_i+1, m_f > 1
          tmrad(evi) = 0.d0
          azalt(evi) = 0.d0
       end if
       
       ! - If the time is < 0, the Moon's SET times (only?) may repeat the previous value on days where
       !   the Moon does not set.
       ! - Don't do this for planets, which may rise/transit/set twice a day!
       ! - Don't do this for twilight (Sun with rsAlt<0), which can cross midnight in summer!
       ! - This may be needed for planets/the Sun near the poles...
       ! 
       ! if(tmRad(evi).lt.0.d0 .and. deq0(rsAlt)) then  ! Pre-2025
       ! if((pl.eq.0.or.pl.eq.3) .and. tmRad(evi).lt.0.d0 .and. deq0(rsAlt)) then  ! Sun and Moon only
       ! if(pl.eq.0 .and. tmRad(evi).lt.0.d0 .and. deq0(rsAlt)) then  ! Moon only
       if(pl.eq.0 .and. tmrad(evi).lt.0.d0) then  ! Moon only - all altitudes
          tmrad(evi) = 0.d0
          azalt(evi) = 0.d0
       end if
       
       lconverge(evi) = iter  ! Number of iterations needed to converge
       
    end do  ! evi: 1-evMax to cycle through transit, rise, set
    
    
    ! Store results:
    tmhr = tmrad * r2h                ! Times radians -> hours
    if(lfor_night) then
       tmhr = tmhr + 12   ! Times relative to noon -> relative to midnight
       where(deq(tmhr,12.d0)) tmhr = 0.d0
    end if
    
    tt = tmhr(1)                       ! Transit time
    rt = tmhr(2)                       ! Rise time
    st = tmhr(3)                       ! Set time
    
    ta = azalt(1) + refract(azalt(1))  ! Transit altitude + refraction
    rh = azalt(2)                      ! Rise azimuth
    sh = azalt(3)                      ! Set azimuth
    
    
    ! Take care of optional output variables:
    if(present(converge)) converge = lconverge  ! Number of iterations needed to converge
    
  subroutine riset(jd,pl,  rt,tt,st, rh,ta,sh,  rsAlt, for_night, ltime, cWarn, converge) …
  end subroutine riset
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Old routine for rise, transit and set times for planets and asteroids - don't use this for Sun and Moon
  !!         Computes three sets of planet positions, and interpolates between them.
  !! 
  !! \param jd         Julian day number.  The rise/transit and set data are computed for the (calendar) DAY of this JD (midnight-midnight),
  !!                   unless for_night=.true. in which case this happens for the NIGHT starting at JD.
  !! \param pl         Planet/object number  (-1 - 9: ES, Moon, Mer-Plu; >10000: asteroids)
  !!
  !! \param rt         Rise time (hours) (output)
  !! \param tt         Transit time (hours) (output)
  !! \param st         Set time (hours) (output)
  !! 
  !! \param rh         Rising wind direction (rad) (output)
  !! \param ta         Transit altitude (rad) (output)
  !! \param sh         Setting wind direction (rad) (output)
  !! 
  !! \param rsAlt      Altitude to return rise/set data for (degrees; 0. is actual rise/set).  rsAlt>90: compute transit only
  !! \param for_night  Compute rise/tranit/set data for the NIGHT starting at the indicated JD, rather than the calendar day.
  !!                   Times between 0h and 12h are for the NEXT calendar day! (optional; default: false)
  !! 
  !! \param ltime      Passed to planet_position(). If .true., include light time, doubling the CPU time while gaining a bit of accur.
  !! \param cWarn      Warn upon convergence failure (optional; default: true)
  !! 
  !! \note
  !! - for rsAlt = 0.d0, rise and set times are computed
  !! - for (rsAlt.ne.0), the routine calculates when alt=rsAlt is reached
  !! - use rsAlt=-6,-12,-18 for the sun for twilight calculations (rsAlt is expressed in degrees).
  !! - transit error for Moon <= ~10s (due to interpolation); typically ~4s?
  !! - code gets confused around midnight - it sometimes returns a small negative tmdy(i), which is on the day before(!)
  !!
  !! \see
  !! - Meeus, Astronomical algorithms, Ch.15, but with geographic longitude east of Greenwich defined as > 0
  
  subroutine riset_ipol(jd,pl, rt,tt,st, rh,ta,sh, rsAlt, for_night, ltime, cWarn)
    use sufr_kinds, only: double
    use sufr_constants, only: pi,pi2, d2r,am2r, enpname, earthr,au
    use sufr_system, only: warn
    use sufr_angles, only: rev, rev2
    use sufr_date_and_time, only: cal2jd,jd2cal
    use sufr_numerics, only: deq0, deq
    
    use thesky_local, only: tz, lat0,lon0,deltat
    use thesky_planets, only: planet_position
    use thesky_planetdata, only: planpos
    use thesky_coordinates, only: refract
    
    implicit none
    integer, intent(in) :: pl
    real(double), intent(in) :: jd, rsAlt
    real(double), intent(out) :: rt,tt,st, rh,ta,sh
    logical, intent(in), optional :: for_night, ltime, cWarn
    
    integer :: evi,iter,yr,mnt,tc,evMax, indic
    real(double) :: dy,day0,  jd0,jd1,jd2, tmdy(3),dTmdy(3), tmhr(3)
    real(double) :: ra0,dec0,ra1,dec1,ra2,dec2,ra,dec, rsa,cosH0,h0,agst0,th0,n,dtm,accur,  ha,alt,azAlt(3)
    character :: event(3)*(13)
    logical :: lfor_night, lltime, lcWarn
    save :: indic
    
    
    ! Take care of optional input variables:
    lfor_night = .true.
    if(present(for_night)) lfor_night = for_night
    
    lltime = .true.
    if(present(ltime)) lltime = ltime
    
    lcwarn = .true.
    if(present(cwarn)) lcwarn = cwarn
    
    if((pl.eq.0.or.pl.eq.3) .and. indic.ne.12345) then
       call warn("riset_ipol():  Don't use this routine for Sun or Moon - use riset() instead.", 0)
       indic = 12345
    end if
    
    rt=0.d0; tt=0.d0; st=0.d0;  rh=0.d0; ta=0.d0; sh=0.d0
    alt = 0.d0;  ha = 0.d0;  dtmdy = 0.d0;  h0 = 0.d0;  dec = 0.d0
    tc = 1        ! 0: geocentric (unused), 1: topocentricl;  see also different rsa for the Moon below
    event = ['Transit time ','Rise time    ','Set time     ']
    accur = 1.d-6          ! Accuracy.  1d-6 ~ 0.1s. Don't make this smaller than 1d-16
    
    
    rsa = -0.5667d0*d2r                     ! Standard rise/set altitude for planets
    ! NOTE: Sun and Moon should not use this routine(!)
    if(pl.eq.3) rsa = rsa - 16.0d0*am2r     ! Compensate for radius of Sun, -16 arcminutes  ! rsa = -0.8333d0*d2r - default for Sun
    if(pl.eq.0) rsa = 0.125d0*d2r           ! Approximate standard altitude for Moon
    
    if(abs(rsalt).gt.1.d-9) rsa = rsalt*d2r
    
    call jd2cal(jd,yr,mnt,dy)
    day0 = dble(floor(dy))-tz/24.d0  ! Midnight local time of the desired day, needed for agst0
    if(lfor_night) day0 = day0 + 0.5d0  ! Noon local time at the start of the desired NIGHT, needed for agst0
    jd0  = cal2jd(yr,mnt,day0-1.d0)  ! Midnight local time, day before
    jd1  = cal2jd(yr,mnt,day0)       ! Midnight local time, needed for agst0
    jd2  = cal2jd(yr,mnt,day0+1.d0)  ! Midnight local time, day after
    
    call planet_position(jd0,pl, lbaccur=1.d-6,raccur=1.d-2, ltime=lltime)
    ra0  = planpos(5+tc*20)
    dec0 = planpos(6+tc*20)
    
    call planet_position(jd1,pl, lbaccur=1.d-6,raccur=1.d-2, ltime=lltime)
    ra1   = planpos(5+tc*20)
    dec1  = planpos(6+tc*20)
    agst0 = planpos(45)       ! AGST for midnight
    
    ! NOTE: Sun and Moon should not use this routine(!)
    if(pl.eq.0) then
       if(tc.eq.0) then  ! Geocentric (unused)
          rsa = asin(earthr/(planpos(4)*au))*0.7275d0-0.5667d0*d2r  ! Exact altitude for Moon, ~0.1-0.2deg, for geocentric coord.
       else  ! tc.eq.1 - topocentric:
          rsa = -0.8333d0*d2r  ! For Moon, in combination with topocentric coordinates
       end if
    end if
    
    call planet_position(jd2,pl, lbaccur=1.d-6,raccur=1.d-2, ltime=lltime)
    ra2  = planpos(5+tc*20)
    dec2 = planpos(6+tc*20)
    
    if(abs(ra1-ra0).gt.2..or.abs(ra2-ra1).gt.2.) then  ! ra flips 2pi <-> 0
       if(ra0.ge.4.5.and.ra1.ge.4.5.and.ra2.le.2.) ra2 = ra2+2*pi
       if(ra0.ge.4.5.and.ra1.le.2.and.ra2.le.2.)   ra0 = ra0-2*pi
       if(ra0.le.2.and.ra1.ge.4.5.and.ra2.ge.4.5)  ra0 = ra0+2*pi
       if(ra0.le.4.5.and.ra1.le.2.and.ra2.ge.4.5)  ra2 = ra2-2*pi
       ! write(6,'(A)') 'RA flipped'
    end if
    
    cosh0 = (sin(rsa)-sin(lat0)*sin(dec2)) / (cos(lat0)*cos(dec2))  ! Cosine of the hour angle of rise/set; Meeus, Eq.15.1
    
    tmdy = 0.d0
    azalt = 0.d0
    
    evmax = 3                    ! 'Maximum' event to compute: transit only: evMax=1, +rise/set: evMax=3
    if(abs(cosh0).gt.1.d0) then  ! Body never rises/sets
       evmax = 1                 ! Compute transit time and altitude only
    else
       h0 = rev(2*acos(cosh0))/2.d0  ! Hour angle of rise/set
    end if
    
    
    tmdy(1) = rev(ra2 - lon0 - agst0)/pi2 + dtmdy(1)    ! Transit time in days;  Meeus Eq. 15.2a, but lon0 > 0 for E
    if(evmax.eq.3) then
       tmdy(2) = rev(tmdy(1)*pi2 - h0)/pi2 + dtmdy(2)   ! Rise time in days;  Meeus Eq.15.2b
       tmdy(3) = rev(tmdy(1)*pi2 + h0)/pi2 + dtmdy(3)   ! Set time in days;  Meeus Eq.15.2c
    end if
    
    evi = 1  ! i=1,evMax
    do while(evi.le.evmax)
       iter = 0
       
       dtm = huge(dtm)
       do while(abs(dtm).ge.accur)
          th0 = agst0 + 6.300388092591991d0*tmdy(evi)  ! Siderial day in radians; Meeus, p.103
          n   = tmdy(evi) + deltat/86400.d0
          ra  = rsipol( ra0, ra1, ra2, n)  ! Interpolate right ascension
          dec = rsipol(dec0,dec1,dec2, n)  ! Interpolate declination
          
          ha = rev2(th0 + lon0 - ra)                                   ! Hour angle;  Meeus p.103
          alt = asin(sin(lat0)*sin(dec) + cos(lat0)*cos(dec)*cos(ha))  ! Altitude; Meeus, Eq.13.6
          
          ! Correction to transit/rise/set times:
          if(evi.eq.1) then  ! Transit
             dtm = -ha/pi2
          else              ! Rise/set
             dtm = (alt-rsa)/(pi2*cos(dec)*cos(lat0)*sin(ha))
          end if
          tmdy(evi) = tmdy(evi) + dtm
          
          iter = iter+1
          if(iter.gt.30) exit  ! do-while loop
       end do  ! do while(abs(dtm).ge.accur)
       
       
       if(iter.gt.30) then  ! Convergence failed
          if(lcwarn .and. (pl.ne.3.or.nint(rsalt).ne.-18)) &
               write(0,'(A,F10.3,A)') '  * WARNING:  riset_ipol():  Riset failed '// &
               'to converge: '//trim(enpname(min(pl,19)))//'  '//trim(event(evi)),rsalt,'d'
          
          tmdy(evi) = 0.d0
          azalt(evi) = 0.d0
       else               ! Result converged, store it
          if(evi.eq.1) then
             azalt(evi) = alt  ! Transit altitude
          else
             azalt(evi) = atan2(sin(ha),(cos(ha)*sin(lat0)-tan(dec)*cos(lat0)))  ! Rise,set hour angle -> azimuth
          end if
       end if
       
       
       ! NOTE: Sun and Moon should not use this routine(!)
       if(pl.eq.0) then
          if(tmdy(evi).lt.0.) then    ! Sometimes, (at least) the Moon doesn't r,t,s twice, because m_i slightly > 1, and rev(m_i) << 1,
             dtmdy(evi) = dtmdy(evi) + 1.d0     ! then initially m_f < 0, though if m_i = m_i + 1,  then 0 < m_f < 1
             
             ! Restart the do loop:
             evi = 1
             tmdy(1) = rev(ra2 - lon0 - agst0)/pi2 + dtmdy(1)  ! Transit; tmdy(1)=m0, tmdy(2)=m1, tmdy(3)=m2 in Meeus, but lon0 > 0 for E
             if(evmax.eq.3) then
                tmdy(2) = rev(tmdy(1)*pi2 - h0)/pi2 + dtmdy(2)   ! Rise
                tmdy(3) = rev(tmdy(1)*pi2 + h0)/pi2 + dtmdy(3)   ! Set
             end if
             cycle  ! i
          end if  ! if(tmdy(evi).lt.0.)
          
          if(tmdy(evi).gt.1.d0) then  ! in case after m_i = m_i+1, m_f > 1
             tmdy(evi) = 0.d0
             azalt(evi) = 0.d0
          end if  ! if(tmdy(evi).gt.1.d0)
       end if  ! if(pl.eq.0)
       
       
       ! NOTE: Sun and Moon should not use this routine(!)
       ! - If the time is < 0, the Moon's SET times (only?) may repeat the previous value on days where
       !   the Moon does not set.
       ! - Don't do this for planets, which may rise/transit/set twice a day!
       ! - Don't do this for twilight (Sun with rsAlt<0), which can cross midnight in summer!
       ! - This may be needed for planets/the Sun near the poles...
       ! 
       ! if(tmdy(evi).lt.0.d0 .and. deq0(rsAlt)) then  ! Pre-2025
       ! if((pl.eq.0.or.pl.eq.3) .and. tmdy(evi).lt.0.d0 .and. deq0(rsAlt)) then
       ! if(pl.eq.0 .and. tmdy(evi).lt.0.d0 .and. deq0(rsAlt)) then  ! Moon only
       if(pl.eq.0 .and. tmdy(evi).lt.0.d0) then  ! Moon only - all altitudes
          tmdy(evi) = 0.d0
          azalt(evi) = 0.d0
       end if
       
       evi = evi+1  ! evi=1,evMax: cycle through transit, rise, set
    end do  ! do while(evi.le.evMax)
    
    
    ! Store results:
    tmhr = tmdy * 24                   ! Times days -> hours
    if(lfor_night) then
       tmhr = tmhr + 12   ! Times relative to noon -> relative to midnight
       where(deq(tmhr,12.d0)) tmhr = 0.d0
    end if
    
    rt = tmhr(2)                       ! Rise time
    tt = tmhr(1)                       ! Transit time
    st = tmhr(3)                       ! Set time
    
    rh = azalt(2)                      ! Rise azimuth
    ta = azalt(1) + refract(azalt(1))  ! Transit altitude + refraction
    sh = azalt(3)                      ! Set azimuth
    
  subroutine riset_ipol(jd,pl, rt,tt,st, rh,ta,sh, rsAlt, for_night, ltime, cWarn) …
  end subroutine riset_ipol
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief Rise, transit and set times routine for an object with fixed ra & dec
  !!
  !! \param jd         Julian day number
  !! \param ra         Right ascension (rad)
  !! \param dec        Declination (rad)
  !!
  !! \param rt         Rise time (hours) (output)
  !! \param tt         Transit time (hours) (output)
  !! \param st         Set time (hours) (output)
  !! 
  !! \param rh         Rising wind direction (rad) (output)
  !! \param ta         Transit altitude (rad) (output)
  !! \param sh         Setting wind direction (rad) (output)
  !!
  !! \param rsAlt      Altitude to return rise/set data for (degrees; 0. is actual rise/set).  rsAlt>90: compute transit only
  !! \param for_night  Compute rise/tranit/set data for the NIGHT starting at the indicated JD, rather than the calendar day.
  !!                     Times between 0h and 12h are for the NEXT calendar day! (optional; default: false)
  !!
  !! \param cWarn      Warn upon convergence failure (optional; default: true)
  !!
  !! for rsAlt = 0., rise and set times are computed
  !! for rsAlt.ne.0, the routine calculates when alt=rsAlt is reached
  !! use rsAlt=-6,-12,-18 for the sun for twilight calculations
  !! (rsAlt is expressed in degrees).
  !!
  !! \see
  !! - Meeus, Astronomical algorithms, Ch.15, but with geographic longitude east of Greenwich defined as > 0
  
  subroutine riset_ad(jd, ra,dec, rt,tt,st,rh,ta,sh, rsAlt, for_night, cWarn)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2, d2r
    use sufr_angles, only: rev, rev2
    use sufr_date_and_time, only: cal2jd, jd2cal
    use sufr_numerics, only: deq0, deq
    
    use thesky_planets, only: planet_position
    use thesky_planetdata, only: planpos
    use thesky_local, only: tz, lat0,lon0
    
    implicit none
    real(double), intent(in) :: jd, ra,dec, rsAlt
    real(double), intent(out) :: rt,tt,st,rh,ta,sh
    logical, intent(in), optional :: for_night, cWarn
    
    integer :: evi,iter,yr,mnt,evMax
    real(double) :: dy,day0,jd0, tmdy(3),tmhr(3),  rsa,cosH0,h0,agst0,th0,dtm,accur,  ha,alt,azAlt(3)  !,n
    character :: event(3)*(13)
    logical :: lfor_night, lcWarn
    
    ! Take care of optional input variables:
    lfor_night = .true.
    if(present(for_night)) lfor_night = for_night
    
    lcwarn = .true.
    if(present(cwarn)) lcwarn = cwarn
    
    alt = 0.d0;  ha = 0.d0;  h0 = 0.d0
    accur = 1.d-6  ! Accuracy.  1d-6 ~ 0.1s. Don't make this smaller than 1d-16
    event = ['Transit time ','Rise time    ','Set time     ']
    
    rsa = rsalt*d2r
    evmax = 3                     ! 'Maximum' event to compute: transit only: evMax=1, +rise/set: evMax=3
    if(rsalt.gt.90.d0) evmax = 1  ! Compute transit time and altitude only
    
    call jd2cal(jd, yr,mnt,dy)
    day0 = dble(floor(dy))-tz/24.d0     ! Midnight local time of the desired day, needed for agst0
    if(lfor_night) day0 = day0 + 0.5d0  ! Noon local time at the start of the desired NIGHT, needed for agst0
    jd0  = cal2jd(yr,mnt,day0)
    
    call planet_position(jd0,3)
    agst0 = planpos(45)            ! AGST for midnight (or noon)
    
    tmdy = 0.d0
    azalt = 0.d0
    if(evmax.eq.3) then
       cosh0 = (sin(rsa)-sin(lat0)*sin(dec))/(cos(lat0)*cos(dec))  ! Cosine of the hour angle of rise/set; Meeus, Eq.15.1
       if(abs(cosh0).gt.1.d0) then  ! Body never rises/sets
          evmax = 1                 ! Compute transit time and altitude only
       else
          h0 = rev(2*acos(cosh0))/2.d0  ! Hour angle of rise/set
       end if
    end if
    
    tmdy(1) = rev(ra - lon0 - agst0)/pi2    ! Transit;  tmdy(1)=m0, tmdy(2)=m1, tmdy(3)=m2 in Meeus, but lon0 > 0 for E
    if(evmax.eq.3) then
       tmdy(2) = rev(tmdy(1)*pi2 - h0)/pi2  ! Rise
       tmdy(3) = rev(tmdy(1)*pi2 + h0)/pi2  ! Set
    end if
    
    
    do evi=1,evmax
       iter = 0
       dtm = huge(dtm)
       do while(abs(dtm).ge.accur)
          th0 = agst0 + 6.300388092591991d0*tmdy(evi)  ! Siderial day in radians; Meeus, p.103
          
          ha = rev2(th0 + lon0 - ra)                                   ! Hour angle;  Meeus p.103
          alt = asin(sin(lat0)*sin(dec) + cos(lat0)*cos(dec)*cos(ha))  ! Altitude;  Meeus Eq. 13.6
          
          ! Correction on transit/rise/set times:
          if(evi.eq.1) then  ! Transit
             dtm = -ha/pi2
          else              ! Rise/set
             dtm = (alt-rsa)/(pi2*cos(dec)*cos(lat0)*sin(ha))
          end if
          tmdy(evi) = tmdy(evi) + dtm
          
          iter = iter+1
          if(iter.gt.30) exit  ! do-while loop
       end do  ! do while(abs(dtm).ge.accur)
       
       
       if(iter.gt.30) then  ! Convergence failed
          if(lcwarn) write(0,'(A,F10.3,A)') '  * WARNING:  riset_ad():  Riset failed to converge: '//trim(event(evi)),rsalt,'d *'
          
          tmdy(evi) = 0.d0
          azalt(evi) = 0.d0
       else               ! Result converged, store it
          if(evi.eq.1) then
             azalt(evi) = alt  ! Transit altitude
          else
             azalt(evi) = atan2(sin(ha),(cos(ha)*sin(lat0)-tan(dec)*cos(lat0)))  ! Rise,set hour angle -> azimuth
          end if
       end if
       
       if(tmdy(evi).lt.0. .and. deq0(rsalt)) then
          tmdy(evi) = 0.d0
          azalt(evi) = 0.d0
       end if
       
    end do  ! evi
    
    
    ! Store results:
    tmhr = tmdy * 24                   ! Times days -> hours
    if(lfor_night) then
       tmhr = tmhr + 12   ! Times relative to noon -> relative to midnight
       where(deq(tmhr,12.d0)) tmhr = 0.d0
    end if
    
    rt = tmhr(2)   ! Rise time
    tt = tmhr(1)   ! Transit time
    st = tmhr(3)   ! Set time
    
    rh = azalt(2)  ! Rise azimuth
    ta = azalt(1)  ! Transit altitude
    sh = azalt(3)  ! Set azimuth
    
  subroutine riset_ad(jd, ra,dec, rt,tt,st,rh,ta,sh, rsAlt, for_night, cWarn) …
  end subroutine riset_ad
  !*********************************************************************************************************************************
  
  
  
  !*********************************************************************************************************************************
  !> \brief Compute the date (m,d) at which an object with ra can be observed best, i.e., it transits at midnight  
  !!
  !! \param  jd00  Julian day?
  !! \param  ra0   Right ascension
  !!
  !! \param m     Month of best observation (output)
  !! \param d     Day of month of best observation (output)
  !!
  !! \note 
  !! - VERY SLOW!!! -> use best_obs_date_ra instead!!!  (but not (yet) in libTheSky...)
  !! - riset() should be faster now (04/2010)
  
  subroutine best_obs_date(jd00,ra0, m,d)
    use sufr_kinds, only: double
    use sufr_constants, only: r2h
    use sufr_angles, only: rev,rv12
    use sufr_date_and_time, only: jd2cal
    use sufr_numerics, only: deq
    
    use thesky_planets, only: planet_position
    use thesky_planetdata, only: planpos, nplanpos
    use thesky_local, only: tz
    use thesky_datetime, only: gettz
    
    implicit none
    real(double), intent(in) :: jd00,ra0
    integer, intent(out) :: m,d
    real(double) :: jd0,oldjd,z(nplanpos)
    real(double) :: rt,stt,ott,st,rh,ta,sh,midnight,hh,dt,dd
    integer :: y,iter
    
    oldjd = -99.d0
    jd0 = jd00
    
    dd = 9999.d0
    iter = 0
    do while(abs(dd).gt.0.5d0)
       call riset(jd0,3, rt,stt,st,rh,ta,sh, 0.d0)  ! Sun transit time
       midnight = mod(stt+36.d0,24.d0)
       
       call planet_position(jd0,3)                  ! Position of the Sun
       z = planpos
       
       hh  = rev(z(44) - ra0)                       ! Hour angle is LST - RA
       ott = -hh*r2h                                ! Time of object transit = -ha
       dt  = rv12(ott-midnight)                     ! Difference between transit time and midnight on this day
       dd  = dt/24.d0*365.25d0                      ! Estimate for difference in days between today and day we're looking for
       jd0 = dble(nint(jd0 + dd + 0.5d0)) - 0.5d0   ! Force midnight UT
       tz  = gettz(jd0)
       jd0 = jd0 - tz/24.d0                         ! Local midnight
       
       ! Consider result 'converged':
       if(deq(jd0,oldjd)) exit
       if(iter.gt.5 .and. abs(dd).lt.0.7d0) exit
       if(iter.gt.10 .and. abs(dd).lt.1.d0) exit
       if(iter.gt.100) exit
       
       oldjd = jd0
       iter = iter + 1
    end do
    
    call jd2cal(jd0, y,m,dd)
    d = floor(dd)
    
  subroutine best_obs_date(jd00,ra0, m,d) …
  end subroutine best_obs_date
  !*********************************************************************************************************************************
  
  
  
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Quadratic interpolation of rise/set times for riset_ipol().
  !!
  !! \param  y1  Coordinate for day 1
  !! \param  y2  Coordinate for day 2
  !! \param  y3  Coordinate for day 3
  !! \param   n  Fraction of day
  !! 
  !! \retval rsIpol  Interpolated coordinate
  
  function rsipol(y1,y2,y3, n)
    use sufr_kinds, only: double
    
    implicit none
    real(double), intent(in) :: y1,y2,y3,n
    real(double) :: rsipol, a,b,c
    
    a = y2-y1
    b = y3-y2
    c = b-a
    
    rsipol = y2 + n/2.d0 * (a + b + c*n)   ! 1 + (1+y1/2)*n + (y1/2 + y2 + y3/2)*n^2
    
  function rsipol(y1,y2,y3, n) …
  end function rsipol
  !*********************************************************************************************************************************
  
  
end module thesky_riset
!***********************************************************************************************************************************