Path: typhoon01.sbc.net!cyclone.swbell.net!cyclone-sf.pbi.net!63.208.208.143!feed2.onemain.com!feed1.onemain.com!news-spur1.maxwell.syr.edu!news.maxwell.syr.edu!newsfeed.skycache.com!Cidera!portc03.blue.aol.com!uunet!dca.uu.net!bbnews1.unisys.com!plnews.pl.unisys.com!not-for-mail
From: "Clive Tooth" <clive@pisquaredoversix.force9.co.uk>
Newsgroups: sci.math
Subject: Re: Area of Sphereical Triangles
Date: Mon, 23 Apr 2001 12:55:38 +0100
Organization: Correct Programmnig Inc
Lines: 144
Message-ID: <9c153j$i59$1@mail.pl.unisys.com>
References: <3AE3CAFB.C375E83E@mediaone.net>
NNTP-Posting-Host: 192.39.129.234
X-Newsreader: Microsoft Outlook Express 4.72.3612.1700
X-MimeOLE: Produced By Microsoft MimeOLE V4.72.3612.1700
Xref: cyclone.swbell.net sci.math:239828

Franck wrote in message <3AE3CAFB.C375E83E@mediaone.net>...

>   I was hoping someone here could help me. If I have 3 Lat Lon Points
>on a globe, where can I find the formulas to compute the area between
>them?

I will address the case where you are prepared to take the Earth to be
spherical. In this case the following Delphi code may be used (use the
function TriangleArea):

######################################################
type

  PointType=record
    px: floating;
    py: floating;
    pz: floating
  end; { PointType }

  LongLatType=record
    Long: floating; { East is positive }
    Lat: floating   { North is positive }
  end; { LongLatType }

function TriangleArea(LL1, LL2, LL3: LongLatType; var Error: boolean):
floating;
const
  EarthMeanRadiusKilometers=6371.0;
  RadiansPerDegree=Pi/180;
var
  a1, a2, a3: floating;
  p1, p2, p3: PointType;
  u1, u2, u3: PointType;

  function  ComputeDihedralAngles(u1, u2: PointType): floating;
forward;
  function  ComputeUnitNormals(p1, p2: PointType): PointType;
forward;
  function  ConvertLongLatToXYZ(const LL: LongLatType): PointType;
forward;
  function  CosD(x: floating): floating;
forward;
  function  SinD(x: floating): floating;
forward;

  function ComputeDihedralAngles(u1, u2: PointType): floating;
  begin
    Result:= ArcCos(-u1.px*u2.px-u1.py*u2.py-u1.pz*u2.pz)
  end; { ComputeDihedralAngles }

  function ComputeUnitNormals(p1, p2: PointType): PointType;
  var
    a: floating;
  begin
    with Result do begin
      px:= p1.py*p2.pz - p1.pz*p2.py;
      py:= p1.pz*p2.px - p1.px*p2.pz;
      pz:= p1.px*p2.py - p1.py*p2.px;
      a:= 1/Sqrt(Sqr(px)+Sqr(py)+Sqr(pz));
      px:= px*a;
      py:= py*a;
      pz:= pz*a;
    end
  end; { ComputeUnitNormals }

  function ConvertLongLatToXYZ(const LL: LongLatType): PointType;
  begin
    with LL, Result do begin
      px:= CosD(Lat)*CosD(Long);
      py:= CosD(Lat)*SinD(Long);
      pz:= SinD(Lat)
    end
  end; { ConvertLongLatToXYZ }

  function CosD(x: floating): floating;
  begin
    Result:= Cos(x*RadiansPerDegree)
  end; { CosD }

  function SinD(x: floating): floating;
  begin
    Result:= Sin(x*RadiansPerDegree)
  end; { SinD }

begin
  Result:= 0;
  Error:= true;
  try
    p1:= ConvertLongLatToXYZ(LL1);
    p2:= ConvertLongLatToXYZ(LL2);
    p3:= ConvertLongLatToXYZ(LL3);

    u1:= ComputeUnitNormals(p2, p3);
    u2:= ComputeUnitNormals(p3, p1);
    u3:= ComputeUnitNormals(p1, p2);

    a1:= ComputeDihedralAngles(u2, u3);
    a2:= ComputeDihedralAngles(u3, u1);
    a3:= ComputeDihedralAngles(u1, u2);

    Result:= (a1+a2+a3-Pi)*Sqr(EarthMeanRadiusKilometers);

    Error:= false
  except
    {}
  end
end; { TriangleArea }

function L(x, y: floating): LongLatType;
begin
  with Result do begin
    Long:= x;
    Lat:= y
  end
end; { L }
######################################################

The function will set Error to true if the three points are ill-conditioned
in some way (co-planar, coincident, etc).

For example to compute the area of a triangle bounded by
  Toronto (in Ontario, Canada) at 79.417 West, 43.667 North
  Tichcax (in Yucatan, Mexico) at 89.700 West, 20.867 North
  Saco do Tamanduá (in Paraiba, Brazil) at 37.467 West, 6.817 South
use the call:

TriangleArea(L(-79.417, 43.667), L(-89.700, 20.867), L(-37.467, -6.817), b)

The result turns out to be about 9,778,439.29 square kilometers.

See also:

For some FORTRAN routines (see REAL FUNCTION AREAS):
http://www.netlib.org/toms/772

For some thoughts on what the mean radius of the Earth is:
http://home.online.no/~sigurdhu/WGS84_Eng.html

--
Clive Tooth
http://www.pisquaredoversix.force9.co.uk/
End of document