*************************************************************************
* SEP - given celestial coordinates for two stars (RA1, Dec1) and (RA2, *
* Dec2), calculate the angular separation rho, position angle, and      *
* components along RA and Dec axes "delta(X)" "delta(Y)". The equations *
* for position angle and separation are rigorous (not approximations),  *
* using the law of cosines for the angles of spherical triangles.  The  *
* equation for "delta(X)" and "delta(Y)" is an approximation.  PA is    *
* measured as North through East as is astronomical custom. Based on    *
* equations from Chapter 2 of "Spherical Astronomy" by Edgar W. Woolard *
* and Gerald M. Clemence (1966) Academic Press, New York.               *
* Eric Mamajek 8/10/2014                                                *
*************************************************************************
      program sep
      implicit double precision (a-z)
      parameter(pi        = 3.141592653589793238d0)
      parameter(rad       = 180.0000000000000000000d0/pi)
      parameter(thrsixhun = 3600.0000000000000000d0)
      parameter(thrsixty  = 360.00000000000000000d0)
      parameter(zero      =   0.00000000000000000d0)
      write(*,*)'Calculate separation, position angle'
      write(*,*)'RA1, Dec1, RA2, Dec2 [degrees]'
      read(*,*)ra1,dec1,ra2,dec2

      deltara = ra2 - ra1
* calculate POSITION ANGLE (THETA; measured N through E)
      x  = dsin(deltara/rad)*dcos(dec2/rad)
      y1 = dcos(dec1/rad)*dsin(dec2/rad)
      y2 = dsin(dec1/rad)*dcos(dec2/rad)*dcos(deltara/rad)
      y  = y1 - y2
      pa = rad*datan2(x,y)
      if (pa .lt. zero) then
         pa = pa + thrsixty
      endif 
      write(*,*)' Position Angle [deg] = ',pa

* calculate delta(X), delta(Y)
      delx = deltara*dcos(dec2/rad)
      dely = dec2 - dec1
      delxasec = delx*thrsixhun
      delyasec = dely*thrsixhun
      write(*,*)'         Delta(X) ["] = ',delxasec
      write(*,*)'         Delta(Y) ["] = ',delyasec
     
* calculate SEPARATION (RHO) 
      f = dsin(dec1/rad)*dsin(dec2/rad) 
      g = dcos(deltara/rad)*dcos(dec1/rad)*dcos(dec2/rad)
      rho = rad*dacos(f+g)
      rhoasec = rho*thrsixhun
      write(*,*)'Separation \rho [deg] = ',rho
      write(*,*)'Separation \rho   ["] = ',rhoasec

      end
