function handles = qskymap(coordinates, data, projection, mapNames, ...
                           unitNames, channelNames, colormapScale)
% QSKYMAP Perform spherical or other projection of sky map data
%
% QSKYMAP displays the specified sky map data for viewing and printing.  Two
% classes of projections are provided: a three dimensional onto a sphere, which
% may be interactively manipulated in the Matlab figure window, and the large
% variety of two dimensional projections provided in the Matlab mapping toolbox.
%
% usage:
%
%   handle = qskymap(coordinates, data, projection, mapNames, ...
%                    unitNames, channelNames, colormapScale);
%
%   coordinates       matrix of sky positions [radians]
%   data              matrix or cell array of sky map data
%   projection        projection type ('sphere', 'mollweid', etc.);
%   mapNames          cell array of map names for figure titles
%   unitNames         cell array of unit names for colorbar labels
%   channelNames      cell array of channel name strings for labels
%   colormapScale     range of data values corresponding to colormap
%
%   handles           figure handles for projected images
%
% The sky positions should be provided as a two column matrix of spherical
% coordinates in the form [theta phi] as returned by QTILESKY.  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).  Additional columns in the coordinate matrix
% are permitted, but are simply ignored.
%
% If data is a matrix or a cell array, a separate figure is generated for each
% column or cell of data.  Each row specifies the data at the sky position of
% the corresponding row in the coordinate matrix.
%
% A projection type of 'sphere' produces a three dimensional plot of the data on
% a sphere, which can be interactively manipulated in Matlab.  Other projection
% types are handled by the Matlab mapping toolbox.  In particular, the
% projection type 'mollweid' produces the standard equal area Mollweide
% projection that is commonly used in astrophysical publications.  See MAPPROJ
% for a list of other available projections.  If no projection type is
% specified, QSKYMAP assumes 'sphere' by default.
%
% The optional cell array of map names strings is used to label the figure
% tiles.  If no map names are specified, the figures are not titled.  Similarly,
% the optional cell array of unit names are used to label the colorbars.  If no
% unit names are specified, colorbars are not displayed.
%
% The optional cell array of channel name strings is used to label the detector
% zenith directions and pairwise detector baseline directions on spherical sky
% map plots.  If no channel names are specified, only the cartesian coordinate
% axes are labelled.  No labels are produced for two dimensional projections.
%
% The specified colormap scale can be used to explicitly specify the range of
% data values covered by the colormap.  Data values outside of the specified
% range are mapped to the corresponding extrema of the colormap.  By default,
% the colormap is autoscaled to the range of values in the data.
%
% See also QTILESKY, QRESPONSE, and MAPPROJ.

%   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: qskymap.m,v 1.3 2007/07/28 20:48:43 shourov Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               image parameters                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% map boundary
mapLeft = 0.14;
mapWidth = 0.80;
mapBottom = 0.14;
mapHeight = 0.72;
mapPosition = [mapLeft mapBottom ...
               mapWidth mapHeight];

% colorbar position
colorbarLeft = mapLeft;
colorbarWidth = mapWidth;
colorbarBottom = 0.12;
colorbarHeight = 0.02;
colorbarPosition = [colorbarLeft colorbarBottom ...
                    colorbarWidth colorbarHeight];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               data for labels                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% labelling parameters
vectorRadius = 1.15;
labelRadius = 1.20;

% cartesian coordinate axis
X.zenith = [1 0 0];
Y.zenith = [0 1 0];
Z.zenith = [0 0 1];

% cartesian direction of detector zeniths
G.zenith = [+0.603505859; +0.104172766; +0.790524328;];
H.zenith = [-0.338402190; -0.599658144; +0.725185541;];
L.zenith = [-0.011751435; -0.861335199; +0.507901150;];
T.zenith = [-0.618003751; +0.527187743; +0.583219039;];
V.zenith = [+0.711661388; +0.131946680; +0.690020392;];

