anchor.diff:  
   C     If we want the result in km/s²
   C     TODO: find out the meaning of the following const, conversion?
  -      fcv = 1.0 / (3.086e+16)
  +      fcv = 1.0 / METERS_IN_PARSEC
         total_force = fcv * abs(total_force)
   
   C     The sign of z_coordinate will be
  
positive.code: C         TODO: find out the meaning of ecinf    
            ecinf = 0.5 * y(2) * y(2)
            call epot(y(1), epotf)
            coordinate_Zcylindr(wd_index) = y(1) / (3.086e+16)
            ww(wd_index) = y(2) + WOSUN
        end do
        end subroutine
  
negative.code:       ypla=ycar
        rpla=sqrt(xpla*xpla+ypla*ypla)
        zsig=z
        z=abs(z/METERS_IN_PARSEC)                  
        vh2=0.5*vh*vh
        ro2=ro*ro
        rc12=rc1*rc1
  