/********************************************************/
/*                         RADec                        */
/* Utility functions for objects with known RA and Dec. */
/********************************************************/

#include <math.h>
#include <stdio.h>

#include "os.h"
#include "coords.h"
#include "menu.h"
#include "sv_header.h"
#include "radec.h"
#include "datime.h"

/* Set declination limit, above which object is consid- */
/* ered so close to the celestial pole that it doesn't  */
/* move during the day:                                 */
#define OBJECT_NOMOVE (REAL)1.57

/*------------------ Private functions -----------------*/
static void phenom_details(
            REAL ra, REAL dec, observerstr *ob_ptr, REAL hourang,
            double *cosazptr, int *hourptr, int *minptr);

/*------------------- Global variables -----------------*/
static double horalt;  /* 'Altitude' of horizon.        */


/*--------------------- Functions ----------------------*/

/*------------------------------------------------------*/
/*  Function to set altitude which object must exceed   */
/*          in order to be regarded as visible.         */
/*------------------------------------------------------*/
  void radec_sethoriz(REAL altitude)
{
  horalt = (double)altitude;
  return;
}

/*------------------------------------------------------*/
/*  Convert from RA & Dec to Alt & Azim following the   */
/*  method of Duffett-Smith p39.  Angles in radians.    */
/*------------------------------------------------------*/
    void radec_altaz(REAL ra, REAL dec, observerstr *ob_ptr, REAL sid,
                     REAL *altptr, REAL *azimptr)
{
  static double altmax = 1.5700;
  static double h;
  static double sind, cosd, sinlat, coslat;
  static double sinalt, alt;
  static double cosazim, azim;

/* Calculate local hour angle of object.                */
  h = (double)sid - (double)ra;

/* Calculate sin and cos of object's declination and    */
/* observer's latitude.                                 */
  sind   = sin((double)dec);
  cosd   = cos((double)dec);
  sinlat = sin((double)ob_ptr->latit);
  coslat = cos((double)ob_ptr->latit);

/* Calculate altitude of object.                        */
  sinalt = sind*sinlat + cosd*coslat*cos(h);
  if (sinalt < -1.0) sinalt = -1.0;
  if (sinalt >  1.0) sinalt =  1.0;
  alt = asin(sinalt);

/* Calculate azimuth of object.                         */
/* Special case when object is very nearly overhead.    */
  if (alt > altmax)
  /* Arbitrarily set azimuth to zero.                   */
    azim = 0.0;
  else
  {
    cosazim = (sind - sinlat*sinalt)/(coslat*cos(alt));
    if (cosazim < -1.0) cosazim = -1.0;
    if (cosazim >  1.0) cosazim =  1.0;
    azim = acos(cosazim);
  }

  *altptr = (REAL)alt;
  *azimptr= (sin(h) < 0.0 ? (REAL)azim: (REAL)2.0*PI - (REAL)azim);

  return;
}


/*------------------------------------------------------*/
/*  Function to calculate whether the selected object   */
/*  rises, sets or culminates today. Angles in radians. */
/*------------------------------------------------------*/
    void radec_phenomena(REAL ra, REAL dec, observerstr *ob_ptr,
                         double *chp, BOOL *riset_ptr, BOOL *cul_ptr)
{
/* Max zenith distance for object to be visible at      */
/* culmination:                                         */
  REAL zd_vis = PIby2 - (REAL)horalt;

  double d_latit  = (double)ob_ptr->latit;
  double d_dec    = (double)dec;


/* If object is at or very near the celestial pole, it  */
/* does not move, and therefore does not rise, set or   */
/* culminate.                                           */
  if (dec >= OBJECT_NOMOVE  || \
      dec <=-OBJECT_NOMOVE)
  {
    *riset_ptr = FALSE;
    *cul_ptr   = FALSE;
    return;
  }

/*              *** Rising and Setting ***              */

/*    Object is deemed to 'rise' & 'set' when at an     */
/*         altitude of 'horalt' radians.                */

/* Calculate cos of hour angle as object crosses this   */
/* altitude.                                            */
  *chp = (sin(horalt) - sin(d_latit)*sin(d_dec)) / \
                         (cos(d_latit)*cos(d_dec));

/* If *chp is not a valid cosine, object does not rise  */
/* or set.                                              */
  if (*chp >=  1.0 || *chp <= -1.0)
    *riset_ptr = FALSE;
  else
    *riset_ptr = TRUE;


/*                 *** Culmination ***                  */

/* Object is visible at culmination if less than zd_vis */
/* radians from the zenith when local hour angle of     */
/* Aries = rt ascension of object.                      */
  *cul_ptr = \
    dec > ob_ptr->latit - zd_vis && \
    dec < ob_ptr->latit + zd_vis;

  return;

}


