program SphereTri

   implicit none

   real(8), parameter :: PI = 3.141592653589793_8
   real(8), parameter :: DEG_PER_RAD = 180._8/PI
   real(8), parameter :: RAD_PER_DEG = PI/180._8
   real, parameter :: REARTH_M = 6370.     ! km
   real, parameter :: REARTH   = 1.        ! unit sphere
   real, parameter :: ninety = 90._8

   type point
     real :: lon
     real :: lat
   end type point
   type point_8
     real(8) :: lon
     real(8) :: lat
   end type point_8

   integer :: n, np1
   real(8)       :: lambda, theta, tmp
   real, pointer :: a, b, c
   real          :: s, E, Area
   real, target  :: dist(3)
   real          :: Angle(3)
   real(8)       :: x(3),y(3)
   real          :: coord(2)
   type(point) :: TriVtx(3)

   a => dist(1) ; b => dist(2) ; c => dist(3)
   write(*,*) ' '
   write(*,*) 'Exact area = ',real(RAD_PER_DEG*sin(RAD_PER_DEG),4)
 
   do
!--------------------------------------------------------
!  enter triangle vertices
!--------------------------------------------------------
     do n = 1,3
       write(*,*) ' '
       write(*,*) 'Enter vertex ',n,' (lon,lat)'
       read(*,*)  coord
       TriVtx(n)%lon = coord(1) ; TriVtx(n)%lat = coord(2)
     enddo

     if( TriVtx(1)%lon == TriVtx(2)%lon .and. TriVtx(1)%lat == TriVtx(2)%lat ) then
       exit
     endif

!--------------------------------------------------------
!  test mapping from polar to x,y
!--------------------------------------------------------
     do n = 1,3
       theta = RAD_PER_DEG*real(Trivtx(n)%lat,4)
       lambda = RAD_PER_DEG*real(Trivtx(n)%lon,4)
       tmp  = cos(theta)
       x(n) = tmp * cos(lambda)
       y(n) = tmp * sin(lambda)
     enddo

     write(*,*) ' '
     write(*,*) 'x = ',real(x(:),4)
     write(*,*) 'y = ',real(y(:),4)
     Area = poly_area( 3, real(x), real(y) )

     do n = 1,3
       if( n < 3 ) then
         np1 = n + 1
       else
         np1 = 1
       endif
       dist(n) = length( TriVtx(n),TriVtx(np1) )
       Angle(n) = distance( TriVtx(n),TriVtx(np1) )
       write(*,*) ' '
       write(*,'(''from ('',f7.2,'','',f7.2,'') to ('',f7.2,'','',f7.2,'')'')') &
                TriVtx(n)%lon,TriVtx(n)%lat,TriVtx(np1)%lon,TriVtx(np1)%lat
       write(*,*) 'distance = ',dist(n)
       write(*,*) 'distance = ',Angle(n)
     enddo
     s = .5*sum(dist(:))
     E = tan(.5*s)*tan(.5*(s-a))*tan(.5*(s-b))*tan(.5*(s-c))
     E = sqrt(E)
     E = 4.*tan(E)
     write(*,*) ' '
     write(*,*) 'Area on unit sphere = ',E
     Area = (REARTH_M**2)*E
     write(*,*) 'Triangle area = ',Area
     write(*,*) 'Equivalent sq = ',sqrt(Area)
   enddo

   CONTAINS

   real function length( pnt1_in, pnt2_in )

   implicit none

   type(point), intent(in) :: pnt1_in
   type(point), intent(in) :: pnt2_in

   real(8) :: c
   real    :: delLon
   real(8) :: delLon_rad
   type(point_8) :: pnt1_rad, pnt2_rad

