  function [Cost,R,G,B,C,nFevals] = Tetra(r,g,b,c)
% r,g,b,c specifies the cost per unit area of red, green, blue, and cyan paint. 
% Cost returns the value of the most expensive tetrahedron that is circumscribed
% by the unit sphere and whose four faces are painted red, green, blue and
% cyan. The areas of the red, green, blue, and cyan faces are reported
% in R,G,B,and C. The number of times that fminsearch calls the objective
% function is returned in nFevals. 

% The tetrahedron has four vertices: P0,P1,P2,P3
% We fix P0 at the north pole and keep P1 on zero longitude meridian.
% Thius leaves 5 unknowns which are used as follows
%        P0: (pi/2,0)
%        P1: (x(1),0)          x(1) the latitude of P1
%        P2: (x(2),x(3))       x(2) and x(3) are the lat and long of P2
%        P3: (x(4),x(5))       x(4) and x(5) are the lat and long of P3

% Initial guess is approximately the Platonic tetrahedron
x0 = [-pi/10 -pi/10   2*pi/3  -pi/10   -2*pi/3]';
% Set fminsearch parameters...
options = optimset('MaxFunEvals',1000,'TolFun',.00000001,'TolX',.00000001);
% The optimization...
[x,Cost,EXITFLAG,OUTPUT] = fminsearch(@(x) TetraCost(x,r,g,b,c),x0,options);
% Output...
P0 = Sphere2Cart(pi/2,0);
P1 = Sphere2Cart(x(1),0);
P2 = Sphere2Cart(x(2),x(3));
P3 = Sphere2Cart(x(4),x(5));
R = TriArea(P0,P1,P2);
G = TriArea(P0,P2,P3);
B = TriArea(P0,P3,P1);
C = TriArea(P1,P2,P3);
Cost = -Cost;
nFevals = OUTPUT.funcCount;



  function P = Sphere2Cart(theta,phi)
% theta: latitude in radians (-pi/2<=theta<=pi/2)
% phi:   longitude in radians (-pi<=phi<=pi)
% (P.x,P.y,P.z) Cartesian coordinates of a point on the unit sphere
% that has latitude theta and longitude phi
P = struct('x',cos(theta)*cos(phi),'y',cos(theta)*sin(phi),'z',sin(theta));

  function A = TriArea(P1,P2,P3)
% Area of the triangle with vertices at points P1, P2, and P3
d12 = sqrt((P1.x-P2.x)^2+(P1.y-P2.y)^2+(P1.z-P2.z)^2);
d23 = sqrt((P2.x-P3.x)^2+(P2.y-P3.y)^2+(P2.z-P3.z)^2);
d31 = sqrt((P3.x-P1.x)^2+(P3.y-P1.y)^2+(P3.z-P1.z)^2);
s = (d12+d23+d31)/2;
A = sqrt(s*(s-d12)*(s-d23)*(s-d31));

  function c = TetraCost(x,r,g,b,c)
% x is a column 5-vector that specifies a tetrahedron
% c is the negative  of the tetrahedron cost
  P0 = Sphere2Cart(pi/2,0);
  P1 = Sphere2Cart(x(1),0);
  P2 = Sphere2Cart(x(2),x(3));
  P3 = Sphere2Cart(x(4),x(5));
  A012 = TriArea(P0,P1,P2);    % The red face
  A023 = TriArea(P0,P2,P3);    % The green face
  A031 = TriArea(P0,P3,P1);    % The blue face
  A123 = TriArea(P1,P2,P3);    % The cyan face
  c = -(r*A012 + g*A023 + b*A031 + c*A123);
  
  % up to -2 for badly chosen initial guess
  % -2 for 8 variables
  % -1 for 6 variables 
  
