Archive for the ‘Java’ Category

Converting UTM to Latitude and Longitude (Or Vice Versa)

Steven Dutch, Natural and Applied Sciences, University of Wisconsin – Green Bay
First-time Visitors: Please visit Site Map and Disclaimer. Use “Back” to return here

Information on the UTM system
New! Javascript-Convert Between Geographic and UTM Coordinates
Spreadsheet For UTM Conversion
Help! My Data Doesn’t Look Like A UTM Grid!

I get enough inquiries on this subject that I decided to create a page for it.

Caution! Unlike latitude and longitude, there is no physical frame of reference for UTM grids. Latitude is determined by the earth’s polar axis. Longitude is determined by the earth’s rotation. If you can see the stars and have a sextant and a good clock set to Greenwich time, you can find your latitude and longitude. But there is no way to determine your UTM coordinates except by calculation.

UTM grids, on the other hand, are created by laying a square grid on the earth. This means that different maps will have different grids depending on the datum used (model of the shape of the earth). I saw US military maps of Germany shift their UTM grids by about 300 meters when a more modern datum was used for the maps. Also, old World War II era maps of Europe apparently used a single grid for all of Europe and grids in some areas are wildly tilted with respect to latitude and longitude.

The two basic references for converting UTM and geographic coordinates are U.S. Geological Survey Professional Paper 1395 and U. S. Army Technical Manual TM 5-241-8 (complete citations below). Each has advantages and disadvantages.

For converting latitude and longitude to UTM, the Army publication is better. Its notation is more consistent and the formulas more clearly laid out, making code easier to debug. In defense of the USGS, their notation is constrained by space, and the need to be consistent with cartographic literature and descriptions of several dozen other map projections in the book.

For converting UTM to latitude and longitude, however, the Army publication has formulas that involve latitude (the quantity to be found) and which require reverse interpolation from tables. Here the USGS publication is the only game in town.

Some extremely tiny terms that will not seriously affect meter-scale accuracy have been omitted.

Converting Between Decimal Degrees, Degrees, Minutes and Seconds, and Radians

(dd + mm/60 +ss/3600) to Decimal degrees (dd.ff)

dd = whole degrees, mm = minutes, ss = seconds

dd.ff = dd + mm/60 + ss/3600

Example: 30 degrees 15 minutes 22 seconds = 30 + 15/60 + 22/3600 = 30.2561

Decimal degrees (dd.ff) to (dd + mm/60 +ss/3600)

For the reverse conversion, we want to convert dd.ff to dd mm ss. Here ff = the fractional part of a decimal degree.

mm = 60*ff

ss = 60*(fractional part of mm)

Use only the whole number part of mm in the final result.

30.2561 degrees = 30 degrees

.2561*60 = 15.366 minutes

.366 minutes = 22 seconds, so the final result is 30 degrees 15 minutes 22 seconds

Decimal degrees (dd.ff) to Radians

Radians = (dd.ff)*pi/180

Radians to Decimal degrees (dd.ff)

(dd.ff) = Radians*180/pi

Degrees, Minutes and Seconds to Distance

A degree of longitude at the equator is 111.2 kilometers. A minute is 1853 meters. A second is 30.9 meters. For other latitudes multiply by cos(lat). Distances for degrees, minutes and seconds in latitude are very similar and differ very slightly with latitude. (Before satellites, observing those differences was a principal method for determining the exact shape of the earth.)

Converting Latitude and Longitude to UTM

Okay, take a deep breath. This will get very complicated, but the math, although tedious, is only algebra and trigonometry. It would sure be nice if someone wrote a spreadsheet to do this. Or better yet, a Javascript page.

P = point under consideration
F = foot of perpendicular from P to the central meridian. The latitude of F is called the footprint latitude.
O = origin (on equator)
OZ = central meridian
LP = parallel of latitude of P
ZP = meridian of P
OL = k0S = meridional arc from equator
LF = ordinate of curvature
OF = N = grid northing
FP = E = grid distance from central meridian
GN = grid north
C = convergence of meridians = angle between true and grid north

Another thing you need to know is the datum being used:

Datum Equatorial Radius, meters (a) Polar Radius, meters (b) Flattening (a-b)/a Use
NAD83/WGS84 6,378,137 6,356,752.3142 1/298.257223563 Global
GRS 80 6,378,137 6,356,752.3141 1/298.257222101 US
WGS72 6,378,135 6,356,750.5 1/298.26 NASA, DOD
Australian 1965 6,378,160 6,356,774.7 1/298.25 Australia
Krasovsky 1940 6,378,245 6,356,863.0 1/298.3 Soviet Union
International (1924) -Hayford (1909) 6,378,388 6,356,911.9 1/297 Global except as listed
Clake 1880 6,378,249.1 6,356,514.9 1/293.46 France, Africa
Clarke 1866 6,378,206.4 6,356,583.8 1/294.98 North America
Airy 1830 6,377,563.4 6,356,256.9 1/299.32 Great Britain
Bessel 1841 6,377,397.2 6,356,079.0 1/299.15 Central Europe, Chile, Indonesia
Everest 1830 6,377,276.3 6,356,075.4 1/300.80 South Asia

Don’t interpret the chart to mean there is radical disagreement about the shape of the earth. The earth isn’t perfectly round, it isn’t even a perfect ellipsoid, and slightly different shapes work better for some regions than for the earth as a whole. The top three are based on worldwide data and are truly global. Also, you are very unlikely to find UTM grids based on any of the earlier projections.

The most modern datums (jars my Latinist sensibilities since the plural of datum in Latin is data, but that has a different meaning to us) are NAD83 and WGS84. These are based in turn on GRS80. Differences between the three systems derive mostly from redetermination of station locations rather than differences in the datum. Unless you are locating a first-order station to sub-millimeter accuracy (in which case you are way beyond the scope of this page) you can probably regard them as essentially identical.

I have no information on the NAD83 and WGS84 datums, nor can my spreadsheet calculate differences between those datums and WGS84.

Formulas For Converting Latitude and Longitude to UTM

These formulas are slightly modified from Army (1973). They are accurate to within less than a meter within a given grid zone. The original formulas include a now obsolete term that can be handled more simply – it merely converts radians to seconds of arc. That term is omitted here but discussed below.


  • lat = latitude of point
  • long = longitude of point
  • long0 = central meridian of zone
  • k0  = scale along long0 = 0.9996. Even though it’s a constant, we retain it as a separate symbol to keep the numerical coefficients simpler, also to allow for systems that might use a different Mercator projection.
  • e = SQRT(1-b2/a2) = .08 approximately. This is the eccentricity of the earth’s elliptical cross-section.
  • e’2 = (ea/b)2 = e2/(1-e2) = .007 approximately. The quantity e’ only occurs in even powers so it need only be calculated as e’2.
  • n = (a-b)/(a+b)
  • rho = a(1-e2)/(1-e2sin2(lat))3/2. This is the radius of curvature of the earth in the meridian plane.
  • nu = a/(1-e2sin2(lat))1/2. This is the radius of curvature of the earth perpendicular to the meridian plane. It is also the distance from the point in question to the polar axis, measured perpendicular to the earth’s surface.
  • p = (long-long0) in radians (This differs from the treatment in the Army reference)

