vsop.f90

Go to the documentation of this file.

!> \file vsop.f90  Core VSOP routines 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 for VSOP
!!
!! \note
!! VSOP87 accuracies 1900-2100 in milliarcseconds (VSOP87 paper):
!! - Me: 1,  Ve: 6,  Ea: 5,  Ma: 23,  Ju: 20,  Sa: 100,  Ur: 16,  Ne: 30
!!
!! \note
!! Accuracy < 1" in JD2000 +- range in kyr (Wikipedia):
!! - Me-Ma: 4, Ju,Sa: 2, Ur,Ne: 6  (i.e., Earth position more accurate than 1" between -2000 and +6000)
!!
!! \see Bretagnon, P.; Francou, G. Planetary theories in rectangular and spherical variables - VSOP 87 solutions 
!!      http://esoads.eso.org/abs/1988A%26A...202..309B
 
module thesky_vsop
  implicit none
  save
  
  
contains
  
  
  !*********************************************************************************************************************************
  !> \brief Calculate true heliocentric ecliptic coordinates l,b,r for planet pl and the mean ecliptic and
  !!        equinox of date using VSOP87D
  !! 
  !! \param  tm   Dynamical time in Julian Millennia after J2000.0, (tau in Meeus)
  !! \param  pl   Planet to compute position for
  !! 
  !! \param lon  Heliocentric longitude (rad) (output)
  !! \param lat  Heliocentric latitude (rad) (output)
  !! \param rad  Heliocentric distance (AU) (output)
  !! 
  !! \param  LBaccur  Desired accuracy of L,B (rad, optional)
  !! \param  Raccur   Desired accuracy of R (AU, optional)
  !!
  !! \see http://esoads.eso.org/abs/1988A%26A...202..309B
  
  subroutine vsop87d_lbr(tm,pl, lon,lat,rad, LBaccur,Raccur)
    use sufr_kinds, only: double
    use sufr_angles, only: rev
    use thesky_planetdata, only: vsopnls, vsopdat, vsopnblk
    
    implicit none
    real(double), intent(in) :: tm
    integer, intent(in) :: pl
    real(double), intent(out) :: lon,lat,rad
    real(double), intent(in), optional :: LBaccur,Raccur
    
    integer :: li, pow, Nli,Nle, var, nTerm, skip
    real(double) :: fac, lbr(3), accur, desired_accuracy(3)
    
    ! desired_accuracy = 1.d-9  ! 5e-9 rad = 1 mas = VSOP87 accuracy for Mercury in 1900-2100
    desired_accuracy = 0.d0  ! Use all available terms - setting LBR accuracy equal to VSOP87 truncation still skips some terms!
    if(present(lbaccur)) desired_accuracy(1:2) = [lbaccur,lbaccur]
    if(present(raccur))  desired_accuracy(3)   = raccur
    
    lbr = 0.d0
    !$omp parallel do private(Nli,Nle,skip,li,pow,fac,nTerm)  ! Loop body seems to small for any benefit
    do var=1,3  ! L,B,R
       nli = 0
       if(var.ge.2) nli = vsopnls(1,pl)
       if(var.eq.3) nli = nli + vsopnls(2,pl)
       nle = nli + vsopnls(var,pl)
       
       skip = -1  ! Switch to skip terms if sufficient accuracy is reached
       do li=nli+1,nle
          pow = nint(vsopdat(1,li,pl))
          
          if(pow.ne.skip) then
             fac = tm**pow * vsopdat(2,li,pl)
             lbr(var) = lbr(var) + fac * cos( vsopdat(3,li,pl) + vsopdat(4,li,pl)*tm )
             
             ! Determine current accuracy (A*T^n):
             if((desired_accuracy(var).gt.1.d-15) .and. (mod(li,3).eq.0)) then  ! Reduce overhead, experimentally optimised to 3
                nterm = li - vsopnblk(pow,1,pl)+1
                accur = 2*sqrt(dble(nterm)) * abs(fac)
                if(accur.lt.desired_accuracy(var)) skip = pow  ! Skip the remaining terms for this power
             end if
          end if
          
       end do  ! li
    end do  ! var
    !$omp end parallel do
    
    lon = lbr(1)
    lat = lbr(2)
    rad = lbr(3)
    
  subroutine vsop87d_lbr(tm,pl, lon,lat,rad, LBaccur,Raccur) …
  end subroutine vsop87d_lbr
  !*********************************************************************************************************************************
  
  
  
end module thesky_vsop
!***********************************************************************************************************************************