/********************************************************/
/*              Planets Module for SkyView              */
/*                                                      */
/*                  (c)1992 N P Hawkes                  */
/*                                                      */
/*                   Displays Planets                   */
/********************************************************/

#include "menu.h"
#include "dbox.h"
#include "bbc.h"
#include "wimpt.h"
#include "dbox.h"
#include "sprite.h"
#include "res.h"
#include "resspr.h"
#include "werr.h"

#include "sv_header.h"
#include "planets.h"
#include "datime.h"
#include "radec.h"
#include "ecliptic.h"
#include "riset.h"

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

/********************************************************/
/*                     Constants                        */
/********************************************************/
#define DISP_NAME "Planets"   /* Entry in Display menu. */
#define SEL_NAME  "Planet"    /* Entry in Select menu.  */

#define MAX_PLANS 10          /* Max no. of planets.    */
#define FILE_NAME "PlanData"  /* Name of data file.     */

#define NAME_LEN  15          /* Length of name field.  */
#define NCOEFFS    4          /* # of coeffs in polynoms*/

#define SPR_NAME "planet"     /* Name of sprite.        */
#define SPR_X0 -12            /* Bounding box of        */
#define SPR_Y0 -12            /* sprite, OS units,      */
#define SPR_X1  14            /* relative to            */
#define SPR_Y1  14            /* plotting position.     */
#define SPR_COL  5            /* Sprite x offset.       */
#define SPR_ROW  4            /* Sprite y offset.       */
#define SPR_MODE MODE_20      /* Sprite mode.           */

/*  Set effective altitude of horizon for calculations  */
/* of rising and setting.                               */
#define HORALT (REAL)0.0088

/* Convergence criteria for times of phenomena:         */
#define MAX_ITER 5
#define CONVERGE 5

/********************************************************/
/*                New types of variable.                */
/********************************************************/
/* Structure describing one planet.                     */
typedef struct {
  char name[NAME_LEN+1];  /* Name of planet.            */
  REAL meanl[NCOEFFS];    /* Coeffs for mean longitude. */
  REAL lperi[NCOEFFS];    /*   " for longit. of perihel.*/
  REAL eccen[NCOEFFS];    /*   " for eccentricity.      */
  REAL inclin[NCOEFFS];   /*   " for inclination.       */
  REAL lanode[NCOEFFS];   /*   " for long. of asc. node.*/
  REAL semimaj;           /* Semi-major axis.           */
  REAL angdiam;           /* Angular diameter at 1 AU.  */
  REAL stanmag;           /* Standard magnitude.        */
  REAL sdist;             /* Apparent dist from Sun.    */
  REAL edist;             /* Apparent dist from Earth.  */
  REAL phase;             /* Phase (1 = full).          */
  REAL ra;                /* Right Ascension of planet. */
  REAL dec;               /* Declination of planet.     */
  int  horiz_id;          /* id in Horiz window.        */
  int  vert_id;           /* id in Vert window.         */
} planstr;

/* Structure containing ecliptic rectangular            */
/* coordinates and velocity components.                 */
typedef struct {
  double x;
  double y;
  double z;
  double xdot;
  double ydot;
  double zdot;
} rect_coords;

/********************************************************/
/*                  Global Variables                    */
/********************************************************/
static BOOL display_flag = TRUE;   /* 'Enabled' flag.   */
static BOOL data_valid   = FALSE;  /* Validity of list. */
static int moduleid;               /* Module ID.        */

static FILE *fileptr;              /* Ptr to data file. */
static char fmt[256];              /* Format string.    */
static planstr plan_data[MAX_PLANS];/* Main data store. */
static int plan_count = 0;         /* No. of planets.   */

static menu plan_menu;    /* For selecting a planet.    */
static sprite_id spr_id;  /* ID of planet sprite.       */

static double T_now;      /* Centuries from 0/1/1900.   */

/********************************************************/
/*                 Function Prototypes                  */
/********************************************************/
static void plans_buildfn(void);
static void plans_selectfn(selectfn_reasoncode reason);
static void selection_details(void);
static BOOL cul_calc(int id, REAL altmin,
                 REAL *altptr, REAL *azimptr, int *hourptr, int *minptr);
