Index: trunk/BNC/src/orbComp/sp3Comp.cpp
===================================================================
--- trunk/BNC/src/orbComp/sp3Comp.cpp	(revision 11023)
+++ trunk/BNC/src/orbComp/sp3Comp.cpp	(revision 11076)
@@ -231,4 +231,8 @@
         stat[prn.toString()]._offset = xx[index];
         if (dr.find(prn) != dr.end()){
+          // signal-in-space difference: radial minus clock, only the
+          // epoch-wise clock offset (datum) removed, i.e. including the
+          // satellite-specific clock offset
+          epochs[ie]->_sis[prn] = dr.find(prn)->second[0] - (dc[prn] + xx[index]);
           dcRed[prn] = dc[prn] - dr.find(prn)->second[0]; // clock minus radial component
           stat[prn.toString()]._dcRedMean += dcRed[prn];
@@ -368,6 +372,6 @@
     out << "!\n! Clock residuals and orbit differences in [m]\n"
            "! ----------------------------------------------------------------------------\n";
-    out << "!\n!  Epoch                PRN  radial   along   out        clk    clkRed   iPRN"
-            "\n! ----------------------------------------------------------------------------\n";
+    out << "!\n!  Epoch                PRN  radial   along   out        clk    clkRed   iPRN      SIS"
+            "\n! -------------------------------------------------------------------------------------\n";
   }
   for (unsigned ii = 0; ii < epochs.size(); ii++) {
@@ -376,4 +380,5 @@
     const map<t_prn, double>&       dc = epochs[ii]->_dc;
     const map<t_prn, double>&    dcRed = epochs[ii]->_dcRed;
+    const map<t_prn, double>&      sis = epochs[ii]->_sis;
     for (map<t_prn, ColumnVector>::const_iterator it = dr.begin(); it != dr.end(); it++) {
       const t_prn&  prn = it->first;
@@ -413,5 +418,20 @@
         }
         if (!_summaryOnly) {
-          out << "    " << setw(2) << int(prn) << "\n";
+          out << "    " << setw(2) << int(prn);
+        }
+        if (sis.find(prn) != sis.end()) {
+          double sisVal = sis.find(prn)->second;
+          stat[prn.toString()]._sisSum += sisVal;
+          stat[all.str()     ]._sisSum += sisVal;
+          stat[prn.toString()]._sisSqr += sisVal * sisVal;
+          stat[all.str()     ]._sisSqr += sisVal * sisVal;
+          stat[prn.toString()]._ns     += 1;
+          stat[all.str()     ]._ns     += 1;
+          if (!_summaryOnly) {
+            out << "  " << setw(7) << setprecision(4) << sisVal;
+          }
+        }
+        if (!_summaryOnly) {
+          out << "\n";
         }
       }
@@ -424,7 +444,16 @@
   out << "!\n! Summary";
   out << "\n! -----------------------------------------------------------------------------------------------------------------\n";
-  out << "!\n!     PRN   radialRMS   alongRMS    outRMS     3DRMS      nOrb    clkRMS    clkRedRMS   clkRedSig   nClk    Offset "
-          "\n!              [mm]       [mm]        [mm]      [mm]      [-]      [ns]        [ns]        [ns]      [-]     [ns]  "
-          "\n! -----------------------------------------------------------------------------------------------------------------\n";
+  out << "!\n!     PRN   radialRMS   alongRMS    outRMS     3DRMS      nOrb    clkRMS    clkRedRMS   clkRedSig   nClk    Offset     SISRMS    SISSig"
+          "\n!              [mm]       [mm]        [mm]      [mm]      [-]      [ns]        [ns]        [ns]      [-]     [ns]       [mm]      [mm]"
+          "\n! -------------------------------------------------------------------------------------------------------------------------------------\n";
+  // SIS standard deviation of a system: pooled about the satellite means
+  map<char, double> sisVarSys;
+  for (map<string, t_stat>::const_iterator it = stat.begin(); it != stat.end(); it++) {
+    const t_stat& st = it->second;
+    stringstream all; all << it->first[0] << 99;
+    if (it->first != all.str() && st._ns > 0) {
+      sisVarSys[it->first[0]] += st._sisSqr - st._sisSum * st._sisSum / st._ns;
+    }
+  }
   for (map<string, t_stat>::iterator it = stat.begin(); it != stat.end(); it++) {
     const string& prn  = it->first;
@@ -464,4 +493,13 @@
           out << setw( 9) << setprecision(2) << stat._offset / t_CST::c * 1e9;
         }
+        else {
+          out << setw( 9) << " ";
+        }
+      }
+      if (stat._ns > 0) {
+        double var = (prn == all.str()) ? sisVarSys[sys] / stat._ns
+                                        : stat._sisSqr / stat._ns - pow(stat._sisSum / stat._ns, 2);
+        out << setw(11) << setprecision(1) << sqrt(stat._sisSqr / stat._ns) * 1e3 << ' '
+            << setw( 9) << setprecision(1) << sqrt(max(0.0, var)) * 1e3;
       }
       out << "\n";