Spherical excess of a triangle

Question

Spherical excess of a triangle

As we all know, the sum of angles of any planar triangle is equal to 180 degrees.

However, for a spherical triangle, the sum of angles is always greater than 180 degrees. The difference between the sum of the spherical triangle angles and 180 degrees is called spherical excess . The task is to compute the spherical excess of a triangle with given vertex coordinates.

Some background

A spherical triangle is a part of the sphere defined by three great circles of the sphere.

Both sides and angles of spherical triangle are measured in the term of angle measure, because each side can be considered as a intersection of the sphere and some planar angle with vertex at the center of the sphere:

Each three distinct great circles define 8 triangles, but we only take proper triangles into consideration, ie. triangles whose angle and side measures satisfy

It's convenient to define vertices of a triangle in terms of geographic coordinate system. To compute the length of an arc of sphere given the longitude λ and latitude Φ of its ends we can use formula:

, where

or more explicitely:

(source: https://en.wikipedia.org/wiki/Haversine_formula)

The two basic formulas that can be used to solve a spherical triangle are:

• the law of cosines:

• the law of sines:

(source: https://en.wikipedia.org/wiki/Spherical_trigonometry#Cosine_rules_and_sine_rules)

Given three sides, it's easy to compute the angles using the cosine rule:

Finally, the spherical excess of a triangle is defined:

What's interesting about the relation between the spherical excess of a triangle and its area:

So on a unit sphere, the excess of a triangle is equal to area of that triangle!

The task

Write a function or a program that will compute spherical excess of a triangle in degrees given the triangle vertices coordinates. The vertex coordinates are provided in terms of geographic coordinate system.

Each vertex should be passed in form [latitude in degrees][N|S][longitude in degrees][E|W]. Longitude and E or W can be skipped when latitude is 90 ie. 90N, 90S, 10N100E, 30S20W are proper vertex descriptions, while 80N or 55S are not.

The latitudes and longitudes are always integer in the test cases.

The answers with error less than one degree will be accepted (as in examples below). The result can be rendered as both real or integer therefore, up to your convenience.

Examples

Input

90N0E
0N0E
0N90E

Output

89.999989

Input

90N
0N0E
0N90E

Output

89.999989

Input

0N0E
0N179E
90N0E

Output

178.998863

Input

10N10E
70N20W  
70N40E

Output

11.969793

In all test cases longitude and latitude are integer numbers. Parsing the vertex coordinates is the part of the task, so a vertex must be passed as single string/literal , it's not allowed to pass 80N20E as four parameters/strings: 80, N, 20, E.

This is guaranteed that vertices are all distinct and neither two of the three vertices make an antipodal points pair.

Scoring

This is code-golf, so the shortest code wins.

Answer 1

Matlab, 288 266 bytes

Here the commented version that should explain what is going on:

                                  %parsing the input
for k=1:3;
    s=input('','s');              %request input
    if sum(s>57)<2;               %if we have only one letter, add arbitrary second coordinate
        s=[s,'0E'];
    end;
    S=1-2*(s(s>57)>80);           %calculate the sign of the coordinates
    s(s>57)=44;                   %replace letters with comma
    L(k,:)=eval(['[',s,']']).*S;  %evaluates string as list and multiply with signs
end;
i=[2,3,1];
                                  %calculate the angular distance between each pair of points
a=arrayfun(@distance,L(:,1),L(:,2),L(i,1),L(i,2))*pi/180;
                                  %evaluate the spherical excess
f=@(a,b,c)sum(acos((cos(a)-cos(b).*cos(c))./(sin(b).*sin(c))))-pi;
disp(f(a,a(i),a([3,1,2]))*180/pi)

Fully golfed (linebreaks can be removed):

for k=1:3;s=input('','s');if sum(s>57)<2;s=[s,'0E'];end;
s(s>57)=44;L(k,:)=eval([91,s,93]).*(1-2*(s(s<48)>80));end;
i=[2,3,1];p=pi/180;a=arrayfun(@distance,L(:,1),L(:,2),L(i,1),L(i,2))*p;
b=a(i);disp((sum(acos((cos(a([3,1,2]))-cos(b).*cos(a))./(sin(b).*sin(a))))-pi)/p)

Answer 2

Ruby, Rev 3 264 255 bytes

Major changes:

New constant r = 180/PI defined and used throughout function. e had to be initialized to +PI, so excess now counts downwards and is negated before returning.

t[] eliminated: Ruby allows the data that was assigned to t[] to be directly assigned to u,v,w.

Single i loop to do the job of two, ?: ternary operator switches between tasks.

Many other minor changes.

include Math
->s{r=180/e=PI
x=y=z=n=[]
9.times{|i|i<6?(u,v,w=eval(?[+s[i%3].gsub(/[NE]/,"/r,").gsub(/[SW]/,"/-r,")+"0]")
i%2<1&&x=y=z=1
n[i/2]=(z*=sin(u))+(y*=cos(v)*w=cos(u))+x*=w*sin(v)):e-=acos((n[i-7]-(c=n[i-6])*d=n[i-8])/sqrt((1-c*c)*(1-d*d)))}
-e*r}

Ruby, Rev 1 283 277 bytes

Requires an array of 3 strings.

include Math 
->s{x=y=z=n=[]
6.times{|i|t=eval(?[+s[i%3].gsub(/[NE]/,k="*PI/180,").gsub(/[SW]/,"*-1"+k)+"0]")
i%2<1&&x=y=z=1
n[i/2]=(z*=sin(u=t[0]))+(y*=cos(u)*cos(v=t[1]))+(x*=cos(u)*sin(v))}
e=-PI
3.times{|i|e+=acos((n[i-1]-n[i]*d=n[i-2])/sqrt((1-n[i]**2)*(1-d**2)))}
e/PI*180}

Overview

The lengths of the triangle sides on the unit sphere are equal to the angles between the vectors describing the two points. But we don't need to know that angle. It suffices to know the cosine of the angle, which is easily obtained from cartesian coordinates using the Dot Product.

Explanation

The input strings are converted into a string representation of an array, which is then evaluated and stored in t, as below. The final zero is not needed if two coordinates are given. If only latitude 90 is given, the zero is interpreted as the longitude.

Example:  70N20W --> [70*PI/180,20*-1*PI/180,0]

The Dot Products are of the form a.b=ax*bx+ay*by+az*bz. As the vectors are all of unit length, the dot product is equal to the cosine of the angle between the vectors.

In order to calculate them a loop is iterated 6 times passing twice through the input data. On even iterations 0,2,4 the variables x,y,z are set to 1 to start a new calculation. On each iteration these variables are multiplied by the x,y and z components of each vector, using the longitude and latitude data stored in t[0],t[1] (which for golfing reasons is also assigned to u,v). The sum of the variables is written to array n (the garbage values on even iterations being overwritten by the correct values on odd iterations) so that at the end n contains the 3 dot products [a.b, c.a, b.c].

For the cosine rule, we need the cosines of the three included angles between vertices, but we also need the sines. These are obtained as sqrt(1-cosine**2). As the sines are multiplied together the expression can be rearranged so that only one call to sqrt is required. The fact that we do not know if the sine was positive or negative does not matter, as the haversine formula always gives the positive sine anyway. The important physical quantity is the distance between the points, which is absolute and therefore always positive.

For each iteration i=0..2, we calculate the value for the angle opposite array element i-1 using the other elements i and i-2. Negative array subscripts like this are legal in Ruby, they just wrap around to the beginning of the array.

Ungolfed in test program

Requires three sets of coordinates on the same line, with spaces between them.

include Math
g=->s{
  n=[]         #array for dot products
  x=y=z=1      #it's required to use these variables once before the loop, for some bizarre reason
  6.times{|i|
    t=eval(?[+s[i%3].gsub(/[NE]/,k="*PI/180,").gsub(/[SW]/,"*-1"+k)+"0]")
    i%2<1&&x=y=z=1
    n[i/2]=(z*=sin(u=t[0]))+(y*=cos(u)*cos(v=t[1]))+(x*=cos(u)*sin(v))
  }

  e=-PI        #set e to -PI and begin accumulating angles
  3.times{|i|
    e+=acos((n[i-1]-n[i]*n[i-2])/sqrt((1-n[i]**2)*(1-n[i-2]**2)))
  }

e/PI*180}      #return value

puts g[gets.split]