PROGRAM SODANO_EXTERNAL

	!A wrapper for the Sodano calculation as used in mesointerp.f from my
	!Ph.D. research.  This is for use in Python programs and is run in its
	!entirety from the command line.

	!*(GFORTRAN) Compile with "gfortran sodano.f -o sodano.exe"
	!*(IFORT) Compile with "ifort-free sodano.f -o sodano.exe"

    IMPLICIT NONE
	DOUBLE PRECISION :: xstart,xend,ystart,yend
	DOUBLE PRECISION :: az12,az21
	DOUBLE PRECISION :: dist

	!Test
	xstart=-83.8887928283697
	xend=-83.8888585738692
	ystart=34.695898032828
	yend=34.6960287578698

	CALL ARGS(xstart,xend,ystart,yend)

	CALL SODANO(ystart,yend,xstart,xend,dist,az12,az21,3)

	!Produce output suitable for immediate on-the-fly reading by Python
	write(*,'(f10.4)')dist

	STOP
	END

	SUBROUTINE SODANO(B1,B2,L1,L2,dist,alpha12,alpha21,accuracy)
!	***********************************************************************!
!	This is a subroutine to calculate the distance (in km) and azimuth (in !
!	degrees) between two points on Earth using E.M. Sodano's inverse       !
!	geodetic calculations, assuming the WGS84 ellipsoid.                   !
!	The ellipticity (flatness, f) of the WGS84 ellipsoid is 1/298.257223563!
!	[(a-b)/a, where a is the semimajor axis (the equatorial radius) and b  !
!	is the semiminor axis (the polar radius)].  The eccentricity of        !
!	the WGS84 ellipsoid is found by sqrt(2f-f^2).  The equatorial          !
!	radius is 6378.1370 km. Direct inverse formulae are extended through   !
!	the f^3 and e'^6 terms for the greatest accuracy (accuracy=3).         !
!
!	REFERENCES:
!	National Imagery and Mapping Agency, 2000: Department of Defense
!		World Geodetic System 1984: Its Definition and
!		Relationships with Local Geodetic Systems. Tech. Rep. 
!		TR8350.2, 3d. ed., amendment 1, 3 January 2000, Bethesda,
!		MD, 175 pp. [Available from National Imagery and Mapping
!		Agency, ATTN: ISDFR, Mail Stop D-17, 4600 Sangamore Road,
!		Bethesda, MD 20816-5003.]
!		See also http://earth-info.nga.mil/GandG/pubs.html
!	Sodano, E.M, 1965: General non-iterative solution of the inverse and
!		direct geodetic problems. Bulletin Geodesique, 75, 69-89.
!
!	VARIABLES:
!	B1		geodetic latitude of point 1 in degrees
!	L1		geodetic longitude of point 1 in degrees
!	B2		geodetic latitude of point 2 in degrees
!	L2		geodetic longitude of point 2 in degrees
!	ao		semimajor axis of spheroid (equatorial radius of earth)
!	bo		semiminor axis of spheroid (polar radius of earth)
!	flat		spheroidal flattening (or ellipticity) = 1-bo/ao
!	ei2		second eccentricity squared = (ao^2-bo^2)/bo^2
!	eps		eccentricity
!	alpha12		initial azimuth from pt 1 to pt 2
!	alpha21		initial azimuth from pt 2 to pt 1
!	dist		distance between the two points
!	accuracy	Desired accuracy from Sodano (1965) - 1=App A, 2=App C, 3=App E
!	***********************************************************************
!	Christopher M. Godfrey				      29 September 2004
!	***********************************************************************
!
!	NOTES:
!	I recommend using accuracy=3 for all calculations.  Results closely match
!	measurements using a Garmin GPSIII Plus.  Sodano's inverse is exceedingly
!	accurate for short to medium length geodesics.
!
!	The formulae for the long geodesics do not appear to work correctly; perhaps
!	the problem is with my interpretation of the formulae in Sodano (1965)?  The
!	definition of a geodesic along the back-side of a spheroid is ambiguous
!	anyway.
!
!	Long geodesics over more than 134 degrees of longitude are inaccurate.  In
!	addition, long meridional arcs are inaccurate. Use the DISTANCE subroutine
!	instead for these long lines.

	IMPLICIT NONE
	DOUBLE PRECISION :: pi,ei2
	DOUBLE PRECISION :: flat,eps,phi
	DOUBLE PRECISION :: rmerid,rzonal,bearing
	DOUBLE PRECISION :: radius
	DOUBLE PRECISION :: B1,L1,B2,L2
	DOUBLE PRECISION :: lat1,lon1,lat2,lon2
	DOUBLE PRECISION :: ao,bo,dist
	DOUBLE PRECISION :: lambda,c,L,a,b,m
	DOUBLE PRECISION :: beta1,beta2
	DOUBLE PRECISION :: azimuth12,azimuth21
	DOUBLE PRECISION :: sinphi,cosphi,n,bminb
	DOUBLE PRECISION :: cosbeta1,cosbeta2
	DOUBLE PRECISION :: sinphinot,cosphinot,phinot
	DOUBLE PRECISION :: betanot,sigma,numdenom
	DOUBLE PRECISION :: Anot,Bnot,Cnot,Dnot
	DOUBLE PRECISION :: L2save,B2save,dist1
	DOUBLE PRECISION :: cotalpha12,cotalpha21
	DOUBLE PRECISION :: alpha12,alpha21
	INTEGER :: signL,signB1,signB2
	INTEGER :: accuracy,i
    DOUBLE PRECISION :: dsind,dtand,dcosd,dacosd,datand,dasind
    REAL :: tand,asind,atand

