BearingAngle Function

public function BearingAngle(pointA, pointB) result(angle)

Uses

Find the bearing angle (radians) between two points in a 2D space, as the angle measured in the clockwise direction from the north line with A as the origin to the line segment AB. Assume radians as input unit for geodetic reference system.

bearing angle on spheroid

Arguments

Type IntentOptional Attributes Name
type(Coordinate), intent(in) :: pointA

origin point

type(Coordinate), intent(in) :: pointB

ending point

Return Value real(kind=float)

bearing angle (radians)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: a1
real(kind=float), public :: a2
real(kind=float), public :: b1
real(kind=float), public :: b2
real(kind=float), public :: deltaLongitude
real(kind=float), public :: x
real(kind=float), public :: y

Source Code

FUNCTION BearingAngle &
!
(pointA, pointB) &
!
RESULT (angle)

USE Units, ONLY: &
!Imported parametes:
pi

IMPLICIT NONE

!Arguments with intent in:
TYPE (Coordinate), INTENT(IN) :: pointA !! origin point 
TYPE (Coordinate), INTENT(IN) :: pointB !! ending point 


!Local declarations:
REAL (KIND = float) :: angle !!bearing angle (radians)
REAL (KIND = float) :: a1, a2, b1, b2, x, y, deltaLongitude

!----------------------end of declarations-------------------------------------
a1 = pointA %  easting !longitude of point A
a2 = pointA %  northing !latitude of point A

b1 = pointB %  easting !longitude of point B
b2 = pointB %  northing !latitude of point B


IF ( pointA %  system % system == geodetic) THEN !!bearing angle on spheroid
    
    deltaLongitude = b1 - a1
    
    y = COS ( a2 ) * SIN (b2) - SIN ( a2 ) * COS ( b2 ) * COS ( deltaLongitude )
    
    x = SIN ( deltaLongitude ) * COS ( b2 )
    
    angle = ATAN2 ( x, y ) !arguments are inverted to find angle of the 
                        !vector with respect to the North
    IF ( angle < 0.) angle = angle + 2. * pi
    
ELSE  !points are in a projected coordinate system
    
    y = b2 - a2
    
    x = b1 - a1
    
    angle = ATAN2 ( x, y )
    IF ( angle < 0.) angle = angle + 2. * pi
    
END IF


RETURN
END FUNCTION BearingAngle