diff --git a/binarys.hpp b/binarys.hpp
index 606a187673468bea0fa66831fa6f58e51d688e7c..8d9019b8d1892ef8791d236b499e3730a4b904ee 100644
--- a/binarys.hpp
+++ b/binarys.hpp
@@ -41,7 +41,6 @@ Type solveKepler(Type eccentricity, Type meanAnomaly) {
   return eccentricAnomaly;
 }
 
-// Campbell parameters method
 template<class Type>
 vector<Type> relpos(Type period, Type timeperiastron, Type eccentricity, Type omega2, Type inclination, Type nodeangle, Type semiMajorAxis, Type massRatio,
                     Type time) {
@@ -52,7 +51,6 @@ vector<Type> relpos(Type period, Type timeperiastron, Type eccentricity, Type om
   Type sinOmg = sin(nodeangle);
   Type cosOmg = cos(nodeangle);
 
-  // StarSystem.computeOrbitComponent
   Type eccentricAnomaly = solveKepler(eccentricity, 2. * M_PI * (time - timeperiastron) / period);
 
   Type v = 2. * atan2(efac * sin(eccentricAnomaly / 2.), cos(eccentricAnomaly / 2.)); // true anomaly
@@ -60,9 +58,8 @@ vector<Type> relpos(Type period, Type timeperiastron, Type eccentricity, Type om
   Type oneminusecosE = 1 - eccentricity * cos(eccentricAnomaly);
   Type semiAmplitudeK2 = (AU_SI*1e-3*2*M_PI/DAYSEC *semiMajorAxis * sinI) /
     (period * sqrt(1 - eccentricity * eccentricity));
-  // OrbitalParams.computeOrbitPos
   Type r = semiMajorAxis * oneminusecosE;
-  // StarSystem.computeOrbitComponent
+
   Type ksi2 = r*(cos(Q)*sinOmg + sin(Q)*cosOmg*cosI);
   Type eta2 = r*(cos(Q)*cosOmg - sin(Q)*sinOmg*cosI);
   // distance offset (in A.U.) along the line of sight of the component wrt the barycenter
@@ -104,6 +101,8 @@ vector<Type> getIADSolution(vector<Type> Ares, vector<Type> dpra_nu, vector<Type
   return(X.vec());
 }
 
+// This is an implementation of the equations in The Hipparcos and Tycho Catalogues (ESA SP-1200), Volume 1, Section 1.5.5, 'Epoch Transformation: Rigorous Treatment'.
+// Adapted from L. Ruiz-Dern, itself adapted from L. Lindegren, itself adapted from D. Priou.
 template<class Type>
 vector<Type> epochpropaFull(Type ra, Type dec, Type plx, Type pmra, Type pmdec, Type RV, Type epoch1, Type epoch2) {
 
@@ -204,7 +203,5 @@ vector<Type> epochpropaFull(Type ra, Type dec, Type plx, Type pmra, Type pmdec,
   res[4] = pmd / MAS2RAD;
   res[5] = zeta / (res[2]/KAPPA) / MAS2RAD;
 
-  // std::cout << res(0) << " " << res(1) << " " << res(2) << " " << res(3) << " " << res(4) << " " << res[5] << std::endl;
-
   return(res);
 }