static BOOL rs_calc(int id, BOOL whichcalc,
                 REAL *altptr, REAL *azimptr, int *hourptr, int *minptr);
static BOOL plans_displayfn(BOOL *enabptr);
static os_error *plans_plotfn(int x, int y, int id);
static void plans_infofn(int id);
static BOOL read_plandata(planstr *pptr);
static char *trim(char *str);
static void plan_altaz_etc(int id, rect_coords *earthptr,
                           REAL *altptr, REAL *azimptr);
static void earth_coords(double T, rect_coords *earthptr);
static void plan_coords(int id, double T, rect_coords *planptr);
static double ang_element(REAL coeff[NCOEFFS], double T);
static double ang_elemen1(REAL coeff[NCOEFFS], double T);
static double element(REAL coeff[NCOEFFS], double T);
static void plan_geo_elatlong(rect_coords *planptr, rect_coords *earthptr,
                              REAL *elatitptr, REAL *elongitptr);
static void plan_geo_elatlong_etc(
                              rect_coords *planptr, rect_coords *earthptr, 
                              REAL *elatitptr, REAL *elongitptr,
                              REAL *sdistptr,  REAL *edistptr,
                              REAL *phaseptr);
static void plan_radec(int id, double jul, REAL *raptr, REAL *decptr);
static void plan_radecfn(int id, observerstr *ob_ptr, int hour, int min,
                         REAL *raptr, REAL *decptr);
static void plan_altaz_any(int id, observerstr *ob_ptr, int hour, int min,
                           REAL *altptr, REAL *azimptr);

/********************************************************/
/*                Initialisation Function               */
/********************************************************/
BOOL plans_initfn(int moduleno, modulestr *plans)
{
  int i;
  sprite_id nametype_id;
  os_error *errptr;

/* Record the module number of this module.             */
  moduleid = moduleno;

/* Open data file.                                      */
  fileptr = res_openfile(FILE_NAME, "r");
  if (fileptr==NULL) return FALSE;

/* Build format string for reading name of planet.      */
  sprintf(fmt, " %%%ic", NAME_LEN);

/* Read data, one planet at a time.  Stop on first      */
/* error (probably EOF).  Stay within array bounds.     */
  while (plan_count < MAX_PLANS && \
         read_plandata(&plan_data[plan_count]))
    plan_count++;

/* Close file.  Fatal error if it won't close.          */
  if (fclose(fileptr) != 0)
    werr(FATAL, "Planets Error: Can't close Data file");

/* Quit if no items have been read.                     */
  if (plan_count == 0) return FALSE;

/* Build menu for Select function.                      */
/* Use menu_new for first item, menu_extend thereafter. */
  plan_menu = menu_new(SEL_NAME ":", plan_data[0].name);
  if (plan_menu == NULL) return FALSE;
  for (i=1; i<plan_count; i++)
    menu_extend(plan_menu, plan_data[i].name);

/* Fill in the fields of the modulestr.                 */
  plans->buildfn  = plans_buildfn;
  plans->selectfn = plans_selectfn;
  plans->dispfn   = plans_displayfn;
  plans->infofn   = plans_infofn;
  plans->initial  = display_flag;
  plans->display_entry = DISP_NAME;
  plans->display_menu  = NULL;
  plans->select_entry  = SEL_NAME;
  plans->select_menu   = plan_menu;

/* Set up the sprite_id structure which identifies the  */
/* planet sprite.                                       */
/* First set up a name-type identifier.                 */
  nametype_id.s.name = SPR_NAME;
  nametype_id.tag    = sprite_id_name;
/* Then get the address of the sprite, and put it in    */
/* the global sprite identifier spr_id.                 */
  errptr = sprite_select_rp(resspr_area(), &nametype_id, &spr_id.s.addr);
  if (errptr != NULL) return FALSE;
/* Finally tag the identifier spr_id as address-type.  */
  spr_id.tag = sprite_id_addr;

  return TRUE;
}