Calculate the Meridional Arc

S is the meridional arc through the point in question (the distance along the earth’s surface from the equator). All angles are in radians.

  • S = A’lat – B’sin(2lat) + C’sin(4lat) – D’sin(6lat) + E’sin(8lat), where lat is in radians and
  • A’ = a[1 – n + (5/4)(n2 – n3) + (81/64)(n4 – n5) …]
  • B’ = (3 tanS/2)[1 – n + (7/8)(n2 – n3) + (55/64)(n4 – n5) …]
  • C’ = (15 tan2S/16)[1 – n + (3/4)(n2 – n3)  …]
  • D’ = (35 tan3S/48)[1 – n + (11/16)(n2 – n3)  …]
  • E’ = (315 tan4S/512)[1 – n  …]

The USGS gives this form, which may be more appealing to some. (They use M where the Army uses S)

  • M = a[(1 – e2/4 – 3e4/64 – 5e6/256 ….)lat
    – (3e2/8 + 3e4/32 + 45e6/1024…)sin(2lat)
    + (15e4/256 + 45e6/1024 + ….)sin(4lat)
    – (35e6/3072 + ….) sin(6lat) + ….)] where lat is in radians

This is the hard part. Calculating the arc length of an ellipse involves functions called elliptic integrals, which don’t reduce to neat closed formulas. So they have to be represented as series.

Converting Latitude and Longitude to UTM

All angles are in radians.

