module piab
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    real(8),parameter::a=6378137.000000D0 !Equatorial radius
    real(8),parameter::b=6356752.314140D0 !Polar radius
end module piab


program f90_DIST
    use piab
    implicit none
    real(8)::d1,d2
    real(8)::lat1,lon1
    real(8)::lat2,lon2
    character(len=20)::dummy1,dummy2,dummy3,dummy4

    call getarg(1,dummy1)
    call getarg(2,dummy2)
    call getarg(3,dummy3)
    call getarg(4,dummy4)
    read(dummy1,*) lon1
    read(dummy2,*) lat1
    read(dummy3,*) lon2
    read(dummy4,*) lat2

    d1=DLA(lon1,lat1,lon2,lat2)
    d2=DHB(lon1,lat1,lon2,lat2)

    write(6,'("d1=",f10.3,"km (Lambert-Andyer)")') d1/1000.0D0
    write(6,'("d2=",f10.3,"km (Hubeny)")') d2/1000.0D0

stop

contains

    real(8) function DLA(lon1,lat1,lon2,lat2)
        !Lambert-Andyer's formula
        real(8),intent(in)::lat1,lon1,lat2,lon2
        real(8)::IA,LA
        real(8)::IB,LB
        real(8)::phiA,phiB,XX,delta
        real(8)::rho
        real(8)::FF

        IA=lat1/180.0*pi
        LA=lon1/180.0*pi
        IB=lat2/180.0*pi
        LB=lon2/180.0*pi

        FF=(a-b)/a
        phiA=atan((1.0-FF)*tan(IA))
        phiB=atan((1.0-FF)*tan(IB))
        XX=acos(sin(phiA)*sin(phiB)+cos(phiA)*cos(phiB)*cos(LA-LB))
        delta=FF/8.0*&
        (&
         (sin(XX)-XX)*(sin(phiA)+sin(phiB))**2.0D0/(cos(0.5D0*XX))**2.0D0&
        -(sin(XX)+XX)*(sin(phiA)-sin(phiB))**2.0D0/(sin(0.5D0*XX))**2.0D0&
        )
        DLA=a*(XX+delta)
    end function DLA


    real(8) function DHB(lon1,lat1,lon2,lat2)
        !Hubeny's formula
        real(8),intent(in)::lat1,lon1,lat2,lon2
        real(8)::x1,y1
        real(8)::x2,y2
        real(8)::dy,dx,my,MM,NN,WW,e

        x1=lon1/180.0*pi
        y1=lat1/180.0*pi
        x2=lon2/180.0*pi
        y2=lat2/180.0*pi

        e=sqrt((a*a-b*b)/a/a)
        my=0.5*(y1+y2)
        WW=sqrt(1.0-e*e*sin(my)*sin(my))
        NN=a/WW
        MM=a*(1.0-e*e)/WW/WW/WW
        dy=y1-y2
        dx=x1-x2
        DHB=sqrt((dy*MM)**2.0D0+(dx*NN*cos(my))**2.0D0)
    end function DHB


end program f90_DIST