function [fplus, fcross, deltat] = qresponse(coordinates, detectors)
% QRESPONSE Antenna response for a given detector and source coordinates
%
% QRESPONSE computes the response of the specified gravitational wave
% detector to both polarizations of gravitational wave signal coming
% from the specified directions on the sky with 
%
%  usage:
%
%    [fplus, fcross, deltat] = qresponse(coordinates, detectors);
%
%    coordinates    matrix of sky positions and polarization angles [radians]
%    detector       cell array of site, detector, or channel names
%
%    fplus          matrix of plus polarization antenna response factors
%    fcross         matrix of cross polarization antenna response factors
%    deltat         matrix of time delays relative to earth center
%
% The sky positions and polarizations should be provided as a three column
% matrix of spherical coordinates in the form [theta phi psi].  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 polarization psi is defined as
% the clockwise angle about the vector pointing towards the source from the the
% negative theta direction to the plus polarization direction of the source.
% The range of theta is [0, pi] and the range of phi is [-pi, pi).  The range of
% psi is [0, pi).  If the coordinate matrix has only two columns, they are
% assumed to be theta and phi values and psi is assumed to be zero.  If the
% coordinate matrix has more than three columns, the additional columns are
% simply ignored.
%
% The list of detectors should be specified as a cell array of site identifiers
% or channel names.  Site identifiers are determined by taking the first
% character of each channel name.  The following sites are recognized.
%
%   G  GEO
%   H  LIGO Hanford
%   L  LIGO Livingston
%   V  Virgo
%
% The resulting fplus, fcross, deltat matrices contain a column for each detector
% and a row for each coordinate.
%
% See also QTILESKY 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: qresponse.m,v 1.3 2007/07/26 20:53:04 shourov Exp $

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

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

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

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

% extract spherical coordinates and polarizations
theta = coordinates(:, 1);
phi = coordinates(:, 2);
psi = coordinates(:, 3);

% validate 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
if any((psi < 0) | (psi >= pi)),
  error('phi outside of range [0, pi)');
end

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

% force one dimensional cell array
detectors = detectors(:);

% number of detectors
numberOfDetectors = length(detectors);

% force uppercase detector names
detector = upper(detectors);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             direction to source                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cartesian source direction
sourceDirection = [sin(theta) .* cos(phi) ...
                   sin(theta) .* sin(phi) ...
                   cos(theta)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             polarization tensors                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute polarization tensors
m1 = sin(phi) .* sin(psi) + cos(phi) .* cos(theta) .* cos(psi);
m2 = -cos(phi) .* sin(psi) + sin(phi) .* cos(theta) .* cos(psi);
m3 = -sin(theta) .* cos(psi);
n1 = sin(phi) .* cos(psi) - cos(phi) .* cos(theta) .* sin(psi);
n2 = - cos(phi) .* cos(psi) - sin(phi) .* cos(theta) .* sin(psi);
n3 = sin(theta) .* sin(psi);
mm = [m1.*m1, m1.*m2, m1.*m3, m2.*m1, m2.*m2, m2.*m3, m3.*m1, m3.*m2, m3.*m3];
mn = [m1.*n1, m1.*n2, m1.*n3, m2.*n1, m2.*n2, m2.*n3, m3.*n1, m3.*n2, m3.*n3];
nm = [n1.*m1, n1.*m2, n1.*m3, n2.*m1, n2.*m2, n2.*m3, n3.*m1, n3.*m2, n3.*m3];
nn = [n1.*n1, n1.*n2, n1.*n3, n2.*n1, n2.*n2, n2.*n3, n3.*n1, n3.*n2, n3.*n3];
eplus = mm - nn;
ecross = mn + nm;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          begin loop over detectors                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialize results
fplus = [];
fcross = [];
deltat = [];

% begin loop over detectors
for detectorNumber = 1 : numberOfDetectors,

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                          detector response data                            %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  switch detector{detectorNumber}(1),

    % handle case of GEO600 detector
    case 'G',

      % cartesian coordinates of detector [light seconds]
      detectorLocation = ...
          [+0.012863270; +0.002223531; +0.016743719];

      % detector response matrix
      detectorResponse = ...
          [-0.096843926; -0.365712231; +0.122125440; ...
           -0.365712231; +0.222996416; +0.249808025; ...
           +0.122125440; +0.249808025; -0.126152489;];

    % handle cases of LIGO Hanford 2k and 4km detectors
    case 'H',

      % cartesian coordinates of detector [light seconds]
      detectorLocation = ...
          [-0.007209704; -0.012791166; +0.015345117;];

      % detector response matrix
      detectorResponse = ...
          [-0.392614702; -0.077612253; -0.247388405; ...
           -0.077612253; +0.319524089; +0.227998294; ...
           -0.247388405; +0.227998294; +0.073090613;];

    % handle case of LIGO Livingston 2k detector
    case 'L',

      % cartesian coordinates of detector [light seconds]
      detectorLocation = ...
          [-0.000247758; -0.018333629; +0.010754964;];

      % detector response matrix
      detectorResponse = ...
          [+0.411281744; +0.140209630; +0.247293475; ...
           +0.140209630; -0.109005943; -0.181616031; ...
           +0.247293475; -0.181616031; -0.302275801;];

    % handle case of Virgo detector
    case 'V',

      % cartesian coordinates of detector [light seconds]
      detectorLocation = ...
          [+0.015165071; +0.002811912; +0.014605361;];

      % detector response matrix
      detectorResponse = ...
          [+0.243874678; -0.099086615; -0.232575796; ...
           -0.099086615; -0.447827872; +0.187828535; ...
           -0.232575796; +0.187828535; +0.203953193;];

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                            detector response                               %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine detector response
  fplus = [fplus eplus * detectorResponse];
  fcross = [fcross ecross * detectorResponse];
  deltat = [deltat -sourceDirection * detectorLocation];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                         end loop over detectors                            %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% end loop over detectors
end

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

% return to calling function
return;