angles.f90

Go to the documentation of this file.

!> \file angles.f90  Procedures to handle angles
 
 
!  Copyright (c) 2002-2025  Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
!   
!  This file is part of the libSUFR package, 
!  see: http://libsufr.sourceforge.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 handle angles
 
module sufr_angles
  use sufr_kinds, only: double
  
  implicit none
  save
  private :: double
  
contains
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in radians between 0 and 2pi
  !!
  !! \param  ang  Input angle (radians)
  !! \retval rev  Angle in radians between 0 and 2pi
  
  pure function rev(ang)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2
    
    implicit none
    real(double), intent(in) :: ang
    real(double) :: rev
    
    rev = ang - dble(floor(ang/pi2)) * pi2
    
  pure function rev(ang) …
  end function rev
  !*********************************************************************************************************************************
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in radians between 0 and 2pi  --  single precision
  !!
  !! \param  ang  Input angle (radians)
  !! \retval rrev  Angle in radians between 0 and 2pi
  
  pure function rrev(ang)
    use sufr_constants, only: rpi2
    
    implicit none
    real, intent(in) :: ang
    real :: rrev
    
    rrev = ang - real(floor(ang/rpi2)) * rpi2
    
  pure function rrev(ang) …
  end function rrev
  !*********************************************************************************************************************************
  
 
 
  !*********************************************************************************************************************************
  !> \brief  Returns angle in degrees between 0 and 360
  !!
  !! \param  ang     Input angle (degrees)
  !! \retval rev360  Angle in degrees between 0 and 360
  
  pure function rev360(ang)
    use sufr_kinds, only: double
    
    implicit none
    real(double), intent(in) :: ang
    real(double) :: rev360
    
    rev360 = ang - dble(floor(ang/360.d0)) * 360.d0
    
  pure function rev360(ang) …
  end function rev360
  !*********************************************************************************************************************************
  
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in radians between -pi and pi
  !!
  !! \param  ang   Input angle (radians)
  !! \retval rev2  Angle in radians between -pi and pi
  
  pure function rev2(ang)
    use sufr_kinds, only: double
    use sufr_constants, only: pi,pi2
    
    implicit none
    real(double), intent(in) :: ang
    real(double) :: rev2
    
    rev2 = ang - dble(floor(ang/pi2)) * pi2
    if(rev2.gt.pi) rev2 = rev2 - pi2
    
  pure function rev2(ang) …
  end function rev2
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in radians between -pi and pi  --  single precision
  !!
  !! \param  ang    Input angle (radians)
  !! \retval rrev2  Angle in radians between -pi and pi
  
  pure function rrev2(ang)
    use sufr_constants, only: rpi,rpi2
    
    implicit none
    real, intent(in) :: ang
    real :: rrev2
    
    rrev2 = ang - real(floor(ang/rpi2)) * rpi2
    if(rrev2.gt.rpi) rrev2 = rrev2 - rpi2
    
  pure function rrev2(ang) …
  end function rrev2
  !*********************************************************************************************************************************
  
  
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in radians between cen-pi and cen+pi
  !!
  !! \param  ang   Input angle (radians)
  !! \param  cen   'Central value' (radians)
  !! \retval revc  Angle in radians between cen-pi and cen+pi
  
  pure function revc(ang,cen)
    use sufr_kinds, only: double
    use sufr_constants, only: pi,pi2
    
    implicit none
    real(double), intent(in) :: ang,cen
    real(double) :: revc
    
    revc = ang - floor(ang/pi2) * pi2
    if(cen.ge.pi) then
       if(revc.lt.cen-pi) revc = revc + pi2
    else
       if(revc.gt.cen+pi) revc = revc - pi2
    end if
    
  pure function revc(ang,cen) …
  end function revc
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns time in hours between 0 and 24
  !!
  !! \param  tm  Input time (hours)
  !! \retval rv  Time in hours between 0 and 24
  
  pure function rv(tm)
    use sufr_kinds, only: double
    implicit none
    real(double), intent(in) :: tm
    real(double) :: rv
    
    rv = tm - floor(tm/24.d0) * 24
    
  pure function rv(tm) …
  end function rv
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns time in hours between -12 and 12
  !!
  !! \param  tm    Input time (hours)
  !! \retval rv12  Time in hours between -12 and 12
  
  pure function rv12(tm)
    use sufr_kinds, only: double
    implicit none
    real(double), intent(in) :: tm
    real(double) :: rv12
    
    rv12 = tm - floor(tm/24.d0) * 24
    if(rv12.gt.12.d0) rv12 = rv12 - 24.d0
    
  pure function rv12(tm) …
  end function rv12
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Returns angle in degrees between -180 and 180
  !!
  !! \param  ang    Input angle (degrees)
  !! \retval rv180  Angle in degrees between -180 and 180
  
  pure function rv180(ang)
    use sufr_kinds, only: double
    implicit none
    real(double), intent(in) :: ang
    real(double) :: rv180
    
    rv180 = ang - floor(ang/360.d0) * 360
    if(rv180.gt.180.d0) rv180 = rv180 - 360.d0
    
  pure function rv180(ang) …
  end function rv180
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Calculate the angular separation between two objects with given longitudes and latitudes.
  !!
  !! \param l1  Longitude of object 1 (rad).
  !! \param l2  Longitude of object 2 (rad).
  !! \param b1  Latitude of object 1 (rad).
  !! \param b2  Latitude of object 2 (rad).
  !! 
  !! \retval asep  Angular separation between the two objects (rad).
  !!
  !! \see Celestial Mechanics in a Nutshell (CMiaNS), Sect.5.3: Angular distance and phase angle between two
  !! celestial objects (CMiaNS.sf.net).
  
  pure function asep(l1,l2, b1,b2)
    use sufr_kinds, only: double
    
    implicit none
    real(double), intent(in) :: l1,l2, b1,b2
    real(double) :: asep,dl,db,bb
    
    dl = rev2(l2-l1)
    asep = acos(min(max( sin(b1)*sin(b2) + cos(b1)*cos(b2)*cos(dl), -1.d0),1.d0))  ! Min/max: catch round-off errors leading to |arg|>1.
    
    if(asep.lt.3.d-3) then  ! for angles < 10'
       bb = rev2((b1+b2)/2.d0)
       db = rev2(b2-b1)
       asep = sqrt((dl*cos(bb))**2 + db**2)
    end if
    
  pure function asep(l1,l2, b1,b2) …
  end function asep
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief Calculates the position angle of object 2 with respect to object 1, COUNTERCLOCKWISE from the north
  !!
  !! \param lon1  Longitude of object 1 - RA if measuring from the north (rad)
  !! \param lon2  Longitude of object 2 - RA if measuring from the north (rad)
  !! \param lat1  Latitude of object 1 - Dec if measuring from the north (rad)
  !! \param lat2  Latitude of object 2 - Dec if measuring from the north (rad)
  !! 
  !! \retval calpa  Position angle of object 2 wrt object 1 (rad)
  
  pure function calpa(lon1,lon2,lat1,lat2)
    implicit none
    real(double), intent(in) :: lon1,lon2,lat1,lat2
    real(double) :: calpa,dlon
    
    dlon  = lon2-lon1
    calpa = atan2( sin(dlon),  cos(lat1)*tan(lat2) - sin(lat1)*cos(dlon) )
    
  pure function calpa(lon1,lon2,lat1,lat2) …
  end function calpa
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a position angle to one of eight English two-letter abbreviations (NE, SW)
  !!
  !! \param  pa        Position angle (radians, N=0)
  !! \retval pastr_en  English two-letter wind-direction abbreviation
  
  pure function pastr_en(pa)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2
    
    implicit none
    real(double), intent(in) :: pa
    character :: pastr_en*(2),pas(9)*(2)
    
    pas = (/'N ','NE','E ','SE','S ','SW','W ','NW','N '/)  ! You can use trim()
    pastr_en = pas(ceiling( rev(pa)/pi2 * 8 + 0.5d0 ))      ! 1-9
    
  pure function pastr_en(pa) …
  end function pastr_en
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a position angle to one of eight full Dutch strings (noordoosten, noorden)
  !!
  !! \param  pa        Position angle (radians, N=0)
  !! \retval pastr_nl  Full Dutch wind-direction string
  
  pure function pastr_nl(pa)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2
    
    implicit none
    real(double), intent(in) :: pa
    character :: pastr_nl*(11),pas(9)*(11)
    
    pas = (/'noorden    ','noordoosten','oosten     ','zuidoosten ','zuiden     ','zuidwesten ','westen     ','noordwesten', &
         'noorden    '/)
    pastr_nl = pas(ceiling( rev(pa)/pi2 * 8 + 0.5d0 ))  ! 1-9
    
  pure function pastr_nl(pa) …
  end function pastr_nl
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a position angle to one of eight Dutch abbreviations (NO,ZW)
  !!
  !! \param  pa         Position angle (radians, N=0)
  !! \retval pastr_nls  Dutch wind-direction abbreviation
  
  pure function pastr_nls(pa)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2
    
    implicit none
    real(double), intent(in) :: pa
    character :: pastr_nls*(2),pas(9)*(2)
    
    pas = (/'N ','NO','O ','ZO','Z ','ZW','W ','NW','N '/)  ! You can use trim()
    pastr_nls = pas(ceiling( rev(pa)/pi2 * 8 + 0.5d0 ))     ! 1-9
    
  pure function pastr_nls(pa) …
  end function pastr_nls
  !*********************************************************************************************************************************
  
  
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 16 English wind directions (northnortheast, westsouthwest)
  !!
  !! \param  wd        Wind direction/azimuth (radians, N=0)
  !! \retval wdstr_en  Full English wind direction
  
  pure function wdstr_en(wd)
    use sufr_kinds, only: double
    use sufr_constants, only: r2d
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_en*(14),wds(0:15)*(14)
    
    wds = [character(len=14) :: 'north','northnortheast','northeast','eastnortheast','east','eastsoutheast','southeast', &
         'southsoutheast','south','southsouthwest','southwest','westsouthwest','west','westnorthwest','northwest','northnorthwest']
    
    wdstr_en = wds(mod( nint(wd*r2d/22.5) + 32, 16 ))  ! Mod, so that -1 -> 15
    
  pure function wdstr_en(wd) …
  end function wdstr_en
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 16 three-letter English wind-direction abbreviations (NNE, WSW)
  !!
  !! \param  wd         Wind direction/azimuth (radians, N=0)
  !! \retval wdstr_ens  Three-letter English wind-direction abbreviation
  
  pure function wdstr_ens(wd)
    use sufr_kinds, only: double
    use sufr_constants, only: r2d
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_ens*(3),wds(0:15)*(3)
    
    wds = (/'N  ','NNE','NE ','ENE','E  ','ESE','SE ','SSE','S  ','SSW','SW ','WSW','W  ','WNW','NW ','NNW'/)
    wdstr_ens = wds(mod( nint(wd*r2d/22.5) + 32, 16 ))  ! Mod, so that -1 -> 15
    
  pure function wdstr_ens(wd) …
  end function wdstr_ens
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 16 three-letter Dutch wind-direction abbreviations (NNO, WZW)
  !!
  !! \param  wd         Wind direction/azimuth (radians, N=0)
  !! \retval wdstr_nls  Three-letter Dutch wind-direction abbreviation
  
  pure function wdstr_nls(wd)
    use sufr_kinds, only: double
    use sufr_constants, only: r2d
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_nls*(3),wds(0:15)*(3)
    
    wds = (/'N  ','NNO','NO ','ONO','O  ','OZO','ZO ','ZZO','Z  ','ZZW','ZW ','WZW','W  ','WNW','NW ','NNW'/)
    wdstr_nls = wds(mod( nint( wd*r2d/22.5) + 32, 16 ))  ! Mod, so that -1 -> 15
    
  pure function wdstr_nls(wd) …
  end function wdstr_nls
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 16 full Dutch wind-direction strings (noordnoordoost, zuidwest)
  !!
  !! \param  wd        Wind direction (radians, N=0)
  !! \retval wdstr_nl  Full Dutch wind-direction string
  
  pure function wdstr_nl(wd)
    use sufr_kinds, only: double
    use sufr_constants, only: r2d
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_nl*(14),wds(0:15)*(14)
    
    wds = (/ &
         'noord         ', 'noordnoordoost', 'noordoost     ', 'oostnoordoost ', &
         'oost          ', 'oostzuidoost  ', 'zuidoost      ', 'zuidzuidoost  ', &
         'zuid          ', 'zuidzuidwest  ', 'zuidwest      ', 'westzuidwest  ', &
         'west          ', 'westnoordwest ', 'noordwest     ', 'noordnoordwest'/)
    
    wdstr_nl = wds(mod( nint( wd*r2d/22.5) + 32, 16 ))  ! Mod, so that -1 -> 15
    
  pure function wdstr_nl(wd) …
  end function wdstr_nl
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 8 two-letter Dutch secondary wind-direction abbreviations (NO, Z).
  !!         Wrapper for pastr_nls().
  !!
  !! \param  wd          Wind direction (radians, N=0)
  !! \retval wdstr_nls2  Dutch two-letter wind-direction abbreviation
  
  pure function wdstr_nls2(wd)
    use sufr_kinds, only: double
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_nls2*(2)
    
    wdstr_nls2 = pastr_nls(wd)
    
  pure function wdstr_nls2(wd) …
  end function wdstr_nls2
  !*********************************************************************************************************************************
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Converts a wind direction/azimuth to one of 8 full Dutch secondary wind-direction (noordoost, zuid).
  !!         Derived from pastr_nl().
  !!
  !! \param  wd         Wind direction (radians, N=0)
  !! \retval wdstr_nl8  String with full Dutch wind-direction
  
  pure function wdstr_nl8(wd)
    use sufr_kinds, only: double
    use sufr_constants, only: pi2
    
    implicit none
    real(double), intent(in) :: wd
    character :: wdstr_nl8*(9),wds(9)*(9)
    
    wds = (/'noord    ','noordoost','oost     ','zuidoost ','zuid     ','zuidwest ','west     ','noordwest', 'noord    '/)
    wdstr_nl8 = wds(ceiling( rev(wd)/pi2 * 8 + 0.5d0 ))  ! 1-9
    
  pure function wdstr_nl8(wd) …
  end function wdstr_nl8
  !*********************************************************************************************************************************
  
  
  
  
  
  
end module sufr_angles
!***********************************************************************************************************************************