%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Flemming Hansen, Danish Space Research Institute, Copenhagen, Denmark
%
%  Version 1.0                                                FH 2003-12-16
%  Version 1.1 Small modification of the eci2ecef.m function HPF 2006-09-20
%  Version 1.2 Enroll in Revision Control System (rcs)       LTC 2006-10-08
%
%  Point of contact for all matters regarding this software:
%  Flemming Hansen, Danish Space Research Institute
%  Phone: +45 3532 5712, e-mail: fh@dsri.dk
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  File:  icrf2itrf.m 
%
%  Description: Converts vector or array of vectors in ECI Earth Centered
%  Inertial (ECI = ICRS) coordinates to cartesian Earth Centered Earth
%  Fixed (ECEF = ITRS) coordinates by rotation around the z-axis
% 
%  Initialization: None required
% 
%  Input:
%    JD       time expressed as JD (Julian day): the numbers of days that
%             has elapsed since Greenwich mean noon(12) on
%             January 1, 4713 B.C (Gregorian calendar) plus a fractional
%             part of the day elapsed since the preceding noon. (nx1 - size)
%    r_ECI    Vectors or positions in ECI/ICRF frame (nx3 - size)
%             If empty, quaternions are returned
%    Aux      Auxiliary information (Earth Rotation Parameters) (struct)
% 
%  Output:
%    r_ECEF   Vectors or positions in ECEF/ITRF frame (nx3 - size)
%    q        Quaternion: ITRF <- ICRF (nx4 - size)
%
%  Call:  
%    r_ECEF = icrf2itrf(JD, r_ECI, Aux);
%    q      = icrf2itrf(JD, [], Aux);
%
%  $Revision: 1.2 $
%

% $Id: icrf2itrf.m,v 1.2 2006/10/09 08:14:49 lastec Exp $

function res = icrf2itrf(JD, r_ECI, Aux);

% Calculate the Modified Julian Date base 2000
tmjd2k = JD - 2400000.5 - 51544;

% Calculate Greenwich Aparent Siderial Time
gast = gast_mjd2000(tmjd2k);

% Compute or do the rotation
if((exist('r_ECI', 'var') == 0) || (numel(r_ECI) == 0))
  % Return quaternions - rotation around z axis is easy!
  res = zeros(numel(gast), 4);
  res(:,3) = sin(0.5 * gast);
  res(:,4) = cos(0.5 * gast);
else
  cos_g = cos(gast);
  sin_g = sin(gast);
  res = [cos_g.*r_ECI(:,1)-sin_g.*r_ECI(:,2) ...
         sin_g.*r_ECI(:,1)+cos_g.*r_ECI(:,2) ...
         r_ECI(:,3)];
end%if

return;

