function [coordinates, solidAngles, probabilities] = ...
    qtilesky(tiling, thetaRange, phiRange)
% QTILESKY Tile the sky with a specified angular resolution
%
% QTILESKY returns a set of spherical coordinates that covers the sky, or a
% specified region of the sky, with (close to) the minimum number of points
% required to ensure an angular resolution consistent with the frequency
% range of the search.
%
% usage:
%
%   [coordinates, solidAngles, probabilities] = ...
%       qtilesky(tiling, thetaRange, phiRange);
%
%   tiling           q transform tiling structure from QTILE
%   thetaRange       range of polar angles to tiles [radians]
%   phiRange         range of azimuthal angles to tile [radians]
%
%   coordinates      matrix of sky positions [radians]
%   solidAngles      vector of solid angles [steradians]
%   probabilities    vector of a priori uniform probabilities []
%
% The sky positions coordinates are provided as a two column matrix of spherical
% coordinates in the form [theta phi] as used by QSKYMAP.  The coordinate theta
% is a geocentric colatitude running from 0 at the North pole to pi at the South
% pole, and the coordinate phi is the geocentric longitude in Earth fixed
% coordinates with 0 on the prime meridian.  The range of theta is [0, pi] and
% the range of phi is [-pi, pi).
%
% QTILESKY also returns the approximate differential solid angle associated with
% each sky position, and its corresponding a priori uniform probability density.
%
% The optional theta and phi range parameters are two element vectors of the
% form [min max] that can be used to restrict the region of the sky that is
% tiled.  By default, the whole sky is tiled.
%
% The angular resolution of the sky tiling is determined from the specified q
% transform tiling structure in order to ensure a worst case phase mismatch of
% pi/4 for signals at the maximum frequency of the analysis.
%
% See also QTILE, QRESPONSE, and QSKYMAP.

%   Shourov Chatterji   shourov@ligo.caltech.edu
%   Albert Lazzarini    lazz@ligo.caltech.edu
%   Antony Searle       antony.searle@anu.edu.au
%   Leo Stein           lstein@ligo.caltech.edu
%   Patrick Sutton      psutton@ligo.caltech.edu
%   Massimo Tinto       massimo.tinto@jpl.nasa.gov

% $Id: qtilesky.m,v 1.4 2007/07/26 19:10:40 shourov Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        process command line arguments                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% verify correct number of input arguments
error(nargchk(1, 3, nargin));

% apply default arguments
if (nargin < 2) || isempty(thetaRange),
  thetaRange = [0 pi];
end
if (nargin < 3) || isempty(phiRange),
  phiRange = [-pi pi];
end

% force column vector ranges
thetaRange = thetaRange(:);
phiRange = phiRange(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       validate command line arguments                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% validate tiling structure
if ~strcmp(tiling.id, 'Discrete Q-transform tile structure'),
  error('input argument is not a discrete Q transform tiling structure');
end

% validate theta range
if (length(thetaRange) ~= 2) | (diff(thetaRange) <= 0),
  error('theta range must be two element vector [min max]');
end
if any((thetaRange < 0) | (thetaRange > pi)),
  error('theta range outside of [0, pi]');
end

% validate phi range
if (length(phiRange) ~= 2) | (diff(phiRange) <= 0),
  error('phi range must be two element vector [min max]');
end
if any((phiRange < -pi) | (phiRange > pi)),
  error('phi range outside of [-pi, pi]');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              angular resolution                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% worst case phase shift due to wrong sky position
phaseResolution = pi / 4;

% required timing resolution at maximum frequency of search
timingResolution = phaseResolution / (2 * pi * tiling.lowPassCutoff);

% maximum possible time delay between earth based detectors
earthRadius = 6378135;
speedOfLight = 299792458;
maximumTimeDelay = 2 * earthRadius / speedOfLight;

% requred angular resolution
angularResolution = pi * timingResolution / maximumTimeDelay;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         construct coordinate matrix                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize result vectors
thetas = [];
phis = [];
solidAngles = [];

% determine number of thetas
numberOfThetas = ceil(diff(thetaRange) / angularResolution);

% determine theta spacing
thetaStep = diff(thetaRange) / numberOfThetas;

% begin loop over thetas
for theta = (min(thetaRange) + thetaStep / 2) : ...
            thetaStep : ...
            (max(thetaRange) - thetaStep / 2),

  % determine number of phis
  numberOfPhis = max(ceil(diff(phiRange) * sin(theta) / angularResolution), 1);

  % determine phi spacing
  phiStep = diff(phiRange) / numberOfPhis;

  % determine differential solid angle
  solidAngle = abs(sin(theta)) * thetaStep * phiStep;

  % append new thetas
  thetas = [thetas theta * ones(1, numberOfPhis)];

  % append new phis
  phis = [phis (min(phiRange) + phiStep / 2) : ...
          phiStep : ...
          (max(phiRange) - phiStep / 2)];

  % append new differential solid angles
  solidAngles = [solidAngles solidAngle * ones(1, numberOfPhis)];

% end loop over thetas
end

% force row vectors
thetas = thetas(:);
phis = phis(:);
solidAngles = solidAngles(:);

% vector of a priori uniform probability densities
probabilities = solidAngles / (diff(thetaRange) * diff(phiRange));

% construct coordinate matrix
coordinates = [thetas phis];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          return to calling function                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% return to calling function
return;