!--------------------------------------------------------
!  switch points?
!--------------------------------------------------------
   if( pnt1_in%lon > pnt2_in%lon ) then
     pnt1_rad%lon = real(pnt2_in%lon,8) ; pnt1_rad%lat = real(pnt2_in%lat,8)
     pnt2_rad%lon = real(pnt1_in%lon,8) ; pnt2_rad%lat = real(pnt1_in%lat,8)
   else
     pnt1_rad%lon = real(pnt1_in%lon,8) ; pnt1_rad%lat = real(pnt1_in%lat,8)
     pnt2_rad%lon = real(pnt2_in%lon,8) ; pnt2_rad%lat = real(pnt2_in%lat,8)
   endif
   pnt1_rad%lat = ninety - pnt1_rad%lat
   pnt2_rad%lat = ninety - pnt2_rad%lat

!--------------------------------------------------------
!  convert from degrees to radians
!--------------------------------------------------------
   pnt1_rad%lon = rad_per_deg*pnt1_rad%lon
   pnt1_rad%lat = rad_per_deg*pnt1_rad%lat
   pnt2_rad%lon = rad_per_deg*pnt2_rad%lon
   pnt2_rad%lat = rad_per_deg*pnt2_rad%lat

   delLon = pnt1_in%lon - pnt2_in%lon
   delLon_rad = pnt1_rad%lon - pnt2_rad%lon
   if( abs(delLon) <= 180. ) then
     c = sin(pnt1_rad%lat)*sin(pnt2_rad%lat)*cos(delLon_rad) &
         + cos(pnt1_rad%lat)*cos(pnt2_rad%lat)
     length = real(acos(c),4)
   else
     Stop 'delLon > 180'
   endif

   end function length

   real function distance( pnt1_in, pnt2_in )

   implicit none

   type(point), intent(in) :: pnt1_in
   type(point), intent(in) :: pnt2_in

   real(8) :: c
   real    :: delLon
   real(8) :: delLon_rad
   type(point_8) :: pnt1_rad, pnt2_rad

!--------------------------------------------------------
!  transfer to local variables
!--------------------------------------------------------
   pnt1_rad%lon = real(pnt1_in%lon,8) ; pnt1_rad%lat = real(pnt1_in%lat,8)
   pnt2_rad%lon = real(pnt2_in%lon,8) ; pnt2_rad%lat = real(pnt2_in%lat,8)

!--------------------------------------------------------
!  convert from degrees to radians
!--------------------------------------------------------
   pnt1_rad%lon = rad_per_deg*pnt1_rad%lon
   pnt1_rad%lat = rad_per_deg*pnt1_rad%lat
   pnt2_rad%lon = rad_per_deg*pnt2_rad%lon
   pnt2_rad%lat = rad_per_deg*pnt2_rad%lat

   delLon = pnt1_in%lon - pnt2_in%lon
   delLon_rad = pnt1_rad%lon - pnt2_rad%lon
   if( abs(delLon) <= 180. ) then
     c = cos(pnt1_rad%lat)*cos(pnt2_rad%lat)*cos(delLon_rad) &
         + sin(pnt1_rad%lat)*sin(pnt2_rad%lat)
     distance = real(acos(c),4)
   else
     Stop 'delLon > 180'
   endif

   end function distance

   real FUNCTION poly_area( nv, x, y )
!---------------------------------------------------------------
!  calculate the area of polynomial with nv vertices
!---------------------------------------------------------------

!---------------------------------------------------------------
!  dummy arguments
!---------------------------------------------------------------
   integer, intent(in) :: nv
   real, intent(in)    :: x(nv)
   real, intent(in)    :: y(nv)

!---------------------------------------------------------------
!  local variables
!---------------------------------------------------------------
   integer :: i, im1, ip1
   real    :: wrk(nv)

   do i = 1,nv
     ip1 = mod( i,nv ) + 1
     im1 = i - 1
     if( im1 == 0 ) im1 = nv
     wrk(i) = (x(ip1) - x(im1))*y(i)
   end do

   poly_area = -.5*sum( wrk(:) )


   END FUNCTION poly_area

   end program SphereTri