y = northing = K1 + K2p2 + K3p4, where

  • K1 = Sk0,
  • K2 = k0 nu sin(lat)cos(lat)/2 = k0 nu sin(2 lat)/4
  • K3 = [k0 nu sin(lat)cos3(lat)/24][(5 – tan2(lat) + 9e’2cos2(lat) + 4e’4cos4(lat)]

x = easting = K4p + K5p3, where

  • K4 = k0 nu cos(lat)
  • K5 = (k0 nu cos3(lat)/6)[1 – tan2(lat) + e’2cos2(lat)]

Easting x is relative to the central meridian. For conventional UTM easting add 500,000 meters to x.

What the Formulas Mean

The hard part, allowing for the oblateness of the Earth, is taken care of in calculating S (or M). So K1 is simply the arc length along the central meridian of the zone corrected by the scale factor. Remember, the scale is a hair less than 1 in the middle of the zone, and a hair more on the outside.

All the higher K terms involve nu, the local radius of curvature (roughly equal to the radius of the earth or roughly 6,400,000 m), trig functions, and powers of e’2 ( = .007 ). So basically they are never much larger than nu. Actually the maximum value of K2 is about nu/4 (1,600,000), K3 is about nu/24 (267,000) and K5 is about nu/6 (1,070,000). Expanding the expressions will show that the tangent terms don’t affect anything.

If we were just to stop with the K2 term in the northing, we’d have a quadratic in p. In other words, we’d approximate the parallel of latitude as a parabola. The real curve is more complex. It will be more like a hyperbola equatorward of about 45 degrees and an ellipse poleward, at least within the narrow confines of a UTM zone. (At any given latitude we’re cutting the cone of latitude vectors with an inclined plane, so the resulting intersection will be a conic section. Since the projection cylinder has a curvature, the exact curve is not a conic but the difference across a six-degree UTM zone is pretty small.)  Hence the need for higher order terms. Now p will never be more than +/-3 degrees = .05 radians, so p2 is always less than .0025 (1/400) and p4 is always less than .00000625 (1/160000). Using a spreadsheet, it’s easy to see how the individual terms vary with latitude. K2p2 never exceeds 4400 and K3p4 is at most a bit over 3. That is, the curvature of a parallel of latitude across a UTM zone is at most a little less than 4.5 km and the maximum departure from a parabola is at most a few meters.

K4 is what we’d calculate for easting in a simple-minded way, just by calculating arc distance along the parallel of latutude. But, as we get farther from the central meridian, the meridians curve inward, so our actual easting will be less than K4. That’s what K5 does. Since p is never more than +/-3 degrees = .05 radians, p3 is always less than .000125 (1/8000). The maximum value of K5p3 is about 150 meters.

That Weird Sin 1″ Term in the Original Army Reference

The Army reference defines p in seconds of arc and includes a sin 1″ term in the K formulas. The Sin 1″ term is a holdover from the days when this all had to be done on mechanical desk calculators (pre-computer) and terms had to be kept in a range that would  retain sufficient precision at intermediate steps. For that small an angle the difference between sin 1″ and 1″ in radians is negligible. If p is in seconds of arc, then (psin 1″) merely converts it to radians.

The sin 1″ term actually included an extra factor of 10,000, which was then corrected by multiplying by large powers of ten afterward.

The logic is a bit baffling. If I were doing this on a desk calculator, I’d factor out as many terms as possible rather than recalculate them for each term. But perhaps in practice the algebraically obvious way created overflows or underflows, since calculators could only handle limited ranges.

In any case, the sin1″ term is not needed any more. Calculate p in radians and omit the sin1″ terms and the large power of ten multipliers.

Converting UTM to Latitude and Longitude

In response to innumerable e-mails, you cannot use UTM grid coordinates without knowing your zone. There are sixty points on the earth’s surface that have the same numerical UTM coordinates, 120 if you consider that northing is duplicated in both hemispheres.

y = northing, x = easting (relative to central meridian; subtract 500,000 from conventional UTM coordinate).

Calculate the Meridional Arc

This is easy: M = y/k0.

Calculate Footprint Latitude

  • mu = M/[a(1 – e2/4 – 3e4/64 – 5e6/256…)
  • e1 = [1 – (1 – e2)1/2]/[1 + (1 – e2)1/2]

footprint latitude fp = mu + J1sin(2mu) + J2sin(4mu) + J3sin(6mu) + J4sin(8mu), where:

  • J1 = (3e1/2 – 27e13/32 ..)
  • J2 = (21e12/16 – 55e14/32 ..)
  • J3 = (151e13/96 ..)
  • J4 = (1097e14/512 ..)

Calculate Latitude and Longitude

  • e’2 = (ea/b)2 = e2/(1-e2)
  • C1 = e’2cos2(fp)
  • T1 = tan2(fp)
  • R1 = a(1-e2)/(1-e2sin2(fp))3/2. This is the same as rho in the forward conversion formulas above, but calculated for fp instead of lat.
  • N1 = a/(1-e2sin2(fp))1/2. This is the same as nu in the forward conversion formulas above, but calculated for fp instead of lat.
  • D = x/(N1k0)

lat = fp – Q1(Q2 – Q3 + Q4), where:

  • Q1 = N1 tan(fp)/R1
  • Q2 = (D2/2)
  • Q3 = (5 + 3T1 + 10C1 – 4C12 -9e’2)D4/24
  • Q4 = (61 + 90T1 + 298C1 +45T12  – 3C12 -252e’2)D6/720

long = long0 + (Q5 – Q6 + Q7)/cos(fp), where:

  • Q5 = D
  • Q6 = (1 + 2T1 + C1)D3/6
  • Q7 = (5 – 2C1 + 28T1 – 3C12 + 8e’2 + 24T12)D5/120

What Do The Formulas Mean?

As the sketch above shows, because of the poleward curve of parallels, the footprint latitude is always greater than the true latitude. Q1 is just a scaling coefficient and is constant for any given fp. The tangent term basically means the closer to the pole you are, the faster the parallels curve. Q2 is a quadratic term in x. Again, as with converting from geographic coordinates to UTM, we approximate the parallel as a parabola and add higher order corrections.

To determine longitude, we make a simple minded approximation that longitude is proportional to easting, but then, since fp is too large, the true longitude is smaller, since it lies on a parallel closer to the the equator. The divisor cos(fp) corrects for the varying length of degrees of longitude as latitude varies.

A Spreadsheet Program

Before linking to the program, note (especially the last item):

  • It is an Excel spreadsheet. If you have Excel on your computer, it will (or should) open when you click the link. Most major spreadsheet programs can read spreadsheets in other formats.
  • A spreadsheet is not an applet or program. In particular, you can’t manually enter data into a cell and preserve any formulas that are there. That’s why some aspects of data entry are clunkier than they might otherwise be with, say, a Visual Basic program.
  • There are three computation pages, one for single conversions, the other two for batch conversions. Other pages contain information on datums and the specific conversion formulas. To use the batch conversions you need to be somewhat proficient in spreadsheets as you will have to copy data and cell formulas.
  • For our mutual peace of mind, run a virus check.
  • You may copy the program for your own non-commercial use and for non-commercial distribution to others, but not for commercial use. Please give appropriate credit when citing your calculations. You may also modify it as needed for your personal use. In Internet Explorer, right click on the link and select Save Target As… to save the spreadsheet to your own disk.
  • There is now a Javascript Page which I hope is free of some of the bugs in the spreadsheet. Since it’s a program, it should handle complex what-ifs a lot better than a spreadsheet can. All the permissions for the spreadsheet apply here, too.

Spreadsheet For UTM Conversion
JavaScript Page to Convert Between Geographic and UTM Coordinates


Snyder, J. P., 1987; Map Projections – A Working Manual. U.S. Geological Survey Professional Paper 1395, 383 p. If you are at all serious about maps you need this book.

Army, Department of, 1973; Universal Transverse Mercator Grid, U. S. Army Technical Manual TM 5-241-8, 64 p.

NIMA Technical Report 8350.2, “Department of Defense World Geodetic System 1984, Its Definition and Relationships with Local Geodetic Systems,” Second Edition, 1 September 1991 and its supplements. The report is available from the NIMA Combat Support Center and its stock number is DMATR83502WGS84. Non-DoD requesters may obtain the report as a public sale item from the U.S. Geological Survey, Box 25286, Denver Federal Center, Denver, Colorado 80225 or by phone at 1-800-USA-MAPS.



/*---------------------------------------------------------*/ /* parameters and functions for UTM<->Lat/Long conversions */ /* UTM conversion is based on WGS84 ellipsoid parameters */ /*---------------------------------------------------------*/

RADIUS = 6378137.0;
FLATTENING = 0.00335281068; /* GRS80 or WGS84 */
K_NOT = 0.9996;     /* UTM scale factor */
DEGREES_TO_RADIANS = 0.01745329252;
FALSE_EASTING = 500000.0;
FALSE_NORTHING = 10000000.0;

/* These Java functions were converted from original      */
/* C fucntions provided by Dana Yoerger, Steve Gegg,      */
/* and Louis Whitcomb at WHOI                             */

function geo_utm(lat, lon, zone)
      /* first compute the necessary geodetic parameters and constants */

      lambda_not = ((-180.0 + zone*6.0) -3.0)/RADIANS ;
      e_fourth = e_squared * e_squared;
      e_sixth = e_fourth * e_squared;
      e_prime_sq = e_squared/(1.0 - e_squared);
      sin_phi = sin(lat);
      tan_phi = tan(lat);
      cos_phi = cos(lat);
      N = RADIUS/sqrt(1.0 - e_squared*sin_phi*sin_phi);
      T = tan_phi*tan_phi;
      C = e_prime_sq*cos_phi*cos_phi;
      M = RADIUS*((1.0 - e_squared*0.25 -0.046875*e_fourth  -0.01953125*e_sixth)*lat-
	      (0.375*e_squared + 0.09375*e_fourth +
				 0.043945313*e_sixth)*sin(2.0*lat) +
	      (0.05859375*e_fourth + 0.043945313*e_sixth)*sin(4.0*lat) -
	      (0.011393229 * e_sixth)*sin(6.0*lat));
      A = (lon - lambda_not)*cos_phi;
      A_sq = A*A;
      A_fourth =  A_sq*A_sq;

      /* now go ahead and compute X and Y */

      x_utm = K_NOT*N*(A + (1.0 - T + C)*A_sq*A/6.0 +
		   (5.0 - 18.0*T + T*T + 72.0*C -

      /* note:  specific to UTM, vice general trasverse mercator.
         since the origin is at the equator, M0, the M at phi_0,
         always equals zero, and I won't compute it   */                                            

       y_utm = K_NOT*(M + N*tan_phi*(A_sq/2.0 +
			    (5.0 - T + 9.0*C + 4.0*C*C)*A_fourth/24.0 +
			    (61.0 - 58.0*T + T*T + 600.0*C -

       /* now correct for false easting and northing */

       if( lat < 0)
          y_utm +=10000000.0;
       x_utm +=500000;
    return true;

function utm_geo(x_utm, y_utm, zone)
      /* first, subtract the false easting */
      x_utm = x_utm - FALSE_EASTING;

      /* compute the necessary geodetic parameters and constants */

      e_fourth = e_squared * e_squared;
      e_sixth = e_fourth * e_squared;
      oneminuse = sqrt(1.0-e_squared);  

      /* compute the footpoint latitude */

      M = y_utm/K_NOT;
      mu =M/(RADIUS*(1.0 - 0.25*e_squared -
                  0.046875*e_fourth - 0.01953125*e_sixth));
      e1 = (1.0 - oneminuse)/(1.0 + oneminuse);
      e1sq =e1*e1;
      footpoint = mu + (1.5*e1 - 0.84375*e1sq*e1)*sin(2.0*mu) +
              (1.3125*e1sq - 1.71875*e1sq*e1sq)*sin(4.0*mu) +
              (1.57291666667*e1sq*e1)*sin(6.0*mu) +

      /* compute the other necessary terms */

      e_prime_sq =  e_squared/(1.0 -  e_squared);
      sin_phi = sin(footpoint);
      tan_phi = tan(footpoint);
      cos_phi = cos(footpoint);
      N = RADIUS/sqrt(1.0 - e_squared*sin_phi*sin_phi);
      T = tan_phi*tan_phi;
      Tsquared = T*T;
      C = e_prime_sq*cos_phi*cos_phi;
      Csquared = C*C;
      denom = sqrt(1.0-e_squared*sin_phi*sin_phi);
      R = RADIUS*oneminuse*oneminuse/(denom*denom*denom);
      D = x_utm/(N*K_NOT);
      Dsquared = D*D;
      Dfourth = Dsquared*Dsquared;

      lambda_not = ((-180.0 + zone*6.0) -3.0) * DEGREES_TO_RADIANS;

      /* now, use the footpoint to compute the real latitude and longitude */

      var lat = footpoint - (N*tan_phi/R)*(0.5*Dsquared - (5.0 + 3.0*T + 10.0*C -
                           4.0*Csquared - 9.0*e_prime_sq)*Dfourth/24.0 +
                           (61.0 + 90.0*T + 298.0*C + 45.0*Tsquared -
                            252.0*e_prime_sq -
      var lon = lambda_not + (D - (1.0 + 2.0*T + C)*Dsquared*D/6.0 +
                         (5.0 - 2.0*C + 28.0*T - 3.0*Csquared + 8.0*e_prime_sq +

   return true;


Calculators and Convertors

Posted: September 3, 2011 in Java, Mapping

Degrees, Minutes, Seconds and Decimal Degrees Latitude/Longitude Conversions

Circle Intersection Solver

Online Math Calculators and Solvers

Calculate distance, bearing and more between two latitude/longitude points

NDSF Coordinate Conversion Utility

Convert Between Geographic and UTM Coordinates

using a Latitude and Longitude value (Point A), I am trying to calculate another Point B, X meters away bearing 0 radians from point A. Then display the point B Latitude and Longitude values.

Example (Pseudo code):

PointA_Lat = x.xxxx; PointA_Lng = x.xxxx; Distance = 3; //Meters bearing = 0; //radians new_PointB = PointA-Distance; 

I was able to calculate the distance between two Points but what I want to find is the second point knowing the distance and bearing.

Preferably in PHP or Javascript.

Thank you


4 down vote accepted

It seems you are measuring distance (R) in meters, and bearing (theta) counterclockwise from due east. And for your purposes (hundereds of meters), plane geometry should be accurate enough. In that case,

dx = R*cos(theta) ; theta measured counterclockwise from due east dy = R*sin(theta) ; dx, dy same units as R 

If theta is measured clockwise from due north (for example, compass bearings), the calculation for dx and dy is slightly different:

dx = R*sin(theta) ; theta measured clockwise from due north dy = R*cos(theta) ; dx, dy same units as R 

In either case, the change in degrees longitude and latitude is:

delta-longitude = dx/(111320*cos(latitude)) ; dx, dy in meters delta-latitude = dy/110540 ; result in degrees long/lat 

The difference between the constants 110540 and 111320 is due to the earth’s oblateness (polar and equatorial circumferences are different).

Here’s a worked example, using the parameters from a later question of yours:

Given a start location at longitude -87.62788 degrees, latitude 41.88592 degrees, find the coordinates of the point 500 meters northwest from the start location.

If we’re measuring angles counterclockwise from due east, “northwest” corresponds to theta=135 degrees. R is 500 meters.

dx = R*cos(theta) = 500 * cos(135 deg) = -353.55 meters dy = R*sin(theta) = 500 * sin(135 deg) = +353.55 meters delta_longitude = dx/(111320*cos(latitude)) = -353.55/(111320*cos(41.88592 deg)) = -.004266 deg (approx -15.36 arcsec) delta-latitude = dy/110540 = 353.55/110540 = .003198 deg (approx 11.51 arcsec) Final longitude = start_longitude + delta_longitude = -87.62788 - .004266 = -87.632146 Final latitude = start_latitude + delta_latitude = 41.88592 + .003198 = 41.889118 


Positional Astronomy:
Spherical trigonometry

A great-circle arc, on the sphere, is the analogue of a straight line, on the plane.

Where two such arcs intersect, we can define the spherical angle
either as angle between the tangents to the two arcs, at the point of intersection,
or as the angle between the planes of the two great circles where they intersect at the centre of the sphere.
(Spherical angle is only defined where arcs of great circles meet.)

A spherical triangle is made up of three arcs of great circles, all less than 180°.
The sum of the angles is not fixed, but will always be greater than 180°.
If any side of the triangle is exactly 90°, the triangle is called quadrantal.

There are many formulae relating the sides and angles of a spherical triangle.
In this course we use only two: the sine rule and the cosine rule.

diagramConsider a triangle ABC on the surface of a sphere with radius = 1.

(See note)

We use the capital letters A, B, C to denote the angles at these corners;
we use the lower-case letters a, b, c to denote the opposite sides.
(Remember that, in spherical geometry, the side of a triangle is the arc of a great circle,
so it is also an angle.)

Turn the sphere so that A is at the “north pole”,
and let arc AB define the “prime meridian”.

Set up a system of rectangular axes OXYZ:
O is at the centre of the sphere;
OZ passes through A;
OX passes through arc AB (or the extension of it);
OY is perpendicular to both.
Find the coordinates of C in this system:
x = sin(b) cos(A)
y = sin(b) sin(A)
z = cos(b)


Now create a new set of axes,
keeping the y-axis fixed and moving the “pole” from A to B
(i.e. rotating the x,y-plane through angle c).
The new coordinates of C are
x’ = sin(a) cos(180-B) = – sin(a) cos(B)
y’ = sin(a) sin(180-B) =    sin(a) sin(B)
z’ = cos(a)

The relation between the old and new systems
is simply a rotation of the x,z-axes through angle c:
x’ = x cos(c) – z sin(c)
y’ = y
z’ = x sin(c) + z cos(c)

That is:
– sin(a) cos(B) = sin(b) cos(A) cos(c) – cos(b) sin(c)
sin(a) sin(B)  = sin(b) sin(A)
cos(a) = sin(b) cos(A) sin(c) + cos(b) cos(c)

These three equations give us the formulae for solving spherical triangles.

The first equation is the transposed cosine rule,
which is sometimes useful but need not be memorised.

The second equation gives the sine rule. Rearrange as:
   sin(a)/sin(A) = sin(b)/sin(B)
sin(b)/sin(B) = sin(c)/sin(C), etc.

So the sine rule is usually expressed as:
   sin(a)/sin(A) = sin(b)/sin(B) = sin(c)/sin(C)

The third equation gives the cosine rule:
cos(a) = cos(b) cos(c) + sin(b) sin(c) cos(A)
and similarly:
cos(b) = cos(c) cos(a) + sin(c) sin(a) cos(B)
cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C)

Here they are together:

sine rule:
   sin(a)/sin(A) = sin(b)/sin(B) = sin(c)/sin(C)

cosine rule:
cos(a) = cos(b) cos(c) + sin(b) sin(c) cos(A)
cos(b) = cos(c) cos(a) + sin(c) sin(a) cos(B)
cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C)

The cosine rule will solve almost any triangle if it is applied often enough.
The sine rule is simpler to remember but not always applicable.

Note that both formulae can suffer from ambiguity:
E.g. if the sine rule yields
sin(x) = 0.5,
then x may be 30° or 150°.
Or, if the cosine rule yields
cos(x) = 0.5,
then x may be 60° or 300° (-60°).
In this case, there is no ambiguity if x is a side of the triangle, as it must be less than 180°,
but there could still be uncertainty if an angle of the triangle was positive or negative.

So, when applying either formula, check to see if the answer is sensible.
If in doubt, recalculate using the other formula, as a check.


Alderney, in the Channel Islands, has longitude 2°W, latitude 50°N.
Winnipeg, in Canada, has longitude 97°W, latitude 50°N.
How far apart are they, in nautical miles, along a great-circle arc?

If you set off from Alderney on a great-circle route to Winnipeg,
in what direction (towards what azimuth) would you head?

Click here for the answer.



Aviation Formulary V1.46

Posted: August 18, 2011 in Android, Java, Mapping

By Ed Williams


This introduction is written for pilots (and others) who are interested in great circle navigation and would like to know how to compute courses, headings and other quantities of interest. These formulae can be programmed into your calculator or spreadsheet. I’ll attempt to include enough information that those familiar with plane trigonometry can derive additional results if required.

It is a well known that the shortest distance between two points is a straight line. However anyone attempting to fly from Los Angeles to New York on the straight line connecting them would have to dig a very substantial tunnel first. The shortest distance, following the earth’s surface lies vertically above the aforementioned straight line route. This route can be constructed by slicing the earth in half with an imaginary plane through LAX and JFK. This plane cuts the (assumed spherical) earth in a circular arc connecting the two points, called a great circle. Only planes through the center of the earth give rise to great circles. Any plane will cut a sphere in a circle, but the resulting little circles are not the shortest distance between the points they connect. A little thought will show that lines of longitude (meridians) are great circles, but lines of latitude, with the exception of the equator, are not.

I will assume the reader is familiar with latitude and longitude as a means of designating locations on the earth’s surface. For the convenience of North Americans I will take North latitudes and West longitudes as positive and South and East negative. The longitude is the opposite of the usual mathematical convention. True course is defined as usual, as the angle between the course line and the local meridian measured clockwise.

The first important fact to realise is that in general a great circle route has a true course that varies from point to point. For instance the great circle route between two points of equal (non-zero) latitude does not follow the line of latitude in an E-W direction, but arcs towards the pole. It is possible to fly between two points using an unvarying true course, but in general the resulting route differs from the great circle route and is called a rhumb line. Unlike a great circle which encircles the earth in a closed curve, a pilot flying a rhumb line would spiral poleward.

Natural questions are to seek the great circle distance between two specified points and true course at points along the route. The required spherical trigonometric formulae are greatly simplified if angles and distances are measured in the appropriate natural units, which are both radians! A radian, by definition, is the angle subtended by a circular arc of unit length and unit radius. Since the length of a complete circular arc of unit radius is 2*pi, the conversion is 360 degrees equals 2*pi radians, or:


Great circle distance can be likewise be expressed in radians by defining the distance to be the angle subtended by the arc at the center of the earth. Since by definition, one nautical mile subtends one minute (=1/60 degree) of arc, we have:


In all subsequent formulae all distances and angles, such as latitudes, longitudes and true courses will be assumed to be given in radians, greatly simplifying them, and in applications the above formulae and their inverses are necessary to convert back and forth between natural and practical units. Examples of this process are given later.

Note: the nautical mile is currently defined to be 1852 meters – which to be consistent with its historical definition implies the earth’s radius to be 1.852 * (180*60/pi) = 6366.71 km, which indeed lies between the currently accepted ( WGS84) equatorial and polar radii of 6378.137 and 6356.752 km, respectively. Other choices of the earth’s radius in this range are consistent with the spherical approximation and may for some specialized purposes be preferred.

For instance, the former FAI standard for aviation records used the “FAI sphere” with radius 6371.0 km.

To use a different spherical radius, use

      distance_km = radius_km * distance_radians
      distance_radians = distance_km/radius_km

to convert between angular and distance units.

A sample implementation of many of these formulae in the form of an Excel spreadsheet can be found here. The formulae are in VBA macros, for readability, so you will need to enable macros to have them work. If you are unable to open the spreadsheet with macros disabled (to check for viruses) etc, then you may need to patch your Excel. Try for Excel 97 SR-2.

If you decide to program up these formulae, you’d be well-advised to look at the implementation notes and check your results against the worked examples and spreadsheets.

Some great circle formulae:

Distance between points

The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by:


A mathematically equivalent formula, which is less subject to rounding error for short distances is:

d=2*asin(sqrt((sin((lat1-lat2)/2))^2 +
Lat/lon given radial and distance

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

     IF (cos(lat)=0)
        lon=lon1      // endpoint a pole

This algorithm is limited to distances such that dlon <pi/2, i.e those that extend around less than one quarter of the circumference of the earth in longitude. A completely general, but more complicated algorithm is necessary if greater distances are allowed:

     lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
     lon=mod( lon1-dlon +pi,2*pi )-pi

Implementation notes:

Notes on mathematical functions

Note: ^ denotes the exponentiation operator, sqrt is the square root function, acos the arc-cosine (or inverse cosine) function and asin is the arc-sine function. If asin or acos are unavailable they can be implemented using the atan2 function:

  acos returns a value in the range 0 <= acos <= pi
  asin returns a value in the range -pi/2 <= asin <= pi/2

Note: Here atan2 has the conventional (C) ordering of arguments, namely atan2(y,x). This is not universal, Excel for instance uses atan2(x,y), but it has asin and acos anyway. Be warned. It returns a value in the range -pi < atan2 <= pi.

Further note: if your calculator/programming language is so impoverished that only atan is available then use:

   acos(x)=2*atan(sqrt((1-x)/(1+x)))       x>=0
          =pi - 2*atan(sqrt((1+x)/(1-x)))  x<0

   atan2(y,x)=atan(y/x)       x>0
   atan2(y,x)=atan(y/x)+pi    x<0, y>=0
   atan2(y,x)=pi/2            x=0, y>0
   atan2(y,x)=atan(y/x)-pi    x<0, y<0
   atan2(y,x)=-pi/2           x=0, y<0
   atan2(0,0) is undefined and should give an error.

Another potential implementation problem is that the arguments of asin and/or acos may, because of rounding error, exceed one in magnitude. With perfect arithmetic this can’t happen. You may need to use “safe” versions of asin and acos on the lines of:


Note on the mod function. This appears to be implemented differently in different languages, with differing conventions on whether the sign of the result follows the sign of the divisor or the dividend. (We want the sign to follow the divisor or be Euclidean. C’s fmod and Java’s % do not work.) In this document, Mod(y,x) is the remainder on dividing y by x and always lies in the range 0 <=mod <x. For instance: mod(2.3,2.)=0.3 and mod(-2.3,2.)=1.7

If you have a floor function (int in Excel), that returns floor(x)= “largest integer less than or equal to x” e.g. floor(-2.3)=-3 and floor(2.3) =2

        mod(y,x) = y - x*floor(y/x)

The following should work in the absence of a floor function- regardless of whether "int" truncates or rounds downward:

    mod=y - x * int(y/x)
    if ( mod < 0) mod = mod + x
Sign Convention

As stated in the introduction, North latitudes and West longitudes are treated as positive, and South latitudes and East longitudes negative. It’s easier to go with the flow, but if you prefer another convention you can change the signs in the formulae.


Vincenty formula for distance between two Latitude/Longitude points


For the benefit of the terminally obsessive (as well as the genuinely needy), Thaddeus Vincenty (‘TV’) devised formulae for calculating geodesic distances between a pair of latitude/longitude points on the earth’s surface, using an accurate ellipsoidal model of the earth.

When I looked up the references (‘Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations’), I discovered to my surprise that while the mathematics is utterly beyond me, it is actually quite simple to program.

Vincenty’s formula is accurate to within 0.5mm, or 0.000015″ (!), on the ellipsoid being used. Calculations based on a spherical model, such as the (much simpler) Haversine, are accurate to around 0.3% (which is still good enough for most purposes, of course).

Note: the accuracy quoted by Vincenty applies to the theoretical ellipsoid being used, which will differ (to varying degree) from the real earth geoid. If you happen to be located in Colorado, 2km above msl, distances will be 0.03% greater. In the UK, if you measure the distance from Land’s End to John O’ Groats using WGS-84, it will be 28m – 0.003% – less than using the Airy ellipsoid, which provides a better fit for the UK. (Since of course no one can travel on the surface of the theoretical ellipsoid, these differences are generally of no relevance at all).

Enter the co-ordinates into the text boxes to try it out (using deg-min-sec suffixed with N/S/E/W, or signed decimal degrees):

Lat 1: Long 1:

Lat 2: Long 2:

Note: nearly-antipodal points may fail to give a solution, in which case NaN is returned.

Vincenty wrote this at a time when computing resources were very expensive, and made it very computationally efficient. The JavaScript should be quite straightforward to translate into other languages, if required.

See the ‘Direct’ formula to calculate the destination point given the start point, bearing and distance.

Vincenty’s formula as it is used in the script:

a, b = major & minor semiaxes of the ellipsoid
f = flattening (ab)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’)
U2 = atan((1−f).tanφ2)
λ = L (first approximation)
iterate until change in λ is negligible (e.g. 10-12 ≈ 0.06mm) {
sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] (14)
cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ (15)
σ = atan2(sinσ, cosσ) (16)
sinα = cosU1.cosU2.sinλ / sinσ (17)
cos²α = 1 − sin²α (trig identity; §6)
cos2σm = cosσ − 2.sinU1.sinU2/cos²α (18)
C = f/16.cos²α.[4+f.(4−3.cos²α)] (10)
λ′ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]} (11)
u² = cos²α.(a²−b²)/b²
A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]} (3)
B = u²/1024.{256+u².[−128+u².(74−47.u²)]} (4)
Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]} (6)
s = b.A.(σ−Δσ) (19)
α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) (20)
α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ) (21)


  • s is the distance (in the same units as a & b)
  • α1 is the initial bearing, or forward azimuth
  • α2 is the final bearing (in direction p1→p2)

Note: Vincenty observes that eqn. (18) becomes indeterminate over equatorial lines (since cos²α → 0); in this case, set cos2σm to 0 and the result is computed correctly. He also points out that the formula may have no solution between two nearly antipodal points; an iteration limit traps this case (Vincenty says “this will occur when λ, as computed by eqn. (11), is greater than π in absolute value”, but this is not always a reliable test).
Note: some implementations of Vincenty’s formula inefficiently use a large number of trig functions; Vincenty devised this solution with an eye for efficiency in implementation, and this one uses just one each of sin, cos, sqrt, and atan2 for each iteration – only 3 or 4 iterations are generally required. [Formulation updated Dec 05 to make it closer to Vincenty’s original and computationally more efficient.]

The most accurate and widely used globally-applicable model for the earth ellipsoid is WGS-84, used in this script. Other ellipsoids offering a better fit to the local geoid include Airy (1830) in the UK, International 1924 in much of Europe, Clarke (1880) in Africa, GRS-67 in South America, and many others. America (NAD83) and Australia (GDA) use GRS-80, functionally equivalent to the WGS-84 ellipsoid.

WGS-84 a = 6 378 137 m (±2 m) b ≈ 6 356 752.314245 m f ≈ 1 / 298.257223563
GRS-80 a = 6 378 137 m b ≈ 6 356 752.314140 m f = 1 / 298.257222101
Airy 1830 a = 6 377 563.396 m b = 6 356 256.910 m f ≈ 1 / 299.3249646
Internat’l 1924 a = 6 378 388 m b ≈ 6 356 911.946 m f = 1 / 297
Clarke mod.1880 a = 6 378 249.145 m b ≈ 6 356 514.86955 m f = 1 / 293.465
GRS-67 a = 6 378 160 m b ≈ 6 356 774.719 m f = 1 / 298.247167

Note that to locate latitude/longitude points on these ellipses, they are associated with specific datums: for instance, OSGB36 for Airy in the UK, ED50 for Int’l 1924 in Europe; WGS-84 defines a datum as well as an ellipse. See my convert coordinates page for converting points between different datums.



Aside for fellow geeks: it is often said that WGS-84 and GRS-80 are ‘functionally equivalent’. Well, I would concur that a 0.1 mm difference in Earth radius is pretty equivalent, but why any difference? Apparently, when defining WGS-84, the normalized second degree zonal harmonic gravitational coefficient C2,0, derived from the GRS-80 value for the dynamic form factor J2, was truncated to 8 significant digits in the normalization process!

Some of the terms involved are explained in Ed Williams’ notes on Spheroid Geometry.

Test case (from Geoscience Australia), using WGS-84:
Flinders Peak 37°57′03.72030″S, 144°25′29.52440″E
Buninyong 37°39′10.15610″S, 143°55′35.38390″E
s 54 972.271 m
α1 306°52′05.37″
α2 127°10′25.07″ (≡ 307°10′25.07″ p1→p2)


  • Trig functions take arguments in radians, so latitude, longitude, and bearings in degrees (either decimal or degrees/minutes/seconds) need to be converted to radians, rad = π.deg/180. When converting radians back to degrees (deg = 180.rad/π), West is negative if using signed decimal degrees. For bearings, values in the range -π to +π [-180° to +180°] need to be converted to 0 to +2π [0°–360°]; this can be done by (brng+2.π)%2.π [brng+360)%360] where % is the modulo operator.
  • The atan2() function used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).
  • If you implement any formula involving atan2 in Microsoft Excel, you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conventional order is atan2(y, x), but Excel uses atan2(x, y). To use atan2 in a (VBA) macro, you can use WorksheetFunction.Atan2(). You will also have to rename the variables A, B / a, b, as VBA is case-insensitive.
  • All bearings are with respect to true north, 0°=N, 90°=E, etc; if you are working from a compass, magnetic north varies from true north in a complex way around the earth, and the difference has to be compensated for by variances indicated on local maps.
  • For miles, divide km by 1.609344
  • For nautical miles, divide km by 1.852