/*------------------------------------------------------*/
/*     Function to read in the data for one planet.     */
/* Data consist of:                                     */
/*   Planet's name,                                     */
/*   Coeffs for polynomials for five orbital elements,  */
/*   Semi-major axis of orbit,                          */
/*   Angular diameter at 1 AU,                          */
/*   Standard magnitude.                                */
/* All must be present.  Function returns TRUE if all   */
/* data read OK, FALSE otherwise.                       */
/* A suitable source of data is Duffett-Smith p133.     */
/*------------------------------------------------------*/
static BOOL read_plandata(planstr *pptr)
{
/* Temp storage for polynomial coefficients:            */
  float c0, c1, c2, c3;
/* Temp storage for other floating-point values:        */
  float f0, f1, f2;

/* Read name.                                           */
  if ( fscanf(fileptr, fmt, pptr->name) != 1 )
    return FALSE;

/* Put a terminating \0 into the name string.  Then     */
/* trim trailing blanks.                                */
  pptr->name[NAME_LEN] = '\0';
  trim(pptr->name);

/* Read coeffs for mean longitude.                      */
  if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
    return FALSE;
  pptr->meanl[0] = CONV * (REAL)c0;
  pptr->meanl[1] =        (REAL)c1;
  pptr->meanl[2] = CONV * (REAL)c2;
  pptr->meanl[3] = CONV * (REAL)c3;

/* Read coeffs for longitude of perihelion.             */
  if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
    return FALSE;
  pptr->lperi[0] = CONV * (REAL)c0;
  pptr->lperi[1] = CONV * (REAL)c1;
  pptr->lperi[2] = CONV * (REAL)c2;
  pptr->lperi[3] = CONV * (REAL)c3;

/* Read coeffs for eccentricity.                        */
  if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
    return FALSE;
  pptr->eccen[0] = (REAL)c0;
  pptr->eccen[1] = (REAL)c1;
  pptr->eccen[2] = (REAL)c2;
  pptr->eccen[3] = (REAL)c3;

/* Read coeffs for inclination.                         */
  if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
    return FALSE;
  pptr->inclin[0] = CONV * (REAL)c0;
  pptr->inclin[1] = CONV * (REAL)c1;
  pptr->inclin[2] = CONV * (REAL)c2;
  pptr->inclin[3] = CONV * (REAL)c3;

/* Read coeffs for longitude of ascending node.         */
  if ( fscanf(fileptr, "%f %f %f %f", &c0, &c1, &c2, &c3) != 4 )
    return FALSE;
  pptr->lanode[0] = CONV * (REAL)c0;
  pptr->lanode[1] = CONV * (REAL)c1;
  pptr->lanode[2] = CONV * (REAL)c2;
  pptr->lanode[3] = CONV * (REAL)c3;

/* Read semimajor axis, angular diam & standard magn.   */
  if ( fscanf(fileptr, "%f %f %f", &f0, &f1, &f2) != 3 )
    return FALSE;
  pptr->semimaj = (REAL)f0;
  pptr->angdiam = (REAL)f1;
  pptr->stanmag = (REAL)f2;

  return TRUE;
}

/*------------------------------------------------------*/
/*  Function which trims trailing blanks from a string. */
/*------------------------------------------------------*/
static char *trim(char *str)
{
  char *ptr;

  for (ptr = str+strlen(str); ptr>str && *(ptr-1)==' '; ptr--)
    ;
  *ptr = '\0';

  return str;
}
/********************************************************/
/*                 List-Building Function               */
/********************************************************/
static void plans_buildfn(void)
{
  int i;                /* Planet ID.                   */
  plotobj planet;       /* Plotobj under construction.  */
  rect_coords earth;    /* Current coords of Earth.     */

  /* Bounding box of plotting symbol, relative to       */
  /* position of symbol:                                */
  static wimp_box size = {SPR_X0, SPR_Y0, SPR_X1, SPR_Y1};

/* If module is disabled, do not rebuild list.  Instead,*/
/* set validity flag to FALSE and return.               */
  if (!display_flag)
  {
    data_valid = FALSE;
    return;
  }

/* Calculate current heliocentric ecliptic rectangular  */
/* coords of Earth.                                     */
  T_now = ob_data.jul/36525.0;
  earth_coords(T_now, &earth);

/* For each planet in turn:                             */
  for (i=0; i<plan_count; i++)
  {
    /* Build plotobj, and record interesting planetary  */
    /* data for possible use by info function.          */
    planet.id = i;
    planet.module = moduleid;
    planet.plotfn  = plans_plotfn;
    plan_altaz_etc(i, &earth, 
                   &planet.alt, &planet.azim);
    planet.size = size;

    /* Offer this to the main program.                  */
    addobj(planet, &plan_data[i].horiz_id, &plan_data[i].vert_id);
  }

/* Set flag to show that data in list is now valid.     */
  data_valid = TRUE;

  return;
}

