/*------------------------------------------------------*/
/*                         Riset                        */
/*    Set of functions for use in the calculation of    */
/*   rising, setting and culmination times of objects   */
/*          with time-dependent RA and Dec.             */
/*------------------------------------------------------*/

#include <math.h>

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

/*------------------------------------------------------*/
/*   Set maximum declination, beyond which object is    */
/*   considered so close to the celestial pole that it  */
/*                   does not move.                     */
/*------------------------------------------------------*/
#define RISET_NOMOVE (REAL)1.57


/*------------------------------------------------------*/
/* Function which attempts to find time of culmination, */
/*               given an initial guess.                */
/*------------------------------------------------------*/
  BOOL riset_cul(observerstr *ob_ptr, int maxit, int converge,
                 radecfntype *radecfn, int id, int *hour2ptr, int *min2ptr)
{
  int hour1, min1;  /* Time1.                           */
  int diff;         /* Difference between time1 & time2.*/
  REAL ra, dec;     /* RA & Dec of object.              */
  int iteration=0;  /* Iteration No.                    */

/* Iterate to find time of culmination.                 */
  do
  {
    iteration +=1;
    hour1 = *hour2ptr;
    min1  = *min2ptr;
  /* Calculate RA and Dec of object at time1.           */
    radecfn(id, ob_ptr, hour1, min1, &ra, &dec);
  /* Get time2, when this bit of sky culminates.        */
    datime_time_today(ra, ob_ptr, hour2ptr, min2ptr);
    diff = 60 * (*hour2ptr - hour1) + *min2ptr - min1;
    if (diff < 0) diff = -diff;
  }
  while (iteration < maxit  &&  diff > converge);

/* Inform caller whether or not the time converged.     */
  return (diff <= converge);
}


/*------------------------------------------------------*/
/*   Function which attempts to find time of rising or  */
/*          setting, given an initial guess.            */
/*------------------------------------------------------*/
  BOOL riset_riset(observerstr *ob_ptr, int maxit, int converge,
                   radecfntype *radecfn, int id, BOOL rise, REAL horalt,
                   int *hour2ptr, int *min2ptr)
{
  int hour1, min1;  /* Time1.                           */
  int diff;         /* Difference between time1 & time2.*/
  REAL ra, dec;     /* RA & Dec of object.              */
  int iteration=0;  /* Iteration No.                    */
  double d_dec;     /* double version of declination.   */
  double coshrang;  /* Cos of hour angle at 'horizon'.  */
  double hourang;   /* Hour angle at 'horizon'.         */
  BOOL cos_valid;   /* Is coshrang a valid cosine?      */
  REAL lst;         /* Local siderial time @ rise/set.  */
  /* double version of observer's latitude:             */
  double d_latit = (double)ob_ptr->latit;


/* Iterate to find time of rising / setting.            */

  do
  {
    iteration +=1;
    hour1 = *hour2ptr;
    min1  = *min2ptr;

  /* Calculate RA and Dec of object at time1.           */
    radecfn(id, ob_ptr, hour1, min1, &ra, &dec);

  /* Give up if declination is too close to 90 degrees. */
    if (dec >= RISET_NOMOVE || dec <=-RISET_NOMOVE)  return FALSE;

  /* Get cosine of hour angle when patch of sky         */
  /* specified by ra and dec rises.                     */
    d_dec = (double)dec;
    coshrang = (sin(horalt) - sin(d_latit)*sin(d_dec)) /
               (cos(d_latit)*cos(d_dec));

  /* Convert cosine to angle.  Note that for rising,    */
  /* true hour angle is actually 2PI - hourang.         */
    if      (coshrang >  1.0)
    {
      if (iteration == maxit) return FALSE;
      hourang = 0.0;
      cos_valid = FALSE;
    }
    else if (coshrang < -1.0)
    {
      if (iteration == maxit) return FALSE;
      hourang = DblPI;
      cos_valid = FALSE;
    }
    else
    {
      hourang = acos(coshrang);
      cos_valid = TRUE;
    }

  /* Get local siderial time when patch of sky has this */
  /* hour angle.                                        */
    if (rise)
      lst = ra - (REAL)hourang;
    else
      lst = ra + (REAL)hourang;

  /* Get time2, corresponding to this lst.              */
    datime_time_today(lst, ob_ptr, hour2ptr, min2ptr);

  /* Calculate difference between time2 & time1.        */
    diff = 60 * (*hour2ptr - hour1) + *min2ptr - min1;
    if (diff < 0) diff = -diff;

  }
  while ( (diff > converge || !cos_valid) && iteration < maxit );

/* Inform caller if time did not converge or cosine of  */
/* hour angle was invalid.                              */
  return (diff <= converge && cos_valid);
}