See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required.

Creative Commons License I offer these formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.

Paypal donation If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.

© 2002-2010 Chris Veness

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2010             */
/*                                                                                                */
/* from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
/*       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975    */
/*                                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

 * Calculates geodetic distance between two points specified by latitude/longitude using 
 * Vincenty inverse formula for ellipsoids
 * @param   {Number} lat1, lon1: first point in decimal degrees
 * @param   {Number} lat2, lon2: second point in decimal degrees
 * @returns (Number} distance in metres between points
function distVincenty(lat1, lon1, lat2, lon2) {
  var a = 6378137, b = 6356752.314245,  f = 1/298.257223563;  // WGS-84 ellipsoid params
  var L = (lon2-lon1).toRad();
  var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
  var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
  var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
  var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
  var lambda = L, lambdaP, iterLimit = 100;
  do {
    var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
    var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
      (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    if (sinSigma==0) return 0;  // co-incident points
    var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    var sigma = Math.atan2(sinSigma, cosSigma);
    var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    var cosSqAlpha = 1 - sinAlpha*sinAlpha;
    var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    if (isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
    var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    lambdaP = lambda;
    lambda = L + (1-C) * f * sinAlpha *
      (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
  } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

  if (iterLimit==0) return NaN  // formula failed to converge

  var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
  var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
  var s = b*A*(sigma-deltaSigma);
  s = s.toFixed(3); // round to 1mm precision
  return s;
  // note: to return initial/final bearings in addition to distance, use something like:
  var fwdAz = Math.atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
  var revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
  return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() };

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */


