anchor.diff:       &        zkr,
        &        zkri,
        &        xc,
  -     &        xs
  +     &        xs,
  +     &        velocity_by_prop_motion
         double precision :: solarGalactocentricDistance,
        &                    ros,
        &                    cos_bgac,
  
positive.code:      &               radial_velocity
        common /lb/ lgac, bgac
  
  C     Calculating galactocentric coordinates (r,l,b)
        do wd_index = 1, numberOfWDs
  C         Galactic coordinate r (Kpc)
  C         TODO: find out the meaning of ros
            ros = solarGalactocentricDistance 
       &          * solarGalactocentricDistance
       &          + coordinate_R(wd_index) * coordinate_R(wd_index) 
  
negative.code:           cos_lgac = real(cos(lgac(wd_index)))
            sin_bgac = dsin(bgac(wd_index))
            cos_bgac = dcos(bgac(wd_index))
  
            zkr = KAPPA * rgac(wd_index) * PC_PER_KPC
            zkri = 1.0 / zkr
            velocity_by_prop_motion = 1.0 / (KAPPA * rgac(wd_index) 
       &                                     * PC_PER_KPC)
  
  C         TODO: find out why dividing by cos_bgac
            longitude_proper_motion(wd_index) = real(
       &        velocity_by_prop_motion
       &        * (-uu(wd_index) * sin_lgac / cos_bgac
       &           + vv(wd_index) * cos_lgac / cos_bgac))
            latitude_proper_motion(wd_index) = real(
       &        velocity_by_prop_motion
       &        * (-uu(wd_index) * cos_lgac * sin_bgac
       &           -vv(wd_index) * sin_bgac * sin_lgac
       &           + ww(wd_index) * cos_bgac))
            radial_velocity(wd_index) = real(
       &        uu(wd_index) * cos_bgac * cos_lgac
       &        + vv(wd_index) * cos_bgac * sin_lgac
       &        + ww(wd_index) * sin_bgac)
  C         FIXME: Looks like this formula is wrong
            properMotion(wd_index) = sqrt(
       &        longitude_proper_motion(wd_index) 
  