/*------------------------------------------------------*/
/*  Function to calculate current alt & azim (and other */
/*  data of interest) for planet No. 'id'.  T_now and   */
/*coords of Earth (pointed to by earthptr) must be valid*/
/*------------------------------------------------------*/
static void plan_altaz_etc(int id, rect_coords *earthptr,
                           REAL *altptr, REAL *azimptr)
{
  rect_coords planet;
  REAL elatit, elongit;
  REAL ra, dec;

/*Get planet's heliocentric rectangular ecliptic coords.*/
  plan_coords(id, T_now, &planet);

/* Derive geocentric ecliptic latit. & longit, plus     */
/* other planetary data of interest.                    */
  plan_geo_elatlong_etc(&planet, earthptr,
           &elatit, &elongit,
           &plan_data[id].sdist, &plan_data[id].edist, &plan_data[id].phase);

/* Convert ecliptic lat & long to RA & Dec.             */
  ecliptic_radec(elatit, elongit, T_now, &ra, &dec);

/* Convert RA and Dec to altitude & azimuth.            */
  radec_altaz(ra, dec, &ob_data, ob_data.sid, altptr, azimptr);

  plan_data[id].ra  = ra;
  plan_data[id].dec = dec;

  return;
}

/*------------------------------------------------------*/
/* Function to calculate the RA and Dec of planet 'id'  */
/* at the instant implied by 'jul'.                     */
/*------------------------------------------------------*/
static void plan_radec(int id, double jul, REAL *raptr, REAL *decptr)
{
  double T;
  REAL elatit, elongit;
  rect_coords earth, planet;

  T = jul/36525.0;

/* Get Earth's heliocentric rectangular ecliptic coords.*/
  earth_coords(T, &earth);

/* Get same for planet.                                 */
  plan_coords(id, T, &planet);

/* Derive geocentric ecliptic latit. and longit.        */
  plan_geo_elatlong(&planet, &earth, &elatit, &elongit);

/* Convert latit. and longit. to RA and Dec.            */
  ecliptic_radec(elatit, elongit, T, raptr, decptr);

  return;
}

/*------------------------------------------------------*/
/* Function to calculate the RA and Dec of planet 'id'  */
/* at a specified instant.  Same as plan_radec above,   */
/* but arguments conform to requirements of 'Riset' fns.*/
/*------------------------------------------------------*/
static void plan_radecfn(int id, observerstr *ob_ptr, int hour, int min,
                         REAL *raptr, REAL *decptr)
{
  plan_radec(id, datime_julian(ob_ptr, hour, min), raptr, decptr);

  return;
}

/*------------------------------------------------------*/
/* Function to calculate the alt and azim of planet 'id'*/
/* at any time on the day designated in the specified   */
/* observerstr (which also defines the location).       */
/*------------------------------------------------------*/
static void plan_altaz_any(int id, observerstr *ob_ptr, int hour, int min,
                             REAL *altptr, REAL *azimptr)
{
  double jul;      /* Days since noon GMT 0/1/1900.     */
  REAL   sid;      /* Corresponding local siderial time.*/
  REAL ra, dec;

  jul = datime_julian(ob_ptr, hour, min);
  sid = datime_siderial(ob_ptr, jul, hour, min);

/* Get RA and Dec of planet.                            */
  plan_radec(id, jul, &ra, &dec);

/* Convert to altitude and azimuth.                     */
  radec_altaz(ra, dec, ob_ptr, sid, altptr, azimptr);

  return;
}