<p>Lat 1: <input name="lat1" value="53 09 02N"> Long 1: <input name="long1" value="001 50 40W"></p>
<p>Lat 2: <input name="lat2" value="52 12 19N"> Long 2: <input name="long2" value="000 08 33W"></p>

<input type="button" value="calculate distance" 
    onClick="result.value = distVincenty(Geo.parseDMS(lat1.value), Geo.parseDMS(long1.value), 
                                         Geo.parseDMS(lat2.value), Geo.parseDMS(long2.value)) + ' m'">
<input name="result" value="">

As we know, Java 1.5 has not been maintained anymore and Java 6, has been hanging around for a while and Java 7 is coming soon. But it doesn’t mean everybody has to move on to Java 1.6. The problem in Ubuntu is they force you to use OpenJDK. Worse yet, they don’t let you downgrade to 1.5.There are lots of legacy systems running on Java 5 and we can’t forget that. It is also about freedom of choice. If you want to install Sun JDK 5 or 6 Ubuntu should not make your life difficult.

Ubuntu 10 and 11 don’t allow us to natively (or easily) install Sun JDK’s via apt-get.

It is so frustrate when we try to install it on Ubuntu 10.10 and have no luck:

sudo apt-get install sun-java5-jdk (or sudo apt-get install sun-java6-jdk)

