TV Infomercials, Zip Codes and The Haversine Formula

TV Infomercials and Zip Codes Covered by TV Stations

I had read an article called "Points within a Given Distance on a Sphere" by Dr. Math which discussed finding the distance between two points on a sphere. In this article I decided to extend the ideas in that article to a stored procedure that I could use to calculate the distance between two television stations. The transmitters of Broadcast TV stations typically have a broadcast radius of about 200 miles and so when assigning 800 numbers to TV commercials that are close you need to see if the are at least 200 miles apart. In the TV Direct response business it is also necessary to determine which TV stations have broadcast radii that contain a given zip code. This is one piece of information that is used in Sourcing Orders or determining which airdate generated a given order from a customer based upon that customer's zip code and other information. Order Sourcing tells us how many orders we received from a given airing of an infomercial on a given TV station.

By calculating the list of TV stations that contain a given zip code of a customer's order, we simplify the process of determining which airdate might have generated that customer's order. I may write another article that presents a total programming solution to Sourcing Orders but in this article I want to focus on the problem of calculating the distance between two points given their latitude and longitude which requires some special programming tricks and which also has widespread application to many problems and situations in addition to Order Sourcing for TV infomercials.

Calculating Distance Given Latitude and Longitude

There are dozens of examples of code in C++ and other languages showing how to calculate the distance between two points given their latitude and longitude using the "Law of Cosines" for spherical trigonometry. For example, all of these code examples calculate the distance using the formula below which produces an error in the value for the distance on computers:

// Computationally Inaccurate Distance Calculation on Computers
a = sin(lat1) * sin(lat2) 
b = cos(lat1) * cos(lat2) * cos(lon2 - lon1) 
c = arccos(a + b) 
d = R * c 

This code is mathematically correct but "Computationally Inaccurate" for calculating distances because the Inverse Cosine function on computers is truncated to "significant digits" which introduces a significant error in the distance calculation.

Programmers must re-write the Law of Cosines using various trigonometric identities to yield an expression that can be accurately calculated by computers without truncating so-called non-significant digits. One solution is to use the Haversine Functions, which, incidentally, were not named after "Joe Haversine." Haversine is the name of one of the several lesser-known trignometric functions that were used by navigators. The VERSINE (or versed sine) of an angle A is 1-cos(A). The HAVERSINE is half of the versine, or (1-cos(A))/2. And with minimum effort you can show that:

hav(A) = (1-cos(A))/2 = sin^2(A/2) // Where ^2 means squared

The function below is "Computationally Accurate" for computers, and is referred to as the Haversine Formula because of its use of sin(x/2)^2, which is half of the "versed sine" of x, or (1 - cos(x)) / 2. And the function, atan2, on computers does not round off values to significant digits, so we get an accurate value for the distance.

// lat1, lon1 double Latitude / Longitude of point 1 
(decimal degrees) 
// lat2, lon2 double Latitude / Longitude of point 2 (decimal degrees) 
// units char S-Statute Miles; N-Nautical Miles; K-Kilometers 

#include <math.h> 
#define pi 3.14159265358979323846 

double DistLatLong(double lat1, double lon1, double lat2, double lon2, char units) { 
   
   double dlon, dlat; 
   dlon = lon2 - lon1; 
   dlat = lat2 - lat1; 
   a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2; 
   c = 2 * atan2( sqrt(a), sqrt(1-a) ); 
   
   // R (Earth Radius) = 3956.0 mi = 3437.7 nm = 6367.0 km

   switch(units) 
   { 
   case 'S': // STATUTE MILES
      R = 3956.0; 
      break; 
   case 'N': // NAUTICAL 
      R = 3437.7; 
      break; 
   case 'K': // KILOMETERS 
      R = 6367.0; 
      break; 
   } 

   return (R * c); 
}

And for those people who prefer using SQL, you can create an SQL version of the Haversine Formula as follows:

CREATE FUNCTION DistLatLong (@lat1 float, @lon1 float, @lat2 float, @lon2 float)
RETURNS float
AS
BEGIN 
   -- Parameters in RADIANS, result in statute miles
   DECLARE @dlon float, @dlat float, @a float, @c float 
   SET @dlon = @lon2 - @lon1 SET @dlat = @lat2 - @lat1 
   SET @a = POWER(SIN(@dlat / 2.0), 2.0) + 
        COS(@lat1) * COS(@lat2) * POWER(SIN(@dlon / 2.0), 2.0) 
   SET @c = 2.0 * ATN2(SQRT(@a), SQRT(1.0 - @a)) 
   RETURN 3956.0 * @c
END
GO
CREATE FUNCTION dms2rad (@deg int, @min int, @sec float)
RETURNS float
AS
BEGIN 
   RETURN RADIANS(CONVERT(float, @deg) + CONVERT(float, @min)/60.0 + @sec/3600.0)
END

You might find this page useful.

If you have any questions about this article or software for airing infomercials, please feel free to contact me. Please remember that formula only gives a rough estimate of distance on a sphere. 

Bill SerGio, The Infomercial King™