Index: /trunk/BNS/glonass.cpp
===================================================================
--- /trunk/BNS/glonass.cpp	(revision 883)
+++ /trunk/BNS/glonass.cpp	(revision 883)
@@ -0,0 +1,53 @@
+/* -------------------------------------------------------------------------
+ * BKG NTRIP Server
+ * -------------------------------------------------------------------------
+ *
+ * Function:   glo_deriv
+ *
+ * Purpose:    Derivative of the state vector of a Galileo satellite
+ *             using its position, velocity and a simplified force model
+ *
+ * Author:     L. Mervart
+ *
+ * Created:    07-Mai-2008
+ *
+ * Changes:    
+ *
+/******************************************************************************/
+
+#include "glonass.h" 
+
+// Derivative of the state vector
+////////////////////////////////////////////////////////////////////////////
+void glo_deriv(double /* tt */, const ColumnVector& yy, 
+               ColumnVector& yp, void* /* pVoid */) {
+
+  // State vector components
+  // -----------------------
+  ColumnVector rr = yy.rows(1,3);
+  ColumnVector vv = yy.rows(4,6);
+
+  // Acceleration 
+  // ------------
+  const static double GM    = 398.60044e12;
+  const static double AE    = 6378136.0;
+  const static double OMEGA = 7292115.e-11;
+  const static double C20   = -1082.63e-6;
+
+  double rho = rr.norm_Frobenius();
+  double t1  = -GM/(rho*rho*rho);
+  double t2  = 3.0/2.0 * C20 * (GM*AE*AE) / (rho*rho*rho*rho*rho);
+  double t3  = OMEGA * OMEGA;
+  double t4  = 2.0 * OMEGA;
+  double z2  = rr(3) * rr(3);
+
+  ColumnVector aa(3);
+  aa(1) = (t1 + t2*(1.0-5.0*z2/(rho*rho)) + t3) * rr(1) + t4*vv(2); 
+  aa(2) = (t1 + t2*(1.0-5.0*z2/(rho*rho)) + t3) * rr(2) - t4*vv(1); 
+  aa(3) = (t1 + t2*(3.0-5.0*z2/(rho*rho))     ) * rr(3);
+
+  // State vector derivative
+  // -----------------------  
+  yp = vv &
+       aa ;
+}
Index: /trunk/BNS/glonass.h
===================================================================
--- /trunk/BNS/glonass.h	(revision 883)
+++ /trunk/BNS/glonass.h	(revision 883)
@@ -0,0 +1,9 @@
+#ifndef GLONASS_H
+#define GLONASS_H
+
+#include <newmat.h>
+
+void glo_deriv(double tt, const ColumnVector& yy, 
+               ColumnVector& yp, void* pVoid = 0);
+
+#endif
Index: unk/BNS/rungekutta4.cpp
===================================================================
--- /trunk/BNS/rungekutta4.cpp	(revision 882)
+++ 	(revision )
@@ -1,94 +1,0 @@
-/* -------------------------------------------------------------------------
- * BKG NTRIP Server
- * -------------------------------------------------------------------------
- *
- * Function:   rungeKutta4
- *
- * Purpose:    Fourth order Runge-Kutta numerical integrator for ODEs
- *
- * Example usage - numerical integration dy/dt = cos(x) :
- *
- *       std::vector<double> dydtcos(double x, vector<double> y) {
- *           std::vector<double> dydt(y)
- *           dydt.at(i)=cos(x);
- *           return dydt;
- *       }
- *       
- *       vector<double> yi(1,0.0);
- *       double x=0.0, dx=0.10;
- *       
- *       while (x<10.0) {
- *           yi=rungeKutta4(yi,x,dx,(*dydtcos));
- *           x+=dx;
- *       }
- *
- * Author:     L. Mervart
- *
- * Created:    07-Mai-2008
- *
- * Changes:    
- *
-/******************************************************************************/
-
-#include "rungekutta4.h" 
-
-using namespace std;
-
-vector<double> rungeKutta4(
-  double xi,         // the initial x-value
-  vector<double> yi, // vector of the initial y-values
-  double dx,         // the step size for the integration
-  vector<double> (*derivatives)(double x, vector<double> y) ) {
-                     // a pointer to a function that returns the derivative 
-                     // of a function at a point (x,y), where y is an STL 
-                     // vector of elements.  Returns a vector of the same 
-                     // size as y.
-
-  //  total number of elements in the vector
-  int n=yi.size();
-  
-  //  first step
-  vector<double> k1;
-  k1=derivatives(xi, yi);
-  for (int i=0; i<n; ++i) {
-    k1.at(i)*=dx;
-  }
-  
-  //  second step
-  vector<double> k2(yi);
-  for (int i=0; i<n; ++i) {
-    k2.at(i)+=k1.at(i)/2.0;
-  }
-  k2=derivatives(xi+dx/2.0,k2);
-  for (int i=0; i<n; ++i) {
-    k2.at(i)*=dx;
-  }
-  
-  //  third step
-  vector<double> k3(yi);
-  for (int i=0; i<n; ++i) {
-    k3.at(i)+=k2.at(i)/2.0;
-  }
-  k3=derivatives(xi+dx/2.0,k3);
-  for (int i=0; i<n; ++i) {
-    k3.at(i)*=dx;
-  }
-  
-  //  fourth step
-  vector<double> k4(yi);
-  for (int i=0; i<n; ++i) {
-    k4.at(i)+=k3.at(i);
-  }
-  k4=derivatives(xi+dx,k4);
-  for (int i=0; i<n; ++i) {
-    k4.at(i)*=dx;
-  }
-  
-  //  sum the weighted steps into yf and return the final y values
-  vector<double> yf(yi);
-  for (int i=0; i<n; ++i) {
-    yf.at(i)+=(k1.at(i)/6.0)+(k2.at(i)/3.0)+(k3.at(i)/3.0)+(k4.at(i)/6.0);
-  }
-  
-  return yf;
-}
Index: unk/BNS/rungekutta4.h
===================================================================
--- /trunk/BNS/rungekutta4.h	(revision 882)
+++ 	(revision )
@@ -1,13 +1,0 @@
-#ifndef RUNGEKUTTA4_H
-#define RUNGEKUTTA4_H
-
-#include <vector>
-
-std::vector<double> rungeKutta4(
-    double xi,
-    std::vector<double> yi,
-    double dx,
-    std::vector<double> (*derivatives)(double x, std::vector<double> y) 
-);
-
-#endif