After installed, the JDK will go to this directory: /usr/lib/jvm/java-6-sun

You can see all JDK’s installed running the following command:

sudo update-java-alternatives -l

After a while, finally got the solution for this:

sudo add-apt-repository “deb jaunty multiverse”

sudo add-apt-repository “deb jaunty-updates multiverse”

sudo apt-get update

sudo apt-get install sun-java5-jdk

Check it out just to confirm:

sudo update-java-alternatives -l

The commands above worked for me in Ubuntu 10 and 11, installing Sun JDK 1.5 and 1.6.
If it doesn’t work, try to add the following repositories and repeat the process: (anonymous suggestion – thanks a lot):

sudo add-apt-repository “deb hardy multiverse”
sudo add-apt-repository “deb hardy-updates multiverse”


Deriving the Haversine Formula

Date: 04/20/99 at 14:27:24
From: Dena Neff
Subject: Longitude and Latitude Calculations

I'm an SAS programmer at The Coca-Cola Company. It's been a long, long 
time since I have had to use trigonometry. I need to write a program 
module to calculate distances given longitude and latitude data. I am 
trying to find an object within a mile's radius of its location. Would 
you have the equation for this scenario?

Thank you.  I hope I am not too old for a response.  :)

Date: 04/20/99 at 17:16:15
From: Doctor Rick
Subject: Re: Longitude and Latitude Calculations

