//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [sang] = CL_gm_inters3dConesSa(alpha,alpha1,alpha2)
// Solid angle of the intersection of two cones
// 
// Calling Sequence
// sang = CL_gm_inters3dConesSa(alpha,alpha1,alpha2)
//
// Description
// <itemizedlist>
// <listitem>
// <p>Computes the solid angle of the intersection of two cones.</p>
// <p>The cones are defined by their half angles <b>alpha1</b> and <b>alpha2</b>, and
// the angular separation between their axes <b>alpha</b>.</p>
// <p><inlinemediaobject><imageobject><imagedata fileref="inters_3d_cones_sa.png"/></imageobject></inlinemediaobject></p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p>- The half angles of the cones should be in [0,pi/2].</p>
// <p>  If an half-angle equals to pi/2 then the cone actually is an half-space.</p>
// <p>- The angular separation between axes of the cones should be in [0,pi].</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// alpha: Angular separation between axes of the cones (1xN or 1x1)
// alpha1: Half angle of first cone (1xN or 1x1)
// alpha2: Half angle of second cone (1xN or 1x1)
// sang: Solid angle of the intersection of the cones (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Intersection of two cones
// alpha = 0.15;
// alpha1 = 0.1;
// alpha2 = 0.2;
// sang = CL_gm_inters3dConesSa(alpha,alpha1,alpha2);

  I = find(alpha1 < 0 | alpha2 < 0 | alpha < 0 | ...
           alpha1 > %pi/2 | alpha2 > %pi/2 | alpha > %pi);
  if (I ~= []) ; CL__error("Check domain of validity for angles"); end;
  
  // Resize inputs arguments if necessary
  [alpha,alpha1,alpha2] = CL__checkInputs(alpha,1, alpha1,1, alpha2,1);
  
  sang = %nan * ones(alpha);
    
  // Sort alpha1, alpha2 so that alpha1 < alpha2
  alpha_1 = min(alpha1, alpha2);
  alpha2 = max(alpha1, alpha2);
  alpha1 = alpha_1;

  ca1 = cos(alpha1);
  ca2 = cos(alpha2);

  // Angular separation between axes of the cones such that smaller cone is included in bigger cone
  I = find(alpha <= alpha2 - alpha1);
  sang(I) = 2*%pi * (1 - ca1(I)); // intersection = smaller cone
    
  // Angular separation such that there is no intersection
  I = find(alpha >= alpha1 + alpha2);
  sang(I) = 0;
  
  // Standard case
  I = find(alpha < alpha1 + alpha2 & alpha > alpha2 - alpha1);
    
  if (I <> [])
    X   = (alpha1(I) + alpha2(I) + alpha(I)) / 2; 
    sina1 = sin(X - alpha1(I)); 
    sina2 = sin(X - alpha2(I)); 
    sina  = sin(X - alpha(I)); 
    R   = sqrt(sina1 .* sina2 .* sina ./ sin(X)); 

    beta1 = 2 * atan(R, sina2); 
    beta2 = 2 * atan(R, sina1); 
    betaa = 2 * atan(R, sina); 

    // Solid angle of a cone is 2*pi*(1-cos(alpha1)).
    // Solid angle of a "portion" of a cone is beta1*(1-cos(alpha1))
    // Solid angle of a triangle of sides alpha,alpha1,alpha2 is: betaa+beta1+beta2-pi
    // So the solid angle of intersection is:

    sang(I) = 2 * (beta1 .* (1 - ca1(I)) + beta2 .* (1 - ca2(I)) - (betaa + beta1 + beta2 - %pi)); 
  end

endfunction