#include <math.h>
#include "watdefs.h"
#include "afuncs.h"

/*
   calc_dist_and_posn_ang() takes two RA/dec points,  and computes their
angular separation and position angle.  If you don't want the PA,  pass
in NULL for the posn_ang argument.

   After a quick check to make sure p1 != p2,  the code uses a haversine
method to get the distance.  This handles small angles correctly.  The
"usual" formula for such a spherical triangle,

cos(a) = cos(b) * cos(c) + sin(b) * sin(c) * cos( alpha)

   has serious precision problems when cos(a) is close to +/-1 (a=0 or
180 degrees).  Instead,  the haversine formula mentioned in Meeus'
_Astronomical Algorithms_,  page 111,  is used:  (use a fixed-size
font to read the following equations:)

hav( dist) = hav( dec1 - dec2) + cos( dec1) * cos( dec2) * hav( delta_ra)

   A quick review of haversines:  it helps to know that if

                     2
x = hav( theta) = sin ( theta / 2) = (1 - cos( theta) / 2)

   then the "inverse haversine" is computed as

theta = inv_haversine( x) = 2 * arcsine( sqrt( x))

   Once it has the distance,  getting the PA is just about trivial,  and
is again done using haversines,  using the spherical triangle formula

   2                              cos( b - c) - cos( a)
sin ( alpha / 2) = hav( alpha) = -----------------------
                                    2 sin( b) sin( c)

   ...which,  after setting b = 90 - dec1, c = dist, a = 90 - dec2,
alpha = PA (from point 1 to point 2),  gives

            sin( dec1 + dist) - sin( dec2)
hav( PA) = --------------------------------
                 2 cos( dec1) sin( dist)

   In using this,  the code first makes sure that sin(b) * sin(c) isn't
zero (which could happen if the points are superimposed, making dist=0, or
if the starting point is at the pole,  dec1 = +/- 90... in either case,
the PA is undefined).  After doing the division,  it makes sure the result
is <= 1;  mathematically,  it never will be, but roundoff happens.

   The resulting PA lies between 0 and 180,  because the above formula
can't distinguish east from west.  So we do a check of sin( ra1 - ra2)
to fix that.

   (14 Apr 2003) Bill Owen kindly pointed out the following haversine
method of dealing with the PA.  As he points out,  it evades the quadrant
check;  perhaps more importantly,  I think it's resistant to the above
"roundoff happens" problem:

   "...As for position angle, the rigorous formula

  PA = atan2 ( sin (RA2-RA1) cos Dec2,
               cos Dec1 sin Dec2 - cos (RA2-RA1) sin Dec1 cos Dec2 )

is somewhat amenable to the haversine treatment too.  Apply the identity
  cos a = 1 - 2 sin^2 (a/2) = 1 - 2 hav a
to cos (RA2-RA1) in the denominator, and I get

  PA = atan2 ( sin (RA2-RA1) cos Dec2,
               sin (Dec2-Dec1) + 2 hav (RA2-RA1) sin Dec1 cos Dec2 )

which obviates the need for a quadrant check; atan2 does that for you.

In the limit, of course, the haversine term goes away, the sines in the
leading term get replaced by the angles, and we're left with

  PA ~ atan2 ( (RA2-RA1) cos Dec, (Dec2-Dec1) )

as we expect.

-- Bill"
*/

#define PI 3.141592653589793238462643383279502884197
#define HALF_PI (PI / 2.)

static double hav( const double angle)
{
   double rval = sin( angle / 2.);

   return( rval * rval);
}

static double inv_hav( const double hav)
{
   return( 2. * asin( sqrt( hav)));
}

int DLL_FUNC calc_dist_and_posn_ang( const double ra1, const double dec1,
                              const double ra2, const double dec2,
                              double *dist, double *posn_ang)
{
   double delta_ra = ra1 - ra2;
   int rval = -1;

   if( delta_ra || dec1 != dec2)
      {
      const double hav_dist =
             hav( dec1 - dec2) + cos( dec1) * cos( dec2) * hav( delta_ra);

      if( hav_dist >= 1.)         /* points are at 'antipodes'.  Roundoff  */
         {                        /* errors can make hav_dist > 1.         */
         *dist = PI;
         if( posn_ang)
            *posn_ang = 0.;
         }
      else
         {
         *dist = inv_hav( hav_dist);
         if( posn_ang)
            {
            double tval = cos( dec1) * sin( *dist);

            if( !tval)              /* dec1 is at a pole... undefined p.a. */
               *posn_ang = 0.;
            else
               {
               const double hav_pa = (sin( dec1 + *dist) - sin( dec2)) / (2. * tval);

               if( hav_pa > 1.)            /* again,  roundoff troubles */
                  *posn_ang = PI;
               else
                  *posn_ang = inv_hav( hav_pa);
               if( sin( delta_ra) < 0.)             /* Fixed 23 aug 94 to */
                  *posn_ang = PI + PI - *posn_ang;  /* run 90=E,  not W   */
               }
            }
         }
      rval = 0;
      }
   else           /* points are identical */
      {
      *dist = 0.;
      if( posn_ang)          /* posn angle is really undefined,  but let's  */
         *posn_ang = 0.;     /* zero it to evade "Not A Number" messages    */
      }
   return( rval);
}