Index: /trunk/BNS/bns.cpp
===================================================================
--- /trunk/BNS/bns.cpp	(revision 801)
+++ /trunk/BNS/bns.cpp	(revision 802)
@@ -235,6 +235,8 @@
 
   ColumnVector xB(4);
-
-  satellitePosition(GPSweeks, ep, xB(1), xB(2), xB(3), xB(4));
+  ColumnVector vv(3);
+
+  satellitePosition(GPSweeks, ep, xB(1), xB(2), xB(3), xB(4), 
+                    vv(1), vv(2), vv(3));
 
   ColumnVector dx = xx - xB;
Index: /trunk/BNS/bnsutils.cpp
===================================================================
--- /trunk/BNS/bnsutils.cpp	(revision 801)
+++ /trunk/BNS/bnsutils.cpp	(revision 802)
@@ -97,5 +97,6 @@
 ////////////////////////////////////////////////////////////////////////////
 void satellitePosition(double GPSweeks, const gpsEph* ep, 
-                       double& X, double& Y, double& Z, double& dt) {
+                       double& X, double& Y, double& Z, double& dt,
+                       double& vX, double& vY, double& vZ) {
 
   const static double omegaEarth = 7292115.1467e-11;
@@ -142,3 +143,27 @@
   dt = ep->clock_bias + ep->clock_drift*tc + ep->clock_driftrate*tc*tc 
        - 4.442807633e-10 * ep->e * sqrt(a0) *sin(E);
+
+  // Velocity
+  // --------
+  double tanv2 = tan(v/2);
+  double dEdM  = 1 / (1 - ep->e*cos(E));
+  double dotv  = sqrt((1.0 + ep->e)/(1.0 - ep->e)) / cos(E/2)/cos(E/2) / (1 + tanv2*tanv2) 
+               * dEdM * n;
+  double dotu  = dotv + (-ep->Cuc*sin2u0 + ep->Cus*cos2u0)*2*dotv;
+  double dotom = ep->OMEGADOT - omegaEarth;
+  double doti  = ep->IDOT + (-ep->Cic*sin2u0 + ep->Cis*cos2u0)*2*dotv;
+  double dotr  = a0 * ep->e*sin(E) * dEdM * n 
+                + (-ep->Crc*sin2u0 + ep->Crs*cos2u0)*2*dotv;
+  double dotx  = dotr*cos(u) - r*sin(u)*dotu;
+  double doty  = dotr*sin(u) + r*cos(u)*dotu;
+
+  vX  = cosom   *dotx  - cosi*sinom   *doty      // dX / dr
+      - xp*sinom*dotom - yp*cosi*cosom*dotom     // dX / dOMEGA
+                       + yp*sini*sinom*doti;     // dX / di
+
+  vY  = sinom   *dotx  + cosi*cosom   *doty
+      + xp*cosom*dotom - yp*cosi*sinom*dotom
+                       - yp*sini*cosom*doti;
+
+  vZ  = sini    *doty  + yp*cosi      *doti;
 }
Index: /trunk/BNS/bnsutils.h
===================================================================
--- /trunk/BNS/bnsutils.h	(revision 801)
+++ /trunk/BNS/bnsutils.h	(revision 802)
@@ -18,4 +18,5 @@
 
 void satellitePosition(double GPSweeks, const gpsEph* ep, 
-                       double& X, double& Y, double& Z, double& dt);
+                       double& X, double& Y, double& Z, double&,
+                       double& vX, double& vY, double& vZ);
 #endif