/*------------------------------------------------------*/
/*  Function to find rectangular heliocentric ecliptic  */
/*     coords and velocity components for the Earth.    */
/* Assumes semimajor axis is 1.0 AU, and inclination is */
/*  zero.  Method based on that in NAO Technical Note   */
/*      No.46, Section 5 (Dec 1978 / Feb 1981), p5.     */
/*------------------------------------------------------*/
static void earth_coords(double T, rect_coords *ptr)
{
  double L;      /* Mean longitude.                     */
  double M;      /* Mean anomaly.                       */
  double ec;     /* Eccentricity.                       */
  double E;      /* Eccentric anomaly.                  */
  double WW;     /* Longitude of perihelion.            */
  double x, y, r, u, v;
  double K = 0.0172021;

/* Get elements of Sun's "orbit".                       */
  ecliptic_sunorbit(T, &L, &M, &ec);

/* Convert mean longit. of Sun to that of Earth.        */
  L += DblPI;

/* Get longitude of perihelion.                         */
  WW = L - M;

/* Solve Kepler's eqn for eccentric anomaly.            */
  ecliptic_kepler(M, ec, &E);

/* Get coords and velocity components.                  */
  x = cos(E) - ec;
  y = sqrt(1.0 - ec*ec) * sin(E);
  r = sqrt(x*x + y*y);
  u = -(K * sin(E))/r;
  v =  (K * sqrt(1.0 - ec*ec) * cos(E))/r;

  ptr->x    = x*cos(WW) - y*sin(WW);
  ptr->y    = x*sin(WW) + y*cos(WW);
  ptr->z    = 0.0;
  ptr->xdot = u*cos(WW) - v*sin(WW);
  ptr->ydot = u*sin(WW) + v*cos(WW);
  ptr->zdot = 0.0;

  return;
}
/*------------------------------------------------------*/
/*  Function to find rectangular heliocentric ecliptic  */
/*    coords and velocity components for planet 'id'.   */
/*  Method based on that in NAO Technical Note No.46,   */
/*        Section 5 (Dec 1978 / Feb 1981), p5.          */
/*------------------------------------------------------*/
static void plan_coords(int id, double T, rect_coords *ptr)
{
  double L;      /* Mean longitude.                     */
  double WW;     /* Longitude of perihelion.            */
  double ec;     /* Eccentricity.                       */
  double i;      /* Inclination.                        */
  double O;      /* Longitude of ascending node.        */
  double a;      /* Semi major axis.                    */
  double M;      /* Mean anomaly.                       */
  double w;      /* Argument of the perihelion.         */
  double E;      /* Eccentric anomaly.                  */
  double x,  y,  r,  u,  v;
  double x1, y1,     u1, v1;
  double cosw, sinw;
  double cosi, sini;
  double cosO, sinO;
  double K = 0.0172021;

/* Get elements of orbit at time T.                     */
  L = ang_elemen1(plan_data[id].meanl, T);
  WW= ang_element(plan_data[id].lperi, T);
  ec=     element(plan_data[id].eccen, T);
  i = ang_element(plan_data[id].inclin,T);
  O = ang_element(plan_data[id].lanode,T);
  a = (double)plan_data[id].semimaj;

  M = L - WW;
  w = WW - O;

/* Solve Kepler's eqn for eccentric anomaly.            */
  ecliptic_kepler(M, ec, &E);

/* Get coords and velocity components.                  */
  x = a * (cos(E) - ec);
  y = a * sqrt(1.0 - ec*ec) * sin(E);
  r = sqrt(x*x + y*y);
  u = -(K * sqrt(a) * sin(E))/r;
  v =  (K * sqrt(a*(1.0 - ec*ec)) * cos(E))/r;

  cosw = cos(w);  cosi = cos(i);  cosO = cos(O);
  sinw = sin(w);  sini = sin(i);  sinO = sin(O);

  x1 = x*cosw - y*sinw;
  y1 = x*sinw + y*cosw;
  u1 = u*cosw - v*sinw;
  v1 = u*sinw + v*cosw;

  ptr->x    = x1*cosO - y1*sinO*cosi;
  ptr->y    = x1*sinO + y1*cosO*cosi;
  ptr->z    = y1*sini;
  ptr->xdot = u1*cosO - v1*sinO*cosi;
  ptr->ydot = u1*sinO + v1*cosO*cosi;
  ptr->zdot = v1*sini;

  return;
}

