PROGRAM gcarcdist
IMPLICIT NONE
REAL(KIND=8), PARAMETER :: pi=3.141592653589793
REAL(KIND=8), PARAMETER :: Rad = 6371.0       !Radius of Earth
REAL(KIND=8), PARAMETER :: ktm = 0.621371192  !km to miles
REAL(KIND=8) :: lat1, lon1, lat2, lon2    
REAL(KIND=8) :: v1, v2, v3
REAL(KIND=8) :: w1, w2, w3
REAL(KIND=8) :: rho, theta, phi, dtr, rtd
REAL(KIND=8) :: dotprod, magW, magV, arg, s, alpha


! Input two latitude and longitude coordinates on Earth Surface
lat1 = 0.0
lon1 = 0.0

lat2 = 0.0
lon2 = 45.0

! Some preliminaries
dtr = pi/180.0   !degree to radian conversion
rtd = 180.0/pi   !radian to degree conversion

! First lets convert the input lat lons to vectors in Cartesian space

!Position #1 (lat1, lon1)
rho = Rad
theta = lon1*dtr
phi = (90.0-lat1)*dtr

v1 = rho*SIN(phi)*COS(theta)
v2 = rho*SIN(phi)*SIN(theta)
v3 = rho*COS(phi)

!Position #2 (lat2, lon2)
rho = Rad
theta = lon2*dtr
phi = (90.0-lat2)*dtr

w1 = rho*SIN(phi)*COS(theta)
w2 = rho*SIN(phi)*SIN(theta)
w3 = rho*COS(phi)

!Now calculate dot product
dotprod = v1*w1 + v2*w2 + v3*w3
magV = SQRT(v1**2 + v2**2 + v3**2)
magW = SQRT(w1**2 + w2**2 + w3**2)

arg = dotprod/(magV*magW)
theta = ACOS(arg)
alpha = (ACOS(arg))*rtd  !angular distance

s = Rad*theta  !distance in km

!Output Results
write(*,*) "========================================================="
write(*,*) "Lat 1: ", lat1, "  Lon 1: ", lon1
write(*,*) "Lat 2: ", lat2, "  Lon 2: ", lon2
write(*,*) " "
write(*,*) "GCARC Distance: ", alpha, " (deg)"
write(*,*) "GCARC Distance: ", s, " (km)"
write(*,*) "GCARC Distance: ", s*ktm, " (miles)"
write(*,*) "========================================================="




END PROGRAM gcarcdist