!	**************************************************
!	DEFINE CONSTANTS FOR WGS84 ELLIPSOID
	!pi=16*datan(1./5.)-4*datan(1./239.)	!Machin's formula for pi
	pi=3.14159265358979323846
	flat=1.0/(298.257223563)
	ao=6378.1370	!km
	bo=(1-flat)*ao	!km
	eps=((2.0*flat)-(flat**2.0))**(1./2.)
	ei2=(ao**2-bo**2)/bo**2		!Second eccentricity squared
	n=(ao-bo)/(ao+bo)
!	**************************************************

	!Error check
	if(B1.lt.-90.0.or.B1.gt.90.0.or.B2.lt.-90.0.or.B2.gt.90.0.or.L1&
	  .lt.-180.0.or.L1.gt.180.0.or.L2.lt.-180.0.or.L2.gt.180.0)then
	  write(*,'("Your latitudes and longitudes appear to be in a non-standard",/,&
	    "format.  Latitudes should be in the range -90 to 90 and longitudes",/,&
	    "should be in the range -180 to 180.  Please make the necessary",/,&
	    "adjustments and try again.")')
	  STOP "STOP 418"
	endif
	if(abs(L1-L2).gt.134.0)then
	  write(*,'("***WARNING***  Sodano''s inverse breaks down for nearly zonal",/,&
            "lines across more than 134 degrees of longitude. Try using the",/,&
	    "DISTANCE subroutine instead.")')
	endif 

	!DEFINE L
	L=(L2-L1)
	if(abs(L).gt.180.and.L.gt.0)then
	  L=(L2-L1)-360D0
	elseif(abs(L).gt.180.and.L.lt.0)then
	  L=(L2-L1)+360D0
	endif

	!DEFINE BETA
	!NOTE: There is no need to define beta using both tangent and cotangent
	!as Sodano uses on p.73. Both are mathematically equivalent.  This is
	!perhaps a consequence of early computer programming methods?
	beta1=datand(tand(B1)*(1.0-flat))
	beta2=datand(tand(B2)*(1.0-flat))

	a=dsind(beta1)*dsind(beta2)
	b=dcosd(beta1)*dcosd(beta2)
	cosphi=a+b*dcosd(L)

	!Appendix A
	if(accuracy.eq.1)then
	  sinphi=((dsind(L)*dcosd(beta2))**2+((dsind(beta2)*dcosd(beta1))-&
	    (dsind(beta1)*dcosd(beta2)*dcosd(L)))**2)**(0.5)
	endif

	!Appendix C
	if(accuracy.eq.2.or.accuracy.eq.3)then
	  if(B1.gt.0) signB1=1
	  if(B1.lt.0) signB1=-1
	  if(B2.gt.0) signB2=1
	  if(B2.lt.0) signB2=-1
	  if(abs(B1).gt.45)then
	    cosbeta1=dsind((90-signB1*B1)+signB1*2.*(dsind(90-signB1*B1))*(n+n**2+n**3)*sin(beta1))
	  else
	    cosbeta1=dcosd(beta1)
	  endif
	  if(abs(B2).gt.45)then
	    cosbeta2=dsind((90-signB2*B2)+signB2*2.*(dsind(90-signB2*B2))*(n+n**2+n**3)*sin(beta2))
	  else
	    cosbeta2=dcosd(beta2)
	  endif
	  bminb=((B2-B1)*pi/180.)+2*dsind(B2-B1)*((n+n**2+n**3)*a-(n-n**2+n**3)*b) !radians
	  sinphi=((dsind(L)*cosbeta2)**2+(dsin(bminb)+2.*cosbeta2*&
	    dsind(beta1)*((dsind(L/2.))**2))**2)**(0.5)
	endif

	!Determine smaller absolute value of sinphi or cosphi for reference angle phi
	if(abs(cosphi).le.abs(sinphi))then
	  phi=dacosd(cosphi)*pi/180D0	!radians
	else
	  phi=dasind(sinphi)*pi/180D0	!radians
	endif

	if(dsin(phi).ne.0)then
	  c=(b*dsind(L))/(dsin(phi))
	else
	  c=0.0
	endif
	m=1.-c**2

	!Appendix E
	if(accuracy.eq.3.and.dsin(phi).ne.0)then
	  lambda=L*pi/180.+c*(&
	    ((flat+flat**2+flat**3)*phi)&
	    +a*(-(((flat**2)/2.)+flat**3)*dsin(phi)-(flat**2+4*flat**3)*(phi**2)*(1./dsin(phi))&
	      +((3*flat**3)/2.)*(phi**3)*(1./dsin(phi))*(1./dtan(phi)))&
	    +m*(-(((5*flat**2)/4.)+3*flat**3)*phi+(((flat**2)/4.)+((flat**3)/2.))*dsin(phi)*dcos(phi)&
	      +(flat**2+4*flat**3)*(phi**2)*(1./dtan(phi))-((flat**3)/2.)*(phi**3)*&
	        ((1./dsin(phi))**2)-(flat**3)*(phi**3)*((1./dtan(phi))**2))&
	    +(a**2)*((flat**3)*phi+((flat**3)/2.)*dsin(phi)*dcos(phi)+(flat**3)*(phi**3)*((1./dsin(phi))**2))&
	    +(m**2)*(((31*flat**3)/16.)*phi-((9*flat**3)/16.)*dsin(phi)*dcos(phi)+&
	        ((flat**3)/2.)*phi*((dcos(phi))**2)&
	      -((9*flat**3)/2.)*(phi**2)*(1./dtan(phi))+((flat**3)/8.)*dsin(phi)*((dcos(phi))**3)&
	      +((flat**3)/2.)*(phi**3)*((1./dsin(phi))**2)+(2*flat**3)*(phi**3)*((1./dtan(phi))**2))&
	    +a*m*((flat**3)*dsin(phi)-((3*flat**3)/2.)*phi*dcos(phi)+((9*flat**3)/2.)*(phi**2)*(1./dsin(phi))&
	      -((flat**3)/2.)*dsin(phi)*((dcos(phi))**2)-((7*flat**3)/2.)*(phi**3)*&
	        (1./dsin(phi))*(1./dtan(phi)))&
	    )					!radians

	  !Phinot is determined using lambda instead of L; components are defined above in Appendix C
	  cosphinot=a+b*dcos(lambda)
	  sinphinot=((dsin(lambda)*cosbeta2)**2+(dsin(bminb)+2.*cosbeta2*&
	    dsind(beta1)*((dsin(lambda/2.))**2))**2)**(0.5)

	  !Determine smaller absolute value of sinphinot or cosphinot for reference angle phinot
	  if(abs(cosphinot).le.abs(sinphinot))then
	    phinot=dacosd(cosphinot)*pi/180D0	!radians
	  else
	    phinot=dasind(sinphinot)*pi/180D0	!radians
	  endif
	endif


	!DISTANCE
	if(accuracy.eq.1.or.accuracy.eq.2.and.dsin(phi).ne.0)then
	  dist=bo*(&
	    ((1+flat+flat**2)*phi)+&
	    a*((flat+flat**2)*dsin(phi)-(((flat**2)/2.)*(phi**2)*(1./dsin(phi))))+&
	    m*(-((flat+flat**2)/2.)*phi-((flat+flat**2)/2.)*dsin(phi)*dcos(phi)+&
	      ((flat**2)/2.0)*(phi**2)*(1./dtan(phi)))+&
	    (a**2)*(-((flat**2)/2.)*dsin(phi)*dcos(phi))+&
	    (m**2)*(((flat**2)/16.)*phi+((flat**2)/16.)*dsin(phi)*dcos(phi)-&
	      ((flat**2)/2.)*(phi**2)*(1./dtan(phi))-((flat**2)/8.)*dsin(phi)*&
	      ((dcos(phi))**3))+&
	    a*m*(((flat**2)/2.)*(phi**2)*(1./dsin(phi))+((flat**2)/2.)*dsin(phi)*((dcos(phi))**2))&
	    )					!km
	elseif(accuracy.eq.1.or.accuracy.eq.2)then
	  dist=0.0 				!km
	endif
	if(accuracy.eq.3.and.dsin(phinot).ne.0)then
	  !Error check to prevent -nan
	  numdenom=(b*dsin(lambda))/(dsin(phinot))
	  if(numdenom.lt.-1.0)then
	    betanot=180D0
	  elseif(numdenom.gt.1.0)then
	    betanot=0D0
	  else
	    betanot=dacosd((b*dsin(lambda))/(dsin(phinot)))
	  endif
	  !Error check to prevent -nan
	  if(betanot.eq.0D0)then
	    sigma=pi/2.0
	  else
	    if((((2*a)/((dsind(betanot))**2))-dcos(phinot)).gt.1D0)then
	      sigma=0D0
	    else
	      sigma=(dacos(((2*a)/((dsind(betanot))**2))-dcos(phinot)))/2.0
	    endif
	  endif
	  Anot=1+(ei2/4.)*((dsind(betanot))**2)-((3*ei2**2)/64.)*((dsind(betanot))**4)&
	    +((5*ei2**3)/256.)*((dsind(betanot))**6)
	  Bnot=(ei2/4.)*((dsind(betanot))**2)-((ei2**2)/16.)*((dsind(betanot))**4)&
	    +((15*ei2**3)/512.)*((dsind(betanot))**6)
	  Cnot=((ei2**2)/128.)*((dsind(betanot))**4)-((3*ei2**3)/512.)*((dsind(betanot))**6)
	  Dnot=((ei2**3)/1536.)*((dsind(betanot))**6)
	  dist=bo*((Anot*phinot)+(Bnot*dsin(phinot)*dcos(2*sigma))-&
	    (Cnot*dsin(2*phinot)*dcos(4*sigma))+(Dnot*dsin(3*phinot)*dcos(6*sigma)))
	elseif(accuracy.eq.3)then
	  dist=0.0
	endif

	!AZIMUTH
	if(dist.ne.0.0)then
	  if(accuracy.eq.1.or.accuracy.eq.2)then
	    lambda=L*pi/180.+c*(&
	      ((flat+flat**2)*phi)&
	      +a*(-((flat**2)/2.)*dsin(phi)-(flat**2)*(phi**2)*(1/dsin(phi)))&
	      +m*(-((5*flat**2)/4.)*phi+((flat**2)/4.)*dsin(phi)*dcos(phi)+&
	        (flat**2)*(phi**2)*(1./dtan(phi)))&
	      )					!radians
	  endif


	  !Meridional arcs
	  if(L1.eq.L2)then
	    cotalpha12=(dsind(beta2)*dcosd(beta1)-dcos(lambda)*dsind(beta1)&
	      *dcosd(beta2))
	    cotalpha21=(dsind(beta2)*dcosd(beta1)*dcos(lambda)-dsind(beta1)&
	      *dcosd(beta2))
	    if(B1.gt.B2)then
	      alpha12=180.0
	      alpha21=0.0
	    endif
	    if(B2.gt.B1)then
	      alpha12=0.0
	      alpha21=180.0
	    endif
	  else
	    if(accuracy.eq.1.and.dsin(phi).ne.0)then	!Appendix A
	      cotalpha12=(dsind(beta2)*dcosd(beta1)-dcos(lambda)*dsind(beta1)&
	        *dcosd(beta2))/(dsin(lambda)*dcosd(beta2))
	      cotalpha21=(dsind(beta2)*dcosd(beta1)*dcos(lambda)-dsind(beta1)&
	        *dcosd(beta2))/(dsin(lambda)*dcosd(beta1))
	    endif

	    if(accuracy.eq.2.or.accuracy.eq.3.and.dsin(phi).ne.0)then	!Appendix C and E
	      cotalpha12=(dsin(bminb)+2*dcosd(beta2)*dsind(beta1)*&
	        ((dsin(lambda/2.))**2))/(dcosd(beta2)*dsin(lambda))
	      cotalpha21=(dsin(bminb)-2*dcosd(beta1)*dsind(beta2)*&
	        ((dsin(lambda/2.))**2))/(dcosd(beta1)*dsin(lambda))
	    endif

	    alpha12=atand(1./cotalpha12)
	    alpha21=atand(1./cotalpha21)
	  endif

	  !Assign quadrants
	  if(L.gt.0.0)then
	    if(cotalpha12.gt.0)then
	      !Quadrant I
	      do i=1,4
	        if(alpha12.lt.0.0) alpha12=alpha12+90.0
	        if(alpha12.ge.90.0) alpha12=alpha12-90.0
	      enddo
	    else
	      !Quadrant II
	      do i=1,5
	        if(alpha12.lt.90.0) alpha12=alpha12+90.0
	        if(alpha12.ge.180.0) alpha12=alpha12-90.0
	      enddo
	    endif
	    if(cotalpha21.gt.0)then
	      !Quadrant III
	      do i=1,6
	        if(alpha21.lt.180.0) alpha21=alpha21+90.0
	        if(alpha21.ge.270.0) alpha21=alpha21-90.0
	      enddo
	    else
	      !Quadrant IV
	      do i=1,7
	        if(alpha21.lt.270.0) alpha21=alpha21+90.0
	      enddo
	    endif
	  endif
	  if(L.lt.0.0)then
	    if(cotalpha12.gt.0)then
	      !Quadrant III
	      do i=1,6
	        if(alpha12.lt.180.0) alpha12=alpha12+90.0
	        if(alpha12.ge.270.0) alpha12=alpha12-90.0
	      enddo
	    else
	      !Quadrant IV
	      do i=1,7
	        if(alpha12.lt.270.0) alpha12=alpha12+90.0
	      enddo
	    endif
	    if(cotalpha21.gt.0)then
	      !Quadrant I
	      do i=1,4
	        if(alpha21.lt.0.0) alpha21=alpha21+90.0
	        if(alpha21.ge.90.0) alpha21=alpha21-90.0
	      enddo
	    else
	      !Quadrant II
	      do i=1,5
	        if(alpha21.lt.90.0) alpha21=alpha21+90.0
	        if(alpha21.ge.180.0) alpha21=alpha21-90.0
	      enddo
	    endif
	  endif
	else
	  alpha12=0.0
	  alpha21=0.0
	endif

	RETURN
	END

	SUBROUTINE ARGS(xstart,xend,ystart,yend)
	implicit none
	character(len=120) :: harg
	integer :: ierr, i, numarg,j
