function [satp] = satposin(t,eph)
%  SATPOSIN   Calculation of X,Y,Z coordinates in an INERTIAL
%	      reference frame at time t for given ephemeris eph

% Written by Kai Borre
% November 15, 1996

   GM = 3.986008e14;		% earth's universal gravitational parameter
				% m^3/s^2

   %  Units are either seconds, meters, or radians
   %  Assigning the local variables to eph
   svprn   =   eph(1);
   af2	   =   eph(2);
   M0	   =   eph(3);
   roota   =   eph(4);
   deltan  =   eph(5);
   ecc	   =   eph(6);
   omega   =   eph(7);
   cuc	   =   eph(8);
   cus	   =   eph(9);
   crc	   =  eph(10);
   crs	   =  eph(11);
   i0	   =  eph(12);
   idot    =  eph(13);
   cic	   =  eph(14);
   cis	   =  eph(15);
   Omega0  =  eph(16);
   Omegadot=  eph(17);
   toe	   =  eph(18);
   af0	   =  eph(19);
   af1	   =  eph(20);
   t0c	   =  eph(21);

   % Procedure for coordinate calculation
   A = roota*roota;
   t = check_t(t-t0c);
   n0 = sqrt(GM/A^3);
   n = n0+deltan;
   M = M0+n*t;
   M = rem(M,2*pi);
   E = M;
   E_old = E+.001;
   while abs(E-E_old) >= 1.e-12
      E_old = E;
      E = M+ecc*sin(E);
   end
   v = atan2(sqrt(1-ecc^2)*sin(E), cos(E)-ecc);
   phi = v+omega;
   phi = rem(phi,2*pi);
   u = phi		      + cuc*cos(2*phi)+cus*sin(2*phi);
   r = A*(1-ecc*cos(E))       + crc*cos(2*phi)+crs*sin(2*phi);
   i = i0+idot*t	      + cic*cos(2*phi)+cis*sin(2*phi);
   Omega = Omega0+Omegadot*t;
   x1 = cos(u)*r;
   y1 = sin(u)*r;
   satp(1,1) = x1*cos(Omega)-y1*cos(i)*sin(Omega);
   satp(2,1) = x1*sin(Omega)+y1*cos(i)*cos(Omega);
   satp(3,1) = y1*sin(i);

%%%%%%%%% end satposin.m %%%%%%%%%