% cartesian location of detector corner stations
G.location = [+3.85631120e+06; +6.66597800e+05; +5.01964060e+06;];
H.location = [-2.16141493e+06; -3.83469518e+06; +4.60035022e+06;];
L.location = [-7.42760419e+04; -5.49628372e+06; +3.22425702e+06;];
T.location = [-3.94640900e+06; +3.36625900e+06; +3.69915100e+06;];
V.location = [+4.54637400e+06; +8.42990000e+05; +4.37857700e+06;];

% detector baselines
GH.baseline = G.location - H.location;
GL.baseline = G.location - L.location;
GT.baseline = G.location - T.location;
GV.baseline = G.location - V.location;
HL.baseline = H.location - L.location;
HT.baseline = H.location - T.location;
HV.baseline = H.location - V.location;
LT.baseline = L.location - T.location;
LV.baseline = L.location - V.location;
TV.baseline = T.location - V.location;
GH.baseline = GH.baseline / sqrt(sum(GH.baseline.^2));
GL.baseline = GL.baseline / sqrt(sum(GL.baseline.^2));
GT.baseline = GT.baseline / sqrt(sum(GT.baseline.^2));
GV.baseline = GV.baseline / sqrt(sum(GV.baseline.^2));
HL.baseline = HL.baseline / sqrt(sum(HL.baseline.^2));
HT.baseline = HT.baseline / sqrt(sum(HT.baseline.^2));
HV.baseline = HV.baseline / sqrt(sum(HV.baseline.^2));
LT.baseline = LT.baseline / sqrt(sum(LT.baseline.^2));
LV.baseline = LV.baseline / sqrt(sum(LV.baseline.^2));
TV.baseline = TV.baseline / sqrt(sum(TV.baseline.^2));

% plot arguments for cartesian coordinates
X.x = [0 1];
X.y = [0 0];
X.z = [0 0];
Y.x = [0 0];
Y.y = [0 1];
Y.z = [0 0];
Z.x = [0 0];
Z.y = [0 0];
Z.z = [0 1];

% plot arguments for detector zenith vectors
G.x = G.zenith(1) * [0 1];
G.y = G.zenith(2) * [0 1];
G.z = G.zenith(3) * [0 1];
H.x = H.zenith(1) * [0 1];
H.y = H.zenith(2) * [0 1];
H.z = H.zenith(3) * [0 1];
L.x = L.zenith(1) * [0 1];
L.y = L.zenith(2) * [0 1];
L.z = L.zenith(3) * [0 1];
T.x = T.zenith(1) * [0 1];
T.y = T.zenith(2) * [0 1];
T.z = T.zenith(3) * [0 1];
V.x = V.zenith(1) * [0 1];
V.y = V.zenith(2) * [0 1];
V.z = V.zenith(3) * [0 1];

% plot arguments for detector baseline vectors
GH.x = GH.baseline(1) * [-1 +1];
GH.y = GH.baseline(2) * [-1 +1];
GH.z = GH.baseline(3) * [-1 +1];
GL.x = GL.baseline(1) * [-1 +1];
GL.y = GL.baseline(2) * [-1 +1];
GL.z = GL.baseline(3) * [-1 +1];
GT.x = GT.baseline(1) * [-1 +1];
GT.y = GT.baseline(2) * [-1 +1];
GT.z = GT.baseline(3) * [-1 +1];
GV.x = GV.baseline(1) * [-1 +1];
GV.y = GV.baseline(2) * [-1 +1];
GV.z = GV.baseline(3) * [-1 +1];
HL.x = HL.baseline(1) * [-1 +1];
HL.y = HL.baseline(2) * [-1 +1];
HL.z = HL.baseline(3) * [-1 +1];
HT.x = HT.baseline(1) * [-1 +1];
HT.y = HT.baseline(2) * [-1 +1];
HT.z = HT.baseline(3) * [-1 +1];
HV.x = HV.baseline(1) * [-1 +1];
HV.y = HV.baseline(2) * [-1 +1];
HV.z = HV.baseline(3) * [-1 +1];
LT.x = LT.baseline(1) * [-1 +1];
LT.y = LT.baseline(2) * [-1 +1];
LT.z = LT.baseline(3) * [-1 +1];
LV.x = LV.baseline(1) * [-1 +1];
LV.y = LV.baseline(2) * [-1 +1];
LV.z = LV.baseline(3) * [-1 +1];
TV.x = TV.baseline(1) * [-1 +1];
TV.y = TV.baseline(2) * [-1 +1];
TV.z = TV.baseline(3) * [-1 +1];

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

