/********************************************************/
/*                   --  Ecliptic  --                   */
/*          Sundry utility functions concerning         */
/*            ecliptic latitude and longitude.          */
/********************************************************/
#include <math.h>

#include "os.h"
#include "coords.h"
#include "menu.h"
#include "sv_header.h"
#include "ecliptic.h"

/* Set (arbitrary) declination beyond which object is   */
/* considered to be at the celestial pole.              */
#define DEC_LIMIT 1.57

/* Set convergence criteria for solving Kepler's eqn.   */
#define MAX_ITER 10         /* Max number of iterations.*/
#define KEP_CONVERGE 1.0E-6 /* Max change in last iter. */

/*------------------------------------------------------*/
/* Function returning the obliquity of the ecliptic at  */
/*  the instant implied by T (No. of Julian centuries   */
/*   from noon GMT on 0/1/1900).  Result in radians.    */
/*------------------------------------------------------*/
  double ecliptic_obliq(double T)
{
  return ( 23.452 - 0.013 * T )*DblCONV;
}

/*------------------------------------------------------*/
/*   Function to calculate the elements of the sun's    */
/*     "orbit" at the instant defined by T (Julian      */
/*    centuries from 0/1/1900).  Angles in radians.     */
/*        Coeffs from P. Duffett-Smith p. 116.          */
/*------------------------------------------------------*/
  void ecliptic_sunorbit(double T, double *Lptr, double *Mptr, double *eptr)
{
  double T2 = T * T;

/* Mean latitude: */
  *Lptr = 2.7969668E2 + (3.600076892E4 * T) + (3.025E-4 * T2);
  *Lptr = fmod(*Lptr, 360.0)*DblCONV;

/* Mean anomaly:  */
  *Mptr = 3.5847583E2 + (3.599904975E4 * T) - (1.54E-4 * T2) -
          (3.3E-6 * T * T2);
  *Mptr = fmod(*Mptr, 360.0)*DblCONV;

/* Eccentricity:  */
  *eptr = 1.675104E-2 - (4.18E-5 * T) - (1.26E-7 * T2);

  return;
}

/*------------------------------------------------------*/
/* Function to solve Kepler's Equation.  Follows method */
/*  given in NAO Technical Note Section 5 (Dec 1978 &   */
/*    Feb 1981) p7.  Returns TRUE if E converged OK.    */
/*------------------------------------------------------*/
  BOOL ecliptic_kepler(double M, double e, double *Eptr)
{
  int iter_count = 0;
  double denom, newE;
  BOOL done;

  *Eptr = M + e * sin(M);
  do
  {
    iter_count += 1;
    denom = 1 - e * cos(*Eptr);
    if (denom == 0.0) return FALSE;
    newE = *Eptr + (M + e*sin(*Eptr) - *Eptr)/denom;
    done = fabs(*Eptr - newE) <= KEP_CONVERGE;
  } while (!done && iter_count < MAX_ITER);

  return (done);
}

/*------------------------------------------------------*/
/* Function to convert ecliptic latitude and longitude  */
/* to right ascension and declination (all in radians). */
/*------------------------------------------------------*/
  void ecliptic_radec(REAL elatit, REAL elongit, double T,
                      REAL *raptr, REAL *decptr)
{
  double beta;      /* Ecliptic latitude.               */
  double lambda;    /* Ecliptic longitude.              */
  double epsilon;   /* Obliquity of the ecliptic.       */
  double alpha;     /* Right ascension.                 */
  double delta;     /* Declination.                     */
  double x, y;
  double arg;

  beta   = (double)elatit;
  lambda = (double)elongit;

  epsilon = ecliptic_obliq(T);
  arg     = sin(beta)*cos(epsilon) + cos(beta)*sin(epsilon)*sin(lambda);
  if (arg < -1.0) arg = -1.0;
  if (arg >  1.0) arg =  1.0;
  delta   = asin(arg);

/* Special case when object is very close to celestial  */
/* poles.                                               */
  if (fabs(delta) > DEC_LIMIT)
  {
/* Set right ascension to an arbitrary value (0.0).     */
    *raptr  = ZERO;
    *decptr = (REAL)delta;
    return;
  }

  y = cos(beta)*sin(lambda)*cos(epsilon) - sin(beta)*sin(epsilon);
  x = cos(beta)*cos(lambda);
  alpha = atan2(y,x);
  if (alpha <  0.0      ) alpha += 2.0*DblPI;
  if (alpha >= 2.0*DblPI) alpha  = 0.0;

  *raptr  = (REAL)alpha;
  *decptr = (REAL)delta;
  return;
}
