44.2 Routines for Geodetic computations (WGS84)

The SEISAN software libaray include two subroutines for geodetic computations using the WGS84 ellipsoid:

vincenty_forward(lon,lat,azi,dist) and vincenty_inverse(lo,la,lo2,la2)

vincenty_forward(lon,lat,azi,dist) computes the lat and lon for a point at an azimuth and distance (meter) from the point given by input lon and lat. Output lon and lat are returned at the position of input lon and lat.

vincenty_inverse(lo,la,lo2,la2) computes the azimuth, back-azumith and distance between the two points lo,la,lo2,la2. Outout azimuth, back-azumith and distance (meters) are returned at the position of lo,la,lo2.

The code is based on the code at http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/ see also
http://www.movable-type.co.uk/scripts/latlong-vincenty.html

Please check for bugs before use.

Example:

C program example that show how to call the Vincenty routines:
      double precision lo,la,azi,dist,lo2,la2
      real lon,lat
      lon=12.8833
      lat=56.54
      azi=0.0D0
      dist=111000.0D0   ! meters
      lo=DBLE(lon)
      la=DBLE(lat)
      lo2=DBLE(lon)
      la2=DBLE(lat)
      write(*,*)"lonlat:",lo,la
      do i=1,3
        write(*,*)"FORWARD:"
        lo=DBLE(lon)
        la=DBLE(lat)
        write(*,*)"Input [lon,lat,azi,dist]:",lo,la,azi,dist
        write(*,*)"Input: lo1 la1 lo2 la2:",lo,la,lo2,la2
        call vincenty_inverse(lo,la,lo2,la2)
        write(*,*)"Output: [azi,baz,dist]",lo,la,lo2
        azi=azi+DBLE(90.0)
        write(*,*)"-----------------------------------"
      enddo
      end
C End of example