/*------------------------------------------------------*/
/*   Function to calculate an angular orbital element   */
/*   from its polynomial coefficients.  Result is in    */
/*   radians and in the range - 2 pi < element < 2 pi.  */
/*------------------------------------------------------*/
static double ang_element(REAL coeff[NCOEFFS], double T)
{
  return fmod(element(coeff, T), DblTWO_PI);
}

/*------------------------------------------------------*/
/*   Function to calculate an orbital element from its  */
/*              polynomial coefficients.                */
/*------------------------------------------------------*/
static double element(REAL coeff[NCOEFFS], double T)
{
  double el;
  double T2 = T * T;

  el  = (double)coeff[0];
  el += (double)coeff[1] * T;
  el += (double)coeff[2] * T2;
  el += (double)coeff[3] * T2 * T;

  return el;
}

/*------------------------------------------------------*/
/* Function to calculate an orbital element from coeffs */
/*      intended for use with a modified T term.        */
/* See Duffett-Smith p.132.                             */
/*------------------------------------------------------*/
static double ang_elemen1(REAL coeff[NCOEFFS], double T)
{
  double B;
  double el;
  double T2 = T * T;

/* Subtract integer rotations from T term, for better   */
/* accuracy.                                            */
  el  = (double)coeff[1] * T;
  B   = DblTWO_PI*(el - floor(el));

  el  = (double)coeff[0];
  el += B;
  el += (double)coeff[2] * T2;
  el += (double)coeff[3] * T2 * T;

  return fmod(el, DblTWO_PI);
}

/*------------------------------------------------------*/
/* Function to calculate geocentric ecliptic lat & long */
/* for a planet, plus other planetary data of interest. */
/*            Allows for light travel time.             */
/*------------------------------------------------------*/
static void plan_geo_elatlong_etc(
                           rect_coords *planptr,
                           rect_coords *earthptr,
                           REAL *elatptr,  REAL *elongptr,
                           REAL *sdistptr, REAL *edistptr, REAL *phaseptr)
{
  double R;           /* True distance from Earth.      */
  double rho;         /* Apparent distance from Earth.  */
  double xp, yp, zp;  /* Apparent geo. rectang. coords. */
  double tau;         /* Light travel time.             */
  double elongit;     /* App. geocent. ecliptic long.   */
  double helongit;    /* App. heliocent. ecliptic long. */
  double phase;       /* App. phase of planet.          */
  double x, y, z;
  double xdot, ydot, zdot;

/* Get true distance from Earth.                        */
  x = planptr->x - earthptr->x;
  y = planptr->y - earthptr->y;
  z = planptr->z - earthptr->z;
  R = sqrt(x*x + y*y + z*z);

/* Get light travel time.                               */
  tau = R / 173.14;

/* Get apparent distance from Earth.                    */
  xdot = planptr->xdot - earthptr->xdot;
  ydot = planptr->ydot - earthptr->ydot;
  zdot = planptr->zdot - earthptr->zdot;
  xp = x - tau*xdot;
  yp = y - tau*ydot;
  zp = z - tau*zdot;
  rho = sqrt(xp*xp + yp*yp + zp*zp);
  *edistptr = (REAL)rho;

/* Get apparent geocentric ecliptic lat & long.         */
   elongit  = atan2(yp, xp);
  *elatptr  = (REAL)asin(zp / rho);
  *elongptr = (REAL)elongit;

/* Get apparent distance from Sun.                      */
  xp = planptr->x  -  tau * planptr->xdot;
  yp = planptr->y  -  tau * planptr->ydot;
  zp = planptr->z  -  tau * planptr->zdot;
  *sdistptr = (REAL)sqrt(xp*xp + yp*yp + zp*zp);

/* Get apparent heliocentric ecliptic longitude.        */
  helongit = atan2(yp, xp);

/* Get phase (see Duffett-Smith p.139).                 */
  phase = 0.5 * (1.0 + cos(elongit - helongit));
  if (phase < 0.0) phase = 0.0;
  if (phase > 1.0) phase = 1.0;
  *phaseptr = (REAL)phase;

  return;
}

