PROGRAM    GAIA_POS
      IMPLICIT   NONE 
      INCLUDE   'astro_constants.i'
      CHARACTER  RA_GAIA_STR*15, DEC_GAIA_STR*15, RA_VLBI_STR*15, DEC_VLBI_STR*15
      REAL*8     RA_GAIA, DEC_GAIA, RA_VLBI, DEC_VLBI, ARC, PHI
      REAL*8     D_RA_GAIA, D_DEC_GAIA, D_RA_VLBI, D_DEC_VLBI, D_ARC, D_PHI
      REAL*8,    EXTERNAL :: ARC_LEN_AD, ATAN_CS
      INTEGER*4  IUER
!
      RA_GAIA_STR  = '12_15_55.598877'
      DEC_GAIA_STR = '+34_48_15.20180'
      D_RA_GAIA    = 0.43/RAD__TO__MAS
      D_DEC_GAIA   = 0.45/RAD__TO__MAS
      RA_VLBI_STR  = '12_15_55.598593'
      DEC_VLBI_STR = '+34_48_15.19900'
      D_RA_VLBI    = 0.049/RAD__TO__MAS
      D_DEC_VLBI   = 0.055/RAD__TO__MAS
!
      CALL HR_TAT ( RA_GAIA_STR,  RA_GAIA, IUER  )
      CALL GR_TAT ( DEC_GAIA_STR, DEC_GAIA, IUER )
!
      CALL HR_TAT ( RA_VLBI_STR,  RA_VLBI, IUER  )
      CALL GR_TAT ( DEC_VLBI_STR, DEC_VLBI, IUER )
!
      ARC = ARC_LEN_AD ( RA_GAIA, DEC_GAIA, RA_VLBI, DEC_VLBI )
      PHI = ATAN_CS    ( (DEC_GAIA - DEC_VLBI), (RA_GAIA - RA_VLBI)*DCOS(DEC_VLBI) )
      D_ARC = ( (DEC_GAIA - DEC_VLBI)*DSQRT(D_DEC_GAIA**2 + D_DEC_VLBI**2) + &
     &          (RA_GAIA  - RA_VLBI )*DSQRT(D_RA_GAIA**2  + D_RA_VLBI**2 )*DCOS(DEC_VLBI) )/ ARC
      D_PHI = ( (DEC_GAIA - DEC_VLBI)*DSQRT(D_RA_GAIA**2  + D_RA_VLBI**2  )*DCOS(DEC_VLBI)  + &
     &          (RA_GAIA  - RA_VLBI )*DSQRT(D_DEC_GAIA**2 + D_DEC_VLBI**2 )*DCOS(DEC_VLBI)    )/ ARC**2
      WRITE ( 6, 110 ) ARC*RAD__TO__MAS, D_ARC*RAD__TO__MAS, PHI/DEG__TO__RAD, D_PHI/DEG__TO__RAD
 110  FORMAT ( 'Gaia - VLBI: ', F6.2, ' -+ ', F6.2, ' mas  Phi: ', F7.1, ' -+ ', F6.1, ' deg' )
      END  !#!