% verify correct number of input arguments
error(nargchk(2, 7, nargin));

% apply default arguments
if (nargin < 3) || isempty(projection),
  projection = 'sphere';
end
if (nargin < 4) || isempty(mapNames),
  mapNames = [];
end
if (nargin < 5) || isempty(unitNames),
  unitNames = [];
end
if (nargin < 6) || isempty(channelNames),
  channelNames = {};
end
if (nargin < 7) || isempty(colormapScale),
  colormapScale = [];
end

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

% validate coordinate matrix
if size(coordinates, 2) < 2,
  error('coordinates must be at least a two column matrix [theta phi]');
end

% extract spherical coordinates
theta = coordinates(:, 1);
phi = coordinates(:, 2);

% number of sky positions
numberOfSkyPositions = length(theta);

% force cell array of data;
if ~iscell(data),
  data = mat2cell(data, size(data, 1), ones(1, size(data, 2)));
end

% number of maps
numberOfMaps = length(data);

% validate data lengths
for mapNumber = 1 : numberOfMaps,
  if length(data{mapNumber}) ~= numberOfSkyPositions,
    error('data length inconsistent with number of sky positions');
  end
end

% validate spherical coordinates
if any((theta < 0) | (theta > pi)),
  error('theta outside of range [0, pi]');
end
if any((phi < -pi) | (phi >= pi)),
  error('phi outside of range [-pi, pi)');
end

% default map names
if isempty(mapNames),
  mapNames = cell(1, numberOfMaps);
  % for mapNumber = 1 : numberOfMaps,
  %   mapNames{mapNumber} = ['Map ' num2str(mapNumber)];
  % end
end

% force cell array of map names
if ~iscell(mapNames),
  mapNames = mat2cell(mapNames, size(mapNames, 1), size(mapNames, 2));
end

% handle case of single map name
if length(mapNames) == 1,
  for mapNumber = 2 : numberOfMaps,
    mapNames(mapNumber) = mapNames(1);
  end
end

% force one dimensional cell array of map names
mapNames = mapNames(:);

% validate number of map names
if length(mapNames) ~= numberOfMaps,
  error('number of map names inconsistent with number of maps');
end

% default unit names
if isempty(unitNames),
  unitNames = cell(1, numberOfMaps);
end

% force cell array of unit names
if ~iscell(unitNames),
  unitNames = mat2cell(unitNames, size(unitNames, 1), size(unitNames, 2));
end

% handle case of single unit name
if length(unitNames) == 1,
  for mapNumber = 2 : numberOfMaps,
    unitNames(mapNumber) = unitNames(1);
  end
end

% force one dimensional cell array of unit names
unitNames = unitNames(:);

% validate number of unit names
if length(unitNames) ~= numberOfMaps,
  error('number of unit names inconsistent with number of maps');
end

% force cell array of channel names
if ~iscell(channelNames),
  channelNames = mat2cell(channelNames, size(channelNames, 1), ...
                          size(channelNames, 2));
end

% force one dimensional cell array of channel names
channelNames = channelNames(:);

% number of channel names
numberOfChannels = length(channelNames);

% force lowercase projection type name
projection = lower(projection);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               determine labels                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% number of labels
numberOfLabels = 3 + numberOfChannels + ...
                 numberOfChannels * (numberOfChannels - 1) / 2;

% initialize labels
labels = cell(numberOfLabels, 1);

% add cartesian coordinate axes to labels
labels{1} = 'X';
labels{2} = 'Y';
labels{3} = 'Z';

% add detector zeniths to labels
labelNumber = 4;
for channelNumber = 1 : numberOfChannels,
  labels{labelNumber} = channelNames{channelNumber}(1);
  labelNumber = labelNumber + 1;