/*------------------------------------------------------*/
/* Function to calculate geocentric ecliptic lat & long */
/* for a planet, allowing for light travel time.        */
/*------------------------------------------------------*/
static void plan_geo_elatlong(
                           rect_coords *planptr,
                           rect_coords *earthptr,
                           REAL *elatptr,  REAL *elongptr)
{
  double R;           /* True distance from Earth.      */
  double rho;         /* Apparent distance from Earth.  */
  double xp, yp, zp;  /* Apparent geo. rectang. coords. */
  double tau;         /* Light travel time.             */
  double elongit;     /* App. geocent. ecliptic long.   */
  double x, y, z;
  double xdot, ydot, zdot;

/* Get true distance from Earth.                        */
  x = planptr->x - earthptr->x;
  y = planptr->y - earthptr->y;
  z = planptr->z - earthptr->z;
  R = sqrt(x*x + y*y + z*z);

/* Get light travel time.                               */
  tau = R / 173.14;

/* Get apparent distance from Earth.                    */
  xdot = planptr->xdot - earthptr->xdot;
  ydot = planptr->ydot - earthptr->ydot;
  zdot = planptr->zdot - earthptr->zdot;
  xp = x - tau*xdot;
  yp = y - tau*ydot;
  zp = z - tau*zdot;
  rho = sqrt(xp*xp + yp*yp + zp*zp);

/* Get apparent geocentric ecliptic lat & long.         */
   elongit  = atan2(yp, xp);
  *elatptr  = (REAL)asin(zp / rho);
  *elongptr = (REAL)elongit;

  return;
}


/********************************************************/
/*              Object-Selecting Function               */
/********************************************************/
static void plans_selectfn(selectfn_reasoncode reason)
{
  int menu_hit;

  switch (reason)
  {
    case new_selection:
      /* Look at list of menu hits to find which        */
      /* planet is to be selected.                      */
      menu_hit = module_submenu_hits[0];

      /* If no hits on Planet menu, flag the failure    */
      /* and return.  Also return if hit is off the end.*/
      if (menu_hit <= 0 || menu_hit > plan_count)
      {
        selection.selected_OK = FALSE;
        return;
      }

      selection.selected_OK = TRUE;
      /* Fill in Selection details for requested        */
      /* planet.  Allow for menu entries being          */
      /* indexed from 1 instead of 0.                   */
      selection.id = menu_hit - 1;
      selection_details();
      break;

    case window_selection:
      /* A selection has been made by clicking on a     */
      /* planet symbol in a window.                     */
    case recalculate_data:
      /* Observer details have changed, and plotting    */
      /* lists have been rebuilt.                       */
      /* */
      /*In either case, the planet ID is in selection.id*/
      selection_details();
      break;

    case timeonly_recalculate:
      /* As recalculate_data, but only the time of day   */
      /* (and possibly the direction of view - this is   */
      /* of no interest) have changed.                   */
      selection.horiz_id    = plan_data[selection.id].horiz_id;
      selection.vert_id     = plan_data[selection.id].vert_id;
      /* Decide if planet can be seen now.               */
      selection.now      = selection.horiz_id != NOWHERE || \
                           selection.vert_id  != NOWHERE;
      break;

    case sel_info_request:
      /* Write info into text buffer.                   */
      /* Use the Info function.                         */
      plans_infofn(selection.id);
      break;

    default:
      /* Unknown option.  Do nothing.                   */
      break;

    }

  return;
}


/*------------------------------------------------------*/
/*  Function to fill in the details of rising, setting  */
/*                   and culminating.                   */
/*------------------------------------------------------*/
static void selection_details(void)
{
  REAL rise_alt;  /*Altitude at calculated rising time. */
  REAL set_alt;   /*Altitude at calculated setting time.*/

  selection.horiz_id = plan_data[selection.id].horiz_id;
  selection.vert_id  = plan_data[selection.id].vert_id;
  selection.now      = selection.horiz_id != NOWHERE || \
                       selection.vert_id  != NOWHERE;

  selection.culminating = cul_calc(selection.id, HORALT,
                              &selection.cul_alt,
                              &selection.cul_azim,
                              &selection.cul_hour,
                              &selection.cul_min);
  selection.rising      = rs_calc(selection.id, RISECALC,
                              &rise_alt,
                              &selection.rise_azim,
                              &selection.rise_hour,
                              &selection.rise_min);
  selection.setting     = rs_calc(selection.id, SETCALC,
                              &set_alt,
                              &selection.set_azim,
                              &selection.set_hour,
                              &selection.set_min);
  return;
}

