anchor.diff:  !!  Starting guess for the iteration variables
   !   lambda1 = big_number
   !   lambda2 = L
  -!	
  +!   
   !   do while ( abs(lambda2-lambda1) > convergence_limit ) 
  -!	  sin_s = sqrt((cos(u2)*sin(lambda2))**2 + (cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(lambda2))**2)
  -!	  cos_s = sin(u1)*sin(u2) + cos(u1)*cos(u2)*cos(lambda2)
  -!	  s = atan2(sin_s,cos_s)
  -!	  sin_a = cos(u1)*cos(u2)*sin(lambda2)/sin_s
  -!	  cos_a_sq = 1.d0 - sin_a**2
  -!	  cos_2sm = cos_s - s*sin(u1)*sin(u2)/cos_a_sq
  +!     sin_s = sqrt((cos(u2)*sin(lambda2))**2 + (cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(lambda2))**2)
  +!     cos_s = sin(u1)*sin(u2) + cos(u1)*cos(u2)*cos(lambda2)
  +!     s = atan2(sin_s,cos_s)
  +!     sin_a = cos(u1)*cos(u2)*sin(lambda2)/sin_s
  +!     cos_a_sq = 1.d0 - sin_a**2
  +!     cos_2sm = cos_s - s*sin(u1)*sin(u2)/cos_a_sq
   !      if (isnan(cos_2sm)) cos_2sm = 0.d0
  -!	  C = (f/16.d0)*cos_a_sq*(4.d0+f*(4.d0-3.d0*cos_a_sq))
  -!	  lambda1 = lambda2
  -!	  lambda2 = L + (1.d0-C)*f*sin_a*(s+C*sin_a*(cos_2sm+C*cos_s*(-1.d0+2.d0*cos_2sm**2)))
  +!     C = (f/16.d0)*cos_a_sq*(4.d0+f*(4.d0-3.d0*cos_a_sq))
  +!     lambda1 = lambda2
  +!     lambda2 = L + (1.d0-C)*f*sin_a*(s+C*sin_a*(cos_2sm+C*cos_s*(-1.d0+2.d0*cos_2sm**2)))
   !   enddo
  -!	
  +!   
   !   u_2 = cos_a_sq*(a**2-b**2)/b**2
   !   big_A = 1.d0+(u_2/16384.d0)*(4096.d0+u_2*(-768.d0+u_2*(320.d0-175.d0*u_2)))
   !   big_B = (u_2/1024.d0)*(256.d0+u_2*(-128.d0+u_2*(74.d0-47.d0*u_2)))
   !   Delta_s = big_B*sin_s*(cos_2sm+(big_B/4.d0)*(cos_s*(-1.d0+2.d0*cos_2sm**2)- &
  -!			 (big_B/6.d0)*cos_2sm*(-3.d0+4.d0*sin_s**2)*(-3.d0+4.d0*cos_2sm**2)))
  +!          (big_B/6.d0)*cos_2sm*(-3.d0+4.d0*sin_s**2)*(-3.d0+4.d0*cos_2sm**2)))
   !   test_dist_vincenty = b*big_A*(s-Delta_s) ;
   !   
   !   
  
positive.code: !  delta returns the angular distance between two points on a sphere given 
  !  the lat and lon of each using the Haversine formula
  !
  	  implicit none
  	  real(rs) :: delta,lat1,lon1,lat2,lon2
  	  logical,optional :: degrees
  	  	  
  	  if (present(degrees)) then
  	     if (degrees) then
  			lat1=lat1*pi/1.8D2 ; lon1=lon1*pi/1.8D2
  			lat2=lat2*pi/1.8D2 ; lon2=lon2*pi/1.8D2
  	     endif
  	  endif
  	  
  	  delta=atan2( sqrt( (cos(lat2)*sin(lon2-lon1))**2 + (cos(lat1)*sin(lat2) - &
  			  sin(lat1)*cos(lat2)*cos(lon2-lon1))**2) , &
  			  sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1))
  
  	  if (present(degrees)) then
           if (degrees) delta = delta * 1.8D2/pi
        endif
  	  
  	  return
  	  
     end function delta
  !==============================================================================
  
  
negative.code:       real(rs),intent(out)  :: dist
        real(rs)             :: a,b,f,lon1,lat1,lon2,lat2,L,u1,u2,lambda1,lambda2,&
                                sin_s,cos_s,s,sin_a,cos_a_sq,cos_2sm,C,u_2,&
                                big_A,big_B,ds
        logical              :: isnan
  !      real,intent(out),optional :: azi,baz
        
  