Hello, Dena, welcome to Ask Dr. Math. We're allowed to answer "old" 
people now and then, if the question is interesting enough! :-)

You can find several expositions of the "Law of Cosines" for spherical
trigonometry in our Dr. Math Archives. For example:

   Using Longitude and Latitude to Determine Distance 

   Distance using Latitude and Longitude   

But a correspondent recently directed our attention to a Web site,

  The Geographic Information Systems FAQ (See Q5.1)   

This site points out that the cosine formula does not give accurate 
results for small distances, and presents a better formula, the 
Haversine Formula. The form of the Haversine Formula that uses the 
atan2() function is probably best for programming, assuming your 
programming language has an atan2() function. I will copy some 
excerpts here for your convenience, but I recommend looking at the 
site for the other information you will find there.

Note that in the description that follows, atan2(y, x) = arctan(y/x). 
However, there does seem to be some confusion about the order of the 
arguments. C/C++ puts the arguments in the order that I assume, as 
indicated by the Unix man page as well as the Visual C++ Programmer's 
Guide. On the other hand, Excel describes the function (not too 
clearly) as "ATAN2(x_num, y_num) returns the arctangent of the 
specified x- and y-coordinates, in radians." I tested it, and indeed it 
takes the denominator x first.


Presuming a spherical Earth with radius R (see below), and that the
locations of the two points in spherical coordinates (longitude and
latitude) are lon1,lat1 and lon2,lat2, then the Haversine Formula 
(from R. W. Sinnott, "Virtues of the Haversine," Sky and Telescope, 
vol. 68, no. 2, 1984, p. 159): 

  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)) 
  d = R * c

will give mathematically and computationally exact results. The 
intermediate result c is the great circle distance in radians. The 
great circle distance d will be in the same units as R.

Most computers require the arguments of trigonometric functions to be
expressed in radians. To convert lon1,lat1 and lon2,lat2 from degrees, 
minutes, and seconds to radians, first convert them to decimal 
degrees. To convert decimal degrees to radians, multiply the number of 
degrees by pi/180 = 0.017453293 radians/degree. 

Inverse trigonometric functions return results expressed in radians. 
To express c in decimal degrees, multiply the number of radians by 
180/pi = 57.295780 degrees/radian. (But be sure to multiply the number 
of RADIANS by R to get d.) 

The historical definition of a "nautical mile" is "one minute of arc 
of a great circle of the earth." Since the earth is not a perfect 
sphere, that definition is ambiguous. However, the internationally 
accepted (SI) value for the length of a nautical mile is (exactly, by 
definition) 1.852 km or exactly 1.852/1.609344 international miles 
(that is, approximately 1.15078 miles - either "international" or 
"U.S. statute"). Thus, the implied "official" circumference is 360 
degrees times 60 minutes/degree times 1.852 km/minute = 40003.2 km. 
The implied radius is the circumference divided by 2 pi: 

R = 6367 km = 3956 mi 

I'll just add that the intermediate result "a" is the square of half 
of the straight-line distance (chord length) between the two points.

  a = (AB/2)^2
    = (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)

If you are interested in a derivation of the formula, I could provide 
it, but the figure would be hard to draw with "character graphics."

- Doctors Rick and Peterson, The Math Forum

Date: 04/21/99 at 09:14:55
From: Dena Neff
Subject: Longitude and Latitude Calculations

Doctor Rick,  

Thank you very much for your quick response. I'm looking forward to 
visiting the Web sites you provided to review the expositions 
of the "Law of Cosines" and the Haversine Formula.

I think I will always have a nature of curiosity about me. Therefore, 
I would LOVE to see a derivation of the formula.


Date: 04/21/99 at 22:42:45
From: Doctor Rick
Subject: Re: Longitude and Latitude Calculations

Hello again, Dena.