/*------------------------------------------------------*/
/*      Function to derive details of culmination.      */
/* If planet culminates, fn finds alt, azim, hour & min */
/* and returns TRUE.  Returns FALSE if no culmination.  */
/*------------------------------------------------------*/
static BOOL cul_calc(int id, REAL altmin,
                REAL *altptr, REAL *azimptr, int *hourptr, int *minptr)
{
  BOOL ok;

/* Guess a time of day.  Take midday.                   */
  *hourptr = 12;
  *minptr  =  0;

/* Try to improve guessed time.                         */
  ok = riset_cul(
       &ob_data, MAX_ITER, CONVERGE, plan_radecfn, id, hourptr, minptr);

/* Give up if time did not converge.                    */
  if (!ok) return FALSE;

/* Check that planet would be higher than 'altmin'.     */
  plan_altaz_any(id, &ob_data, *hourptr, *minptr, altptr, azimptr);
  return (*altptr >= altmin);
}

/*------------------------------------------------------*/
/*   Function to derive details of rising or setting.   */
/* If phenomenon occurs, fn finds alt, azim, hour & min */
/*      and returns TRUE.  Returns FALSE otherwise.     */
/*------------------------------------------------------*/
static BOOL rs_calc(int id, BOOL whichcalc,
                    REAL *altptr, REAL *azimptr, int *hourptr, int *minptr)
{
  BOOL ok;

/* Guess a time of day.  Take midday.                   */
  *hourptr = 12;
  *minptr  =  0;

/* Try to improve this guess.                           */
  ok = riset_riset(&ob_data, MAX_ITER, CONVERGE, plan_radecfn, id,
                   whichcalc, HORALT, hourptr, minptr);

/* Give up if time did not converge to a valid value.   */
  if (!ok) return FALSE;

/* Check that planet is indeed visible at this time.       */
  plan_altaz_any(id, &ob_data, *hourptr, *minptr, altptr, azimptr);
  return (*altptr > ZERO);
}


/********************************************************/
/*                   Display Function                   */
/********************************************************/
static BOOL plans_displayfn(BOOL *enabptr)
{
/* Toggle Enable/Disable flag.                          */
  display_flag = !display_flag;

/* Inform main program of new status.                   */
  *enabptr = display_flag;

/* If module is now enabled, but plotting list is not   */
/* valid, call list-building function.                  */
  if (display_flag && !data_valid)
    plans_buildfn();

/* Windows will always need updating.                   */
  return TRUE;
}


/********************************************************/
/*                   Plotting Function                  */
/********************************************************/
static os_error *plans_plotfn(int x, int y, int id)
{
  return sv_plotsprite(&spr_id, SPR_MODE, x, y, SPR_COL, SPR_ROW);
}


/********************************************************/
/*                     Info Function                    */
/********************************************************/
static void plans_infofn(int id)
{
  char p_text[RADEC_TEXTLEN]; /*Text buff for RA & Dec.*/
  REAL app_diam;  /* Apparent angular diameter.        */
  REAL logarg;
  REAL app_mag;   /* Apparent magnitude (see           */
                  /* Duffett-Smith p.139).             */

/* Get RA and Dec in text form.                        */
  radec_text(plan_data[id].ra, plan_data[id].dec, p_text);

/* Get apparent angular diameter and magnitude.        */
  app_diam = plan_data[id].angdiam / plan_data[id].edist;
  logarg = (plan_data[id].edist * plan_data[id].sdist) /
           (REAL)sqrt((double)plan_data[id].phase);
  app_mag  = (REAL)5.0*(REAL)log10((double)logarg) + plan_data[id].stanmag;

  sprintf(infoptr,
  "%s%s\n%s%.1f\n%s%.1f%s\n%s%.2f%s%.1f%c\n%s",
    "Planet: ", plan_data[id].name,
    "Approx visual mag. ", (double)app_mag,
    "Distance ", (double)plan_data[id].edist, " AU",
    "Phase ", (double)plan_data[id].phase,
    ", Ang diam ", (double)app_diam, '"',
    p_text);

  return;
}
