[back home]
Octavio Miramontes' codex page
The Physics Institute, UNAM
Center for Complexity Sciences (C3)
Mexico City

All code here is given as such. They are related to my research activities and are shared as copyleft. Feel free to use them, however I do not provide any kind of support and by using them, the user accepts all responsability. All programs have been compiled using GNU-Pascal.



Haversine formula

This is the Haversine formula code. It gives an accurate estimation of the distance between two points on the Earth when the coordinates are given as signed decimals. I use this code to calculate distances registered as successive logs of coordinates coming from a Holux M-241 GPS in a project on Human Mobility. Notice that the distances are sensitive to the value of the Earth radius. The radius is not constant along the latitude. The equatorial radius may be used but a better estimation is obtained when a more appropriate value is given for the correct latitude. Please see here for more info. Versions  of this formula in other programming languages can be found here. A paper on human mobility using this formula can be found here.

Program Haversine;
(*
 *  This program calculates the distance between two points on
 *  Earth. Its uses the Haversine formula and gives results in
 *  meters (this is easy to change however). Pascal has not the
 *  atan2 function implemented so it was to be defined explicitly. Here
 *  we have used the one provided by Paul Schlyter
 *  http://stjarnhimlen.se/comp/tutorial.html#Pcode
 *
 *  More info: http://en.wikipedia.org/wiki/Haversine_formula
  *  *)

Const
          pi      = 3.14159265358979323846;
          half_pi = 1.57079632679489661923;
          rad     = 3.14159265358979323846/180;{conversion factor for degrees
                                                into radians}
  
Function atan2(y,x: extended): extended;
    begin
      if x = 0.0 then
        begin
          if y = 0.0 then
            (* Error! Give error message and stop program *)
          else if y > 0.0 then
            atan2 := half_pi
          else
            atan2 := -half_pi
        end
      else
        begin
          if x > 0.0 then
              atan2 := arctan(y/x)
          else if x < 0.0 then
            begin
              if y >= 0.0 then
                atan2 := arctan(y/x)+pi
              else
                atan2 := arctan(y/x)-pi
            end;
        end;
    end; {atan2}

Var
   radius,lon1,lon2,lat1,lat2,dlon,dlat,a,distance :extended;

Begin {of the main program}

   radius:=6378.137; {Earth equatorial radius in Km; as used in most GPS}

   {Input here coordinates of the two points:}

   lon1:=-99.133333; {The Zocalo square, Mexico City}
   lat1:=19.432778;
  
   lon2:=-99.1862;   {The National Museum of Anthropology, Mexico City}
   lat2:=19.4260;

 {The Haversine formula}

   dlon:= (lon2-lon1)*rad;
   dlat:= (lat2-lat1)*rad;
   a:= sqr(sin(dlat/2)) + cos(lat1*rad)*cos(lat2*rad)*sqr(sin(dlon/2));
   distance:= radius*(2*atan2(sqrt(a),sqrt(1-a)));
   writeln;
   writeln("Distance: ");
   writeln (round(distance*1000), "   Meters");
   writeln (distance:4:3, "  Km");
   End.