Index: trunk/BNC/src/pppModel.cpp
===================================================================
--- trunk/BNC/src/pppModel.cpp	(revision 7245)
+++ trunk/BNC/src/pppModel.cpp	(revision 7246)
@@ -44,5 +44,4 @@
 
 #include "pppModel.h"
-#include "bncutils.h"
 
 using namespace BNC_PPP;
@@ -387,18 +386,100 @@
 ///////////////////////////////////////////////////////////////////////////
 t_iono::t_iono() {
-  _vTec = 0;
-}
-
-t_iono::~t_iono() {
-  delete _vTec;
-}
-
-void t_iono::setTecData(const t_vTec* vTec) {
-   delete _vTec;
-  _vTec = new t_vTec(*vTec);
-}
-
-double t_iono::vtec() {
-
-  return 0.0;
-}
+  _psiPP = _phiPP = _lambdaPP = _lonS = 0.0;
+}
+
+t_iono::~t_iono() {}
+
+double t_iono::stec(const t_vTec* vTec, double signalPropagationTime,
+      const ColumnVector& rSat, const bncTime& epochTime,
+      const ColumnVector& xyzSta) {
+
+  // Latitude, longitude, height are defined with respect to a spherical earth model
+  // -------------------------------------------------------------------------------
+  ColumnVector geocSta;
+  xyz2geoc(xyzSta.data(), geocSta.data());
+
+  // satellite position rotated to the epoch of signal reception
+  // -----------------------------------------------------------
+  ColumnVector xyzSat;
+  double omegaZ = t_CST::omega * signalPropagationTime;
+  xyzSat[0] = rSat[0] * cos(omegaZ) + rSat[1] * sin(omegaZ);
+  xyzSat[1] = rSat[1] * cos(omegaZ) - rSat[0] * sin(omegaZ);
+  xyzSat[2] = rSat[2];
+
+  // elevation and azimuth with respect to a spherical earth model
+  // -------------------------------------------------------------
+  ColumnVector rhoV = xyzSat - xyzSta;
+  double rho = rhoV.norm_Frobenius();
+  ColumnVector neu(3);
+  xyz2neu(geocSta.data(), rhoV.data(), neu.data());
+  double sphEle = acos( sqrt(neu[0]*neu[0] + neu[1]*neu[1]) / rho );
+  if (neu[2] < 0) {
+    sphEle *= -1.0;
+  }
+  double sphAzi = atan2(neu[1], neu[0]);
+
+  double epoch = fmod(epochTime.gpssec(), 86400.0);
+
+  double stec = 0.0;
+  for (unsigned ii = 0; ii < vTec->_layers.size(); ii++) {
+    double layerHeight = vTec->_layers[ii]._height * 1000.0; // m
+    piercePoint(layerHeight, epoch, geocSta.data(), sphEle, sphAzi);
+    double vtec = vtecSingleLayerContribution(vTec->_layers[ii]);
+    stec += vtec * sin(sphEle * _psiPP);
+  }
+  return stec;
+}
+
+double t_iono::vtecSingleLayerContribution(const t_vTecLayer& vTecLayer) {
+
+  double vtec = 0.0;
+  int N = vTecLayer._C.Nrows()-1;
+  int M = vTecLayer._C.Ncols()-1;
+  double fac;
+
+  for (int n = 0; n <= N; n++) {
+    for (int m = 0; m <= min(n, M); m++) {
+      double pnm = associatedLegendreFunction(n, m, sin(_phiPP));
+      double a = double(factorial(n - m));
+      double b = double(factorial(n + m));
+      if (m == 0) {
+        fac = sqrt(2.0 * n + 1);
+      }
+      else {
+        fac = sqrt(2.0 * (2.0 * n + 1) * a / b);
+      }
+      pnm *= fac;
+      double Cnm_mlambda = vTecLayer._C[n][m] * cos(m * _lonS);
+      double Snm_mlambda = vTecLayer._S[n][m] * sin(m * _lonS);
+      vtec += (Snm_mlambda + Cnm_mlambda) * pnm;
+    }
+  }
+
+  if (vtec < 0.0) {
+    return 0.0;
+  }
+  return vtec;
+}
+
+void t_iono::piercePoint(double layerHeight, double epoch, const double* geocSta,
+    double sphEle, double sphAzi) {
+
+  double q = (t_CST::rgeoc + geocSta[2]) / (t_CST::rgeoc + layerHeight);
+
+  _psiPP = M_PI / 2 - sphEle - asin(q * cos(sphEle));
+
+  _phiPP = asin(sin(geocSta[0]) * cos(_psiPP) + cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi));
+
+  if (( (geocSta[0]*180.0/M_PI > 0) && (  tan(_psiPP) * cos(sphAzi)  > tan(M_PI/2 - geocSta[0])) )  ||
+      ( (geocSta[0]*180.0/M_PI < 0) && (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI/2 + geocSta[0])) ))  {
+    _lambdaPP = geocSta[1] + M_PI - asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
+  } else {
+    _lambdaPP = geocSta[1]        + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
+  }
+
+  _lonS = fmod((_lambdaPP + (epoch - 50400) * (M_PI / 43200)), 2*M_PI);
+
+  return;
+}
+
Index: trunk/BNC/src/pppModel.h
===================================================================
--- trunk/BNC/src/pppModel.h	(revision 7245)
+++ trunk/BNC/src/pppModel.h	(revision 7246)
@@ -7,4 +7,5 @@
 #include "t_prn.h"
 #include "satObs.h"
+#include "bncutils.h"
 
 namespace BNC_PPP {
@@ -63,11 +64,15 @@
   t_iono();
   ~t_iono();
-  void setTecData(const t_vTec* vTec);
-  double vtec();
+  double stec(const t_vTec* vTec, double signalPropagationTime,
+      const ColumnVector& rSat, const bncTime& epochTime,
+      const ColumnVector& xyzSta);
  private:
-  t_vTec*   _vTec;
-  int       _layer;
-  double    _layerHeight;
-  double    _roverHeight;
+  double vtecSingleLayerContribution(const t_vTecLayer& vTecLayer);
+  void piercePoint(double layerHeight, double epoch, const double* geocSta,
+      double sphEle, double sphAzi);
+  double _psiPP;
+  double _phiPP;
+  double _lambdaPP;
+  double _lonS;