anchor.diff:  c "c" is the scale factor, in km/radian
   
         IMPLICIT none
  -      REAL lon,lat,x,y,c,pi,r2d,d2r
  -      PARAMETER (pi=3.14159265,r2d=180/pi,d2r=pi/180)
  +      REAL lon,lat,x,y,c,pi,r2d
  +      PARAMETER (pi=3.14159265,r2d=180.0/pi)
   
         lon = x/c
         lat = 2*atan(exp(y/c))-pi/2
  
positive.code: C (Double precision)
  C----
        IMPLICIT none
        REAL*8 pi,r2d,d2r
        PARAMETER (pi=4.0d0*datan(1.0d0),r2d=1.8d2/pi,d2r=pi/1.8d2)
        REAL*8 lon1,lat1,lon2,lat2,colat1,colat2,dlon,dlat,a
        REAL*8 dist,az
  C----
  C Check if points are polar opposite
  C----
        if (lat1.eq.-lat2.and.mod(lon1-lon2,1.8d2).eq.0) then
            dist = pi
            az = 0.0d0
            goto 11
  
negative.code: C----
        a = sin(dlat/2.0e0)*sin(dlat/2.0e0) +
       1   cos(lat1*d2r)*cos(lat2*d2r)*sin(dlon/2.0e0)*sin(dlon/2.0e0)
        if (a.ge.1.0) then
            dist = 0.0e0
        else
            dist = 2.0e0*atan2(sqrt(a),sqrt(1.0-a))
  