anchor.diff:       &        + latitude_proper_motion(wd_index) 
        &          * latitude_proper_motion(wd_index))
   
  -C         Calculating right ascension and the declination      
  -C         Calculating the declination
  -C         TODO: find out the meaning of zz
  -          zz = SIN_DELTA_G * sin_bgac 
  -     &         + COS_DELTA_G * cos_bgac * dcos(THETA - lgac(wd_index))
  -          declination(wd_index) = dasin(zz)
  +C         More info at the link above
  +          declination(wd_index) = dasin(
  +     &        SIN_DELTA_G * sin_bgac + COS_DELTA_G * cos_bgac 
  +     &                                 * dcos(THETA - lgac(wd_index)))
   
  -C         Calculating the right ascension
  -C         TODO: find out the meaning of xs
  +C         TODO: give better names for xs and xc
  +C         These vars are used for conversion from galactic to equatorial 
  +C         coordinates. More info at the link above
             xs = real((cos_bgac * dsin(THETA - lgac(wd_index))) 
        &              / dcos(declination(wd_index)))
  -C         TODO: find out the meaning of xc
             xc =  real((COS_DELTA_G * sin_bgac 
        &                - SIN_DELTA_G * cos_bgac 
        &                  * cos(THETA - lgac(wd_index)))
        &               / cos(declination(wd_index)))
   
   C         Looking at the sign that corresponds to the right ascension
  +C         TODO: find out what is going on here
             if (xs >= 0.0) then
                 if (xc >= 0.0) then 
                     rightAscension(wd_index) = asin(xs) + ALPHA_G
  
positive.code: C     Distances in Kpc
  C     Velocities in Km/s
  C     Proper motions in  arcsec/yr
        implicit none
  
        integer, parameter :: MAX_STARS_COUNT = 6000000
  
negative.code:       double precision :: solarGalactocentricDistance,
       &                    ros,
       &                    cos_bgac,
       &                    sin_bgac,
       &                    zzx
        double precision :: coordinate_R(MAX_STARS_COUNT),
  