function htt = makett(hij, Internal)

% MAKETT - projects metric perturbation into TT gauge
%
% htt = makett(hij, nDet);
%
% hij       metric perturbation for gravitational wave source in 3x3xtime
%           format
% Internal  struct describing unit vector in direction of wave propagation
%           (i.e., from source to detector)
%   .x      cos angle with respect to source Z-axis
%   .phi    azimuth angle from source x axis
%
% htt       metric perturbation projected into TT gauge in 2 x time format
%       htt(1,:) - plus polarization
%       htt(2,:) - cross polarization
%   q.v. LIGO-T040020-Z for + and x conventions
%
% $Id: makett.m,v 1.17 2005/04/12 21:23:48 stuver Exp $
%
% SEE ALSO: GRAVEN

% INITIALIZE output variable
% CHECK inputs for errors
% DEFINE angles
% DEFINE rotation matrix
% DEFINE plus and cross bases
% PROJECT metric pertubation onto bases
% RETURN pluss and cross polarization timeseries

% Check for errors
errorcheck(hij, Internal);

% Direction *from* source: i.e., -Internal
cosTheta = -Internal.x;
sinTheta = sqrt(1-cosTheta^2);
phi = Internal.phi + pi;

% Form X, Y coordinate axes so that [X, Y, -Internal] is right-handed
% coordinate system. 
%
% GOTCHA: Note the -Internal above! The orientation of the X, Y that span the
% polarization plane must be the same here as in detproj. 

X = [cos(phi)*cosTheta, sin(phi)*cosTheta, -sinTheta];
Y = [-sin(phi), cos(phi), 0];
  
% Form the plus and cross bases
% eplus =  X'*X - Y'*Y;
% ecross = X'*Y + Y'*X;

% For each timestep, project the metric perturbation onto the plus and
% cross basis
% h+ = 0.5*(h^{ij}eplus_{ij})
% hx = 0.5*(h^{ij}ecross_{ij})
% * N.B. Multiply h^{ij}e_{ij} by 0.5 since dot(e_{ij},e_{ij}) = 2

% htt = zeros([2,size(hij,3)]);
% for k = 1:size(hij,3)
%     htt(1,k) = 0.5*sum(sum(hij(:,:,k).*eplus));
%     htt(2,k) = 0.5*sum(sum(hij(:,:,k).*ecross));
% end

% Let's be clever and use reshape and permute

nSamp = size(hij,3);
htt = zeros([2,nSamp]);
hij = reshape(hij,[3,3*nSamp]);
xh = permute(reshape(X*hij,[3,nSamp]),[2,1]);
yh = permute(reshape(Y*hij,[3,nSamp]),[2,1]);
htt(1,:) = (xh*X'-yh*Y')'/2;
htt(2,:) = (xh*Y'+yh*X')'/2;


return

% ~~~ END MAIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function errorcheck(hij, Internal)

% ERRORCHECK - checks makeh's input data for errors and stops the code if
%              an error is found
%
% errorcheck(hij, Internal)
%
% hij       metric perturbation for gravitational wave source
% Interanl       Direction of wave propagation
%   .x       cos co-declination
%   .phi     source longitude


% DETERMINE size of hij
% CHECK all inputs for errors

dims = size(hij);
if (dims(1) ~= 3 || dims(2) ~= 3)
    % size(hij) must be in [3, 3, nSamp]
    if (isempty(hij))
      error('makett: hij empty');
    else
      errstr = sprintf('%d',dims(1));
      if (length(dims) > 1)
        errstr = [errstr, sprintf(', %d',dims(2:end))];
      end
      error('makett: size(hij) = [%s] should be [3, 3, nSamp]',errstr);
    end
end

% nDet must contain fields x, phi and psi; these must be in limits
if ~(isfield(Internal, 'x') && -1 <= Internal.x && Internal.x <= 1)
    error('detproj: scalar nSrc.x must satisfy -1 <= nSrc.x <= 1')
end
if ~(isfield(Internal, 'phi') && -pi <= Internal.phi && Internal.phi <= pi)
    error('detproj: scalar nSrc.phi must satisfy -pi <= nSrc.phi <= pi');
end
return