/*------------------------------------------------------*/
/*    Function to calculate details of culmination.     */
/*               All angles in radians.                 */
/*------------------------------------------------------*/
    void radec_cul_details(REAL ra, REAL dec, observerstr *ob_ptr,
                           REAL *al, REAL *az, int *hourptr, int *minptr)
{
  REAL angdiff = dec - ob_ptr->latit;

/* Direction to look to see object at culmination.      */
  if (angdiff >= (REAL)0.0)
  {
    *al = PIby2 - angdiff;
    *az = (REAL)0.0;
  }
  else
  {
    *al = PIby2 + angdiff;
    *az = PI;
  }

/*                 Time of culmination.                 */
/* Object culminates when local siderial time = ra.     */

/* Calculate (one) civil time today at which the local  */
/* siderial time has the required value.                */
  datime_time_today(ra, ob_ptr, hourptr, minptr);

  return;
}


/*------------------------------------------------------*/
/*      Function to calculate details of rising.        */
/*               All angles in radians.                 */
/*------------------------------------------------------*/
    void radec_rise_details(REAL ra, REAL dec, observerstr *ob_ptr,
                            double ch, REAL *az, int *hourptr, int *minptr)
{
  double cosazim;      /* Cosine of azimuth of rising.  */

/* Get hour angle of object on rising.                  */
  REAL hourang = -(REAL)acos(ch);

/* Get hour, min and azimuth of rising.                 */
  phenom_details(ra, dec, ob_ptr, hourang, &cosazim, hourptr, minptr);
  *az = (REAL)acos(cosazim);

  return;
}


/*------------------------------------------------------*/
/*      Function to calculate details of setting.       */
/*               All angles in radians.                 */
/*------------------------------------------------------*/
    void radec_set_details(REAL ra, REAL dec, observerstr *ob_ptr,
                           double ch, REAL *az, int *hourptr, int *minptr)
{
  double cosazim;  /* Cosine of azimuth of setting.     */

/* Get hour angle of object on setting.                 */
  REAL hourang = (REAL)acos(ch);

/* Get hour, min and azimuth of setting.                */
  phenom_details(ra, dec, ob_ptr, hourang, &cosazim, hourptr, minptr);
  *az = TWO_PI - (REAL)acos(cosazim);

  return;
}

/*------------------------------------------------------*/
/* Function to calculate details of rising and setting. */
/* Based on function in Duffett-Smith, p106.            */
/*------------------------------------------------------*/
static void phenom_details(
            REAL ra, REAL dec, observerstr *ob_ptr, REAL hourang,
            double *cosazptr, int *hourptr, int *minptr)
{
  double d_latit  = (double)ob_ptr->latit;
  double d_dec    = (double)dec;
  REAL lst;        /* Local siderial time of phenomenon.*/

/* Get local siderial time of phenomenon.               */
  lst = ra + hourang;

/* Find civil time corresponding to this siderial time. */
  datime_time_today(lst, ob_ptr, hourptr, minptr);

/* Find cosine of azimuth of phenomenon.                */
  *cosazptr = (sin(d_dec) - sin(d_latit)*sin(horalt)) / \
                           (cos(d_latit)*cos(horalt));

/* Make sure *cosazptr is a valid cosine.               */
  if (*cosazptr >  1.0) *cosazptr =  1.0;
  if (*cosazptr < -1.0) *cosazptr = -1.0;

  return;
}


/*------------------------------------------------------*/
/*   Function to build a string giving RA and Dec in    */
/*          text form in conventional units.            */
/* String pointed to by textptr is assumed long enough. */
/*------------------------------------------------------*/
  void radec_text(REAL ra, REAL dec, char *textptr)
{
  REAL ram;    /* Right Ascension in minutes.          */
  REAL decm;   /* Declination in minutes.              */
  int  hrai;   /* Hours of RA (integer).               */
  int  mrai;   /* Minutes of RA (integer).             */
  int  ddeci;  /* Degrees of Dec (integer).            */
  int  mdeci;  /* Minutes of Dec (integer).            */
  int ns;      /* North or South declination.          */

/* Get RA in conventional units.  Round to nearest min.*/
  ram  = ra * (REAL)720.0 / PI + (REAL)0.5;
  if (ram < ZERO) ram = ZERO;
  hrai = (int)ram / 60;
  mrai = (int)ram % 60;
  if (hrai > 23) {hrai = 0;
                  mrai = 0;}

/* Get Dec in conventional units. Round to nearest min.*/
  decm = dec * (REAL)60.0 / CONV;
  if (decm < ZERO)
  {
    decm = -decm;
    if (decm < ZERO) decm = ZERO;
    ns = 'S';
  }
  else
    ns = 'N';
  decm += (REAL)0.5;
  ddeci = (int)decm / 60;
  mdeci = (int)decm % 60;
  if (ddeci > 89) {ddeci = 90;
                   mdeci =  0;}

  sprintf(textptr, "RA %02ih %02im   Dec %02i %02i' %c",
          hrai, mrai, ddeci, mdeci, ns);

  return;
}
