function [xs,ys,zs] = EarthCoord( q1,q2,Ielem,Welem,Oelem )
% [xs,ys,zs] = EarthCoord(q1,q2,Ielem,Welem,Oelem)
%
% q1     in-plane position or velocity component to perigee
% q2     in-plane position or velocity component perpendicular
% Ielem  inclination of the orbital plane
% Welem  argument of perigee
% Oelem  Right ascension of the ascending node
% xs     inertial x position or x velocity component
% ys     inertial y position or y velocity component
% zs     inertial z position or z velocity component 
c1 = -cos(Ielem)*sin(Welem)*sin(Oelem) + cos(Welem)*cos(Oelem); 
c2 = -cos(Ielem)*cos(Welem)*sin(Oelem) - sin(Welem)*cos(Oelem);
c3 =  cos(Ielem)*sin(Welem)*cos(Oelem) + cos(Welem)*sin(Oelem);
c4 =  cos(Ielem)*cos(Welem)*cos(Oelem) - sin(Welem)*sin(Oelem);
c5 =  sin(Ielem)*sin(Welem);
c6 =  sin(Ielem)*cos(Welem);
xs = q1*c1 + q2*c2;
ys = q1*c3 + q2*c4;
zs = q1*c5 + q2*c6;
%
% Be careful, Seeber has some typos in the rotated equations