You called my bluff! Good - now I have a reason to write up my 
derivation of the haversine formula. I had derived the cosine formula 
for great-circle distance, but not the haversine formula. I was 
reluctant to send the formula to you without having checked it out, so 
I derived it before I answered you. It turned out to be relatively 

First, a word about the meaning of "haversine." This is the name of 
one of several lesser-known trigonometric functions that were once 
familiar to navigators and the like. The VERSINE (or versed sine) of 
angle A is 1-cos(A). The HAVERSINE is half the versine, or 
(1-cos(A))/2. With a little math, you can show that

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

where ^2 means squared. You see the form on the right twice in the
haversine formula.

In proving the formula, we will work with a sphere of radius 1. The
distance we get should be multiplied by R, the radius of the earth.

Let's start with a formula for the length of a chord subtending an 
angle theta in the unit circle (circle of radius 1).

      / | \
   1 /  |  \ 1
    /   |   \
   /    |    \
 A      C      B

O is the center of the circle, and A and B are on the circle, so AB is
the chord. If angle AOB = theta, then angle AOC = theta/2 where OC is
perpendicular to the chord AB. C bisects AB, since triangle AOB is
isosceles. Segment AC = sin(theta/2), so the length of chord AB is

Now for the 3-dimensional figure. Here is a sphere of radius 1, with a
wedge removed to show the center, labeled O. The points we are
interested in are A (lat1,lon1) and B (lat2,lon2). The lines of
longitude lon1 and lon2 are shown as curved single lines, meeting at 
the north pole, N; the equator is shown as a double line (equal 

                        ****     N     ****
                    ****     /   |  \      ****
                 ***       /     |    \        ***
               **        /       |      \         **
             **         /        |       \          **
            *          /         |        \           * 
           *          /          |         \           *
          *          A----------------------D           *
         *          /|\          |           \           *
         *         | |  \        |            |           *
        *          | |   \       |            |           *
       *          |  |    \      |             |            *
       *          |  |     \     |             |            *
      *           C----------\-----------------B            *
      *           |  |   -   _\  |             |            *
      *           |  |          _O_            |            * 
      * =        |   |     _          _         |         = *
       *  ===    |   G_                    _    |      === *
       *     ====E______________________________F =====    *
        *         =========              =========        *
         *                 ==============                *
         *                                               *
          *                                             * 
           *                                           *
            **                                       **
              **                                   **
                ***                             ***
                   ****                     ****
                       ****             ****

A and B are opposite vertices of an isosceles trapezoid, ACBD, with
additional vertices C (lat2,lon1) and D (lat1,lon2). I haven't tried 
to draw the straight lines (chords) connecting A to C and B to D, but 
I show the chords connecting A to D and B to C. I won't bother to 
prove that the trapezoid ACBD is planar.

If you join A to O and C to O, the angle AOC is the difference between
lat1 and lat2, or dlat. The length of chord AC is therefore
2*sin(dlat/2), using the chord-length formula derived above. Chord BD 
is the same length.

Points E and F are the points where the latitude lines lat1 and lat2
meet the equator. The length of EF is 2*sin(dlon/2) by the chord-
length formula.

Points A and D are on a circle of constant latitude lat1. The radius 
of this circle is cos(lat1). You can see this by dropping a 
perpendicular from A to OE, meeting OE at G. The angle EOA is lat1, 
so OG = cos(lat1). But this is equal to the radius of the latitude 
circle at A.

Therefore the length of chord AD is 2*sin(dlon/2)*cos(lat1). 
Similarly, the length of chord CB is 2*sin(dlon/2)*cos(lat2).

Now let's go back to 2 dimensions, and find the length of a diagonal 
of an isosceles trapezoid.

      /| *              \
     / |     *           \
    /  |         *        \
   /   |             *     \
 C     H                     B

The length of CH, where AH is perpendicular to CB, is (CB-AD)/2. By 
the Pythagorean theorem,

  AH^2 = AC^2 - CH^2
       = AC^2 - (CB-AD)^2/4

The length of HB is (CB+AD)/2. Using Pythagoras again, we find that

  AB^2 = AH^2 + HB^2
       = AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4
       = AC^2 + CB*AD

We're ready to plug in the lengths of chords AC, AD, and CB from our
work with the sphere:

  AB^2 = 4*(sin^2(dlat/2) + cos(lat1)*cos(lat2)*sin^2(dlon/2))

The intermediate result a is the square of half the chord length AB:

  a = (AB/2)^2
    = (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)

The final step is to find the central angle AOB that corresponds to 
this chord length. The arctan formula can be found from this figure:

            / | \
           /  |  \
          /   |   \
       1 /    |    \ 1
        /     |     \
       /      |      \
      /       |       \
     /        |        \
    /         |         \
   / sqrt(1-a)|          \
  /           |           \
A   sqrt(a)   C             B

If AC = sqrt(a), then Pythagoras says that

  OC = sqrt(OA^2 - AC^2)
     = sqrt(1-a)

Therefore tan(<AOC) = AC/OC = sqrt(a)/sqrt(1-a), or

  c = 2 * arctan(sqrt(a)/sqrt(1-a))

where c is the angle AOB. 



This uses the ‘haversine’ formula to calculate the great-circle distance between two points – that is, the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points (ignoring any hills, of course!).

Haversine formula: a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
where R is earth’s radius (mean radius = 6,371km);
note that angles need to be in radians to pass to trig functions!
var R = 6371; // km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
var lat1 = lat1.toRad();
var lat2 = lat2.toRad();

var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
var d = R * c;

The haversine formula1 ‘remains particularly well-conditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. The ‘versed sine’ is 1-cosθ, and the ‘half-versed-sine’ (1-cosθ)/2 = sin²(θ/2) as used above. It was published by R W Sinnott in Sky and Telescope, 1984, though known about for much longer by navigators. (For the curious, c is the angular distance in radians, and a is the square of half the chord length between the points). A (surprisingly marginal) performance improvement can be obtained, of course, by factoring out the terms which get squared.



Intersection of two circles

Posted: July 3, 2011 in Java

Intersection of two circles

The following note describes how to find the intersection point(s) between two circles on a plane, the following notation is used. The aim is to find the two points P3 = (x3, y3) if they exist.

First calculate the distance d between the center of the circles. d = ||P1 – P0||.

  • If d > r0 + r1then there are no solutions, the circles are separate.
  • If d < |r0 – r1| then there are no solutions because one circle is contained within the other.
  • If d = 0 and r0 = r1then the circles are coincident and there are an infinite number of solutions.

Considering the two triangles P0P2P3 and P1P2P3we can write

a2 + h2 = r02 and b2 + h2 = r12Using d = a + b we can solve for a,

a = (r02 – r12 + d2) / (2 d)

It can be readily shown that this reduces to r0 when the two circles touch at one point, ie: d = r0 + r1

Solve for h by substituting a into the first equation, h2 = r02 – a2


P2 = P0 + a ( P1 – P0 ) / dAnd finally, P3 = (x3,y3) in terms of P0 = (x0,y0), P1 = (x1,y1) and P2 = (x2,y2), is

x3 = x2 +- h ( y1 – y0) / d

y3 = y2 -+ h ( x1 – x0 ) / d