end

% add detector baselines to labels
for primaryChannelNumber = 2 : numberOfChannels,
  for secondaryChannelNumber = 1 : primaryChannelNumber - 1,
    labels{labelNumber} = [channelNames{primaryChannelNumber}(1) ...
                           channelNames{secondaryChannelNumber}(1)];
    labelNumber = labelNumber + 1;
  end
end

% force uppercase labels
labels = upper(labels);

% alphabetically order baseline labels
for labelNumber = 1 : numberOfLabels,
  labels{labelNumber} = sort(labels{labelNumber});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             begin loop over maps                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize handles
handles = [];

% begin loop over maps
for mapNumber = 1 : numberOfMaps,

  % initialize figure window
  if numberOfMaps > 1,
    figure(mapNumber);
  end
  clf;
  set(gca, 'FontSize', 18);

  % add figure handle to vector of figure handles
  handles = [handles gca];

  % determine colormap scale
  if isempty(colormapScale),
    colormapScale = [min(data{mapNumber}) max(data{mapNumber})];
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                     begin switch on projection type                        %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % begin switch on projection type
  switch projection,

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                         spherical projection                             %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % handle case of spherical projection
    case 'sphere',

      % cartesian coordinates of sphere
      x = sin(theta) .* cos(phi);
      y = sin(theta) .* sin(phi);
      z = cos(theta);

      % create triangular mesh
      triangulation = convhulln([x, y, z]);

      % display map
      trisurf(triangulation, x, y, z, data{mapNumber});

      % add labels
      hold on;
      for labelNumber = 1 : numberOfLabels,
        eval(['plot3(vectorRadius * ' labels{labelNumber} '.x, ' ...
              '      vectorRadius * ' labels{labelNumber} '.y, ' ...
              '      vectorRadius * ' labels{labelNumber} '.z, ''k-'');']);;
        eval(['textHandle = text(labelRadius * ' labels{labelNumber} '.x, ' ...
              '                  labelRadius * ' labels{labelNumber} '.y, ' ...
              '                  labelRadius * ' labels{labelNumber} '.z, ' ...
              '                  labels{labelNumber});']);;
        set(textHandle, 'HorizontalAlignment', 'center');
      end
      hold off;

      % enable three dimensional rotation
      rotate3d on;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                      two dimensional projection                          %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % otherwise, perform two dimensional projection
    otherwise,

      % perform projection
      mapstruct = defaultm(projection);
      mapstruct.angleunits = 'radians';
      mapstruct.origin = [0 0 0];
      mapstruct.falseeasting = 0;
      mapstruct.falsenorthing = 0;
      mapstruct.scalefactor = 1;
      [x, y] = mfwdtran(mapstruct, pi / 2 - theta, phi);
      z = zeros(size(x));

      % create triangular mesh
      triangulation = delaunay(x, y);

      % display map
      trisurf(triangulation, x, y, z, data{mapNumber});

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                      end switch on projection type                         %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % end switch on projection type
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                           set image properties                             %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % set colormap scale and enable interpolation
  colormap('default');
  caxis(colormapScale);
  shading interp;

  % set axis position
  set(gca, 'Position', mapPosition);

  % equalize axis scales and turn off axis labelling
  axis image;
  axis vis3d;
  axis off;

  % set viewpoint
  view([0 0 1]);

  % set title
  % title(strrep(mapNames{mapNumber}, '_', '\_'));
  title(mapNames{mapNumber});

  % display colorbar
  if ~isempty(unitNames{mapNumber}),
    subplot('position', colorbarPosition);
    set(gca, 'FontSize', 16);
    colorbarmap = linspace(min(colormapScale), max(colormapScale), 100);
    imagesc(colorbarmap, 1, colorbarmap, colormapScale);
    set(gca, 'YTick',[])
    set(gca, 'TickDir', 'out')
    % xlabel(strrep(unitNames{mapNumber}, '_', '\_'));
    xlabel(unitNames{mapNumber});
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              end loop over maps                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% end loop over maps
end

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

% return to calling function
return;