%=====================================================================
%  Compute the surface area of a spherical triangle given the
%     great circle route side lengths. All units in km or radians. 
%=====================================================================

arclength_a = 80;
arclength_b = 100;
arclength_c = 130;

R = 6378;


%------ Compute planar surface area

p = ( arclength_a + arclength_b + arclength_c ) / 2;

PlanarArea = sqrt( p * (p-arclength_a) * (p-arclength_b) * (p-arclength_c) );


%------ Convert surface arc length to interior spherical angles

a = arclength_a / R;
b = arclength_b / R;
c = arclength_c / R;

%------ Compute spherical surface area using interior angles

s = (a + b + c) / 2;

tanE14 = sqrt( tan(s/2) * tan((s-a)/2) * tan((s-b)/2) * tan((s-c)/2) );

E1 = 4 * atan( tanE14 );

SphericalArea1 = R * R * E1;


%------ Compute surface spherical angles between triangle sides

A1 = acos( ( cos(a) - cos(b) * cos(c) ) / ( sin(b) * sin(c) ) );
B1 = acos( ( cos(b) - cos(c) * cos(a) ) / ( sin(c) * sin(a) ) );
C1 = acos( ( cos(c) - cos(a) * cos(b) ) / ( sin(a) * sin(b) ) );

%------ Compute spherical surface area using surface spherical angles

E2 = A1 + B1 + C1 - pi;

SphericalArea2 = R * R * E2;

%------ Display results

disp(PlanarArea);
disp(SphericalArea1);
disp(SphericalArea2);