!
!  Spherical astronomy module
!
!  Copyright © 1996 - 2013, 2015-8 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  Munipack 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
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.



module astrosphere


  implicit none

  integer, parameter, private :: db = selected_real_kind(15)
  real(db), parameter, private :: rad =  57.29577951308232286464772_db

contains

  function gmst(jd)

    real(db) :: gmst
    real(db), intent(in) :: jd


!  Greenwich sidereal time in hours
!
!  jd is a full Julian date
!
!  The precision is better than 1 second.
!  According to Astronomical Almanac 2000.


    real(db) :: tu,t

    tu = (jd - 2451545.0_db)/36525.0_db
    t = 24110.54841_db + tu*(8640184.812866_db + tu*(0.093104_db-6.2e-6_db*tu))
    gmst = mod(t/3600.0_db + 24.0_db*(jd - aint(jd)) + 12.0_db,24.0_db)

  end function gmst


  function lmst(jd,longitude)

    real(db) :: lmst
    real(db), intent(in) :: jd, longitude


!  local sidereal time in hours
!
!  jd is a full Julian date
!  lambda is a longitude in degrees: -west ... +east
!
!  The precision is better than 1 second.
!  According to Astronomical Almanac 2000.

    lmst = mod(gmst(jd) + longitude/15.0_db,24.0_db)

  end function lmst



  function hangle(lmst,ra)

!   hour angle in degrees

    real(db) :: hangle
    real(db), intent(in) :: lmst, ra

    hangle = mod(lmst - ra,360.0_db)

  end function hangle


  subroutine eq2hor(ha, dec, latitude, az, elev)

    real(db), intent(in) :: ha,dec,latitude
    real(db), intent(out) :: az, elev


!
! equatorial to horizontal coordinates
!
!  all arguments in degrees
!

    real(db) :: sinh, cosh, sind, cosd, sinl, cosl, x,y,z,r

    sinh = sin(ha/RAD)
    cosh = cos(ha/RAD)
    sind = sin(dec/RAD)
    cosd = cos(dec/RAD)
    sinl = sin(latitude/RAD)
    cosl = cos(latitude/RAD)

    x = -cosh*cosd*sinl + sind*cosl
    y = -sinh*cosd
    z = cosh*cosd*cosl + sind*sinl

    r = sqrt(x**2 + y**2)
    if( abs(r) > epsilon(r) )then
       az = RAD*atan2(y,x)
    else
       az = 0.0_db
    end if
    if( az < 0_db ) az = az + 360.0_db
    elev = RAD*atan2(z,r)

  end subroutine eq2hor


  subroutine hor2eq(az, elev, latitude, ha,dec)

    real(db), intent(in) :: az,elev,latitude
    real(db), intent(out) :: ha, dec

!
!  horizontal to equatorial coordinates
!
!  all arguments in degrees
!

    real(db) :: sina, cosa, sine, cose, sinl, cosl, x, y, z, r

    sina = sin(az/RAD)
    cosa = cos(az/RAD)
    sine = sin(elev/RAD)
    cose = cos(elev/RAD)
    sinl = sin(latitude/RAD)
    cosl = cos(latitude/RAD)

    x = -cosa*cose*sinl + sine*cosl
    y = -sina*cose
    z = cosa*cose*cosl + sine*sinl

    r = sqrt(x**2 + y**2)
    if( abs(r) > epsilon(r) )then
       ha = RAD*atan2(y,x)
    else
       ha = 0.0_db
    endif
    dec = RAD*atan2(z,r)

  end subroutine hor2eq


  function refract(z)

    real(db) :: refract
    real(db), intent(in) :: z

!
!     compute refraction angle in degrees
!
!    Smart: Textbook on spherical astronomy
!
!    constants for pressure 760mmHg, 10deg C with
!    suffucient accuracy for z < 75 deg
!

    real(db) :: tanz

    tanz = tan(z/RAD)
    refract = (58.16_db*tanz - 0.067_db*tanz*tanz*tanz)/3600.0_db

  end function refract


  function airmass(z)

    real(db) :: airmass
    real(db), intent(in) :: z

    !
    !     compute airmass,
    !
    !      young&irvine: aj,72,945,(1967)
    !
    ! the airmass is limited on the given range of zenit distances

    real(db) :: secz

    if( 0 <= z .and. z < 86.5 ) then
       secz = 1.0_db/cos(z/RAD)
       airmass = secz*(1.0_db - 1.2e-3_db*(secz**2 - 1.0_db))
    else
       airmass = -1
    end if

  end function airmass

  function xairmass(jd,long,lat,ra,dec)

    real(db) :: xairmass
    real(db), intent(in) :: jd,long,lat,ra,dec
    real(db) :: t,h,ha,a

    t = lmst(jd,long)
    ha = hangle(15.0_db*t,ra)
    call eq2hor(ha,dec,lat,a,h)
    xairmass = airmass(90.0_db - h)

  end function xairmass


!  function longsun(jd,y)
  function longsun(d)

!    use trajd

!    real(db), intent(in) :: jd,y
    real(db), intent(in) :: d   ! days since 1. january
    real(db) :: longsun

   !  approx (!!!!) of length of the Sun
    longsun = mod(279.465 + 0.985647*d,360.0_db)
!    longsun = 279.465 + 0.985647*(jd - datjd(y,1.0_db,1.0_db))

  end function longsun


  function helcor(alpha,delta,ls)

    real(db), intent(in) ::  alpha,delta,ls
    real(db) ::  helcor

    ! heliocentric correction in days, angles in degrees

    helcor = 0.9174077_db*sin(alpha/rad)*cos(delta/rad) + &
         0.3979486_db*sin(delta/rad)
    helcor = helcor*sin(ls/rad) + cos(ls/rad)*cos(alpha/rad)*cos(delta/rad)
    helcor = -0.0057755_db*helcor

  end function helcor

  function phase(jd,min0,per)

    real(db), intent(in) ::  jd,min0,per
    real(db) ::  phase

    phase = mod(jd - min0,per) / per

    ! Phase is negative for jd < min0.

  end function phase

  subroutine propercoo(jd0,jd,a,d,pma,pmd,alpha,delta)

    real(db), intent(in) :: jd0,jd
    real(db), intent(in) :: a,d,pma,pmd
    real(db), intent(out) :: alpha,delta
    real(db) :: dt

    dt = (jd - jd0)/365.25_db
    alpha = a + dt*pma
    delta = d + dt*pmd

  end subroutine propercoo


end module astrosphere