!	integer, external :: iargc
	logical :: output
	double precision :: xstart,xend,ystart,yend

	numarg = iargc()
   
	i = 1
	j = 1
   
	do while ( i <= numarg) 
	  call getarg(i, harg)
	  !print*, 'harg = ', trim(harg)

	  if (harg == "-h") then
	    call help
	  elseif (j == 1) then
	    read(harg,*)xstart
	    j = j + 1
	    i = i + 1
	  elseif (j == 2) then
	    read(harg,*)xend
	    j = j + 1
	    i = i + 1
	  elseif (j == 3) then
	    read(harg,*)ystart
	    j = j + 1
	    i = i + 1
	  elseif (j == 4) then
	    read(harg,*)yend
	    j = j + 1
	    i = i + 1
	  else
	    j=j+1
	    i=i+1
	  endif
	enddo

	if(j>5)then
	  write(*,'(/,"TOO MANY ARGUMENTS!!")')
	  call help
	endif

	if(j<5)then
	  write(*,'(/,"NOT ENOUGH ARGUMENTS!!")')
	  call help
	endif

	RETURN
	END

	SUBROUTINE HELP
	implicit none
	character(len=120) :: cmd
	call getarg(0, cmd)

	write(*,'(/,"Usage: ", A, " xstart xend ystart yend")') trim(cmd)
	write(*,'(8x, "-h       : Print this help message and exit.")')
	write(*,'(8x, "xstart   : Longitude of starting point.")')
	write(*,'(8x, "xend     : Longitude of ending point.")')
	write(*,'(8x, "ystart   : Latitude of starting point.")')
	write(*,'(8x, "yend     : Latitude of ending point.")')
	STOP

	RETURN
	END

	DOUBLE PRECISION FUNCTION dsind(degs)
	DOUBLE PRECISION :: degs,pi
    pi=3.14159265358979323846
	dsind=dsin(degs*(pi/180.))
	RETURN
	END

	DOUBLE PRECISION FUNCTION dcosd(degs)
	DOUBLE PRECISION :: degs,pi
    pi=3.14159265358979323846
	dcosd=dcos(degs*(pi/180.))
	RETURN
	END

	DOUBLE PRECISION FUNCTION datand(arg)
	DOUBLE PRECISION :: arg,pi
    !Returns the double precision inverse tangent of the argument in degrees
    pi=3.14159265358979323846
	datand=datan(arg)*180./pi
	RETURN
	END

	DOUBLE PRECISION FUNCTION dtand(degs)
	DOUBLE PRECISION :: degs,pi
    pi=3.14159265358979323846
	dtand=dtan(degs*(pi/180.))
	RETURN
	END

	REAL FUNCTION tand(degs)
	DOUBLE PRECISION :: degs,pi
    pi=3.14159265358979323846
	tand=tan(degs*(pi/180.))
	RETURN
	END

	DOUBLE PRECISION FUNCTION dacosd(arg)
	DOUBLE PRECISION :: arg,pi
    !Returns the inverse cosine of the argument in degrees
    pi=3.14159265358979323846
	dacosd=dacos(arg)*180./pi
	RETURN
	END

	DOUBLE PRECISION FUNCTION dasind(arg)
	DOUBLE PRECISION :: arg,pi
    !Returns the inverse sine of the argument in degrees
    pi=3.14159265358979323846
	dasind=dasin(arg)*180./pi
	RETURN
	END

	REAL FUNCTION atand(arg)
	DOUBLE PRECISION :: arg,pi
    !Returns the real inverse tangent of the argument in degrees
    pi=3.14159265358979323846
	atand=atan(arg)*180./pi
	RETURN
	END

	REAL FUNCTION asind(arg)
	DOUBLE PRECISION :: arg,pi
    !Returns the real inverse sine of the argument in degrees
    pi=3.14159265358979323846
	asind=asin(arg)*180./pi
	RETURN
	END