/********************************************************************
sweephe4.c
access structures and functions for ephemeris file ep4_
a fast precomputed ephemeris
*********************************************************************/
/* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland.  All rights reserved.

  License conditions
  ------------------

  This file is part of Swiss Ephemeris.

  Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND.  No author
  or distributor accepts any responsibility for the consequences of using it,
  or for whether it serves any particular purpose or works at all, unless he
  or she says so in writing.  

  Swiss Ephemeris is made available by its authors under a dual licensing
  system. The software developer, who uses any part of Swiss Ephemeris
  in his or her software, must choose between one of the two license models,
  which are
  a) GNU Affero General Public License (AGPL)
  b) Swiss Ephemeris Professional License

  The choice must be made before the software developer distributes software
  containing parts of Swiss Ephemeris to others, and before any public
  service using the developed software is activated.

  If the developer choses the AGPL software license, he or she must fulfill
  the conditions of that license, which includes the obligation to place his
  or her whole software project under the AGPL or a compatible license.
  See https://www.gnu.org/licenses/agpl-3.0.html

  If the developer choses the Swiss Ephemeris Professional license,
  he must follow the instructions as found in http://www.astro.com/swisseph/ 
  and purchase the Swiss Ephemeris Professional Edition from Astrodienst
  and sign the corresponding license contract.

  The License grants you the right to use, copy, modify and redistribute
  Swiss Ephemeris, but only under certain conditions described in the License.
  Among other things, the License requires that the copyright notices and
  this notice be preserved on all copies.

  Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl

  The authors of Swiss Ephemeris have no control or influence over any of
  the derived works, i.e. over software or services created by other
  programmers which use Swiss Ephemeris functions.

  The names of the authors or of the copyright holder (Astrodienst) must not
  be used for promoting any software, product or service which uses or contains
  the Swiss Ephemeris. This copyright notice is the ONLY place where the
  names of the authors can legally appear, except in cases where they have
  given special permission in writing.

  The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
  for promoting such software, products or services.
*/


# include "swephexp.h"
# include "sweephe4.h"
# include <string.h>

# define INVALID_BASE	2000000000L
# define EPBS	(2 * NDB)	/* buffer size is 20 days */
# define EP_MIN_IX	2	/* load buffer when index below this */
# define EP_MAX_IX	(EPBS - 4)	/* load buffer when index above this */

FILE *ephfp = NULL;
const int qod[EP_NP] = {5,5,5,5,5,3,3,3,3,3,3,5,3,3,3};

static void inpolq_l(int n, int o, double p, centisec *x,
		      centisec *axu, centisec *adxu);
static int inpolq(int n, int o, double p, double *x,
		      double *axu, double *adxu);
static int ephe4_unpack(int jdl, int pflag, centisec lon[][EPBS], int i0,
			 char *errs);
static int ephe4_unpack_d(int jdl, int pflag, double lon[][EPBS], int i0,
			 char *errs);
static char *my_makepath(char *d, char *s);

# ifdef INTEL_BYTE_ORDER
/********************************************************************/
void shortreorder (UCHAR *p, int n) 
			   /* p points to memory filled with 16-bit values; for
                           each of the values the seqeuence of the two bytes
                           has to be reversed, to translate HP-UX and VAX
			   ordering to MSDOS/Turboc ordering */
{
  int i;
  unsigned char c0;
  for (i = 0; i < n; i += 2, p += 2) {
    c0 = *p;
    *p = *(p+1);
    *(p + 1) = c0;
  }
}
# endif


/****************************************************
  read ephe file and return a pointer to normalized positions
  for planets specified by pflag at julday jdl;
  If the reading failes, NULL is returned and
  an error text of max. 79 char in errtext, except when errtext = NULL.
  If calc is used and succeeds, a message is put into errtext, otherwise
  errtext is empty.
  Attention: jd is an absolute Julian date
  flag bits defined in sweephe4.h, onl
****************************************************/
centisec *ephread(double jd, int plalist, int flag, char *errtext)
{
  static int jdbase =  INVALID_BASE;
  static int lastplalist = 0;
  static centisec lon[EP_NP][EPBS];	/* buffer for 20 days unpacked ephe */
  static centisec out[2 * EP_NP];	/* buffer for return longitude
					   and return speed */
  int p, pf;
  int ix, jdlong, iflagret;
  centisec clp;
  double jfract;
  double x[6];
  if (errtext != NULL)
    *errtext = '\0';
  if (plalist == 0)
    plalist = EP_ALL_BITS;		/* default: all logitudes, no speeds */
  /*
   * we must determine when to reload the lon buffer: if the contents do
   * not allow immediate interpolation or if the plalist selection has
   * changed since the last call.
   */
  if ((plalist & lastplalist) != plalist) {	/* new set is not contained in old */
    jdbase = INVALID_BASE;
  }
  lastplalist = plalist;
  jdlong = floor(jd - 0.5);
  ix = jdlong - jdbase;
  if (ix < EP_MIN_IX || ix >= EPBS) {	/* must reload full buffer */
    jdbase = ((jdlong - EP_MIN_IX) / NDB) * NDB;		/* new base */
    if (jdbase > jdlong - EP_MIN_IX) jdbase -= NDB;	/* fix bug for neg. */
    if (ephe4_unpack (jdbase, plalist, lon, 0, errtext) != OK)
      goto err_exit;
    if (ephe4_unpack (jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
      goto err_exit;
    ix = jdlong - jdbase;
  } else if (ix > EP_MAX_IX) {		/* must shift upper half down
					   and reload upper half of buffer */
    jdbase +=  NDB;		/* new base */
    for (p = 0; p < EP_NP; p++)
      memcpy (&lon[p][0], &lon[p][NDB], NDB * sizeof(centisec));
    if (ephe4_unpack (jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
      goto err_exit;
    ix = jdlong - jdbase;
  }
  jfract = jd - 0.5 - jdlong;        
  /*
   * we use the interpolator even for jfract = 0, because it delivers
   * the speed term. The computation overhead is unimportant
   * in any case.
   */
  for (p = 0, pf = 1; p < EP_NP; p++, pf = pf << 1)
    if ((plalist & pf) != 0) {
      inpolq_l((int) ix, qod[p], jfract, &(lon[p][0]), &(out[p]), &clp);
      if (p <= PLACALC_CHIRON) {	/* normalize all except ecl and nut */
	if (out[p] < 0)
	  out[p] += DEG360;
	else if (out[p] >= DEG360)
	  out[p] -= DEG360;
      }
#ifdef DEBUG
      fprintf(stderr,"ephread p=%d, lon=%.3lf\n", p, out[p] * CS2DEG);
#endif
      if (flag & EP_BIT_SPEED)
	out[p+EP_NP] = clp;
    }
  return out;
err_exit:
  jdbase = INVALID_BASE;
  lastplalist = 0;
  if ((flag & EP_BIT_MUST_USE_EPHE) == 0) {	/* try using calc */
    int sweflag = 0;
    char serr[AS_MAXCH];
    if (flag & EP_BIT_SPEED)
      sweflag = SEFLG_SPEED;
    if (errtext != NULL)
      sprintf(errtext,"ephread failed for jd=%f; used swe_calc().", jd);
    for (p = 0, pf = 1; p < PLACALC_CALC_N; p++, pf = pf << 1) {
      if ((plalist & pf) != 0) {
	if ((iflagret = swe_calc(jd, ephe_plac2swe(p), sweflag, x, serr)) != ERR) {
	  out[p] = swe_d2l(x[0] * DEG);
	  if (flag & EP_BIT_SPEED)
	    out[p + EP_NP] = swe_d2l(x[3] * DEG);
	  if (out[p] < 0)
	    out[p] += DEG360;
	  else if (out[p] >= DEG360)
	    out[p] -= DEG360;
	} else {
	  swe_close();
          if (errtext != NULL)
	    strcat(errtext," calc failed too.");
          return NULL;
	}
      }
    }
    if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
      swe_close();
      sprintf(errtext, "error in swe_calc() %s\n", serr);
      return NULL;
    }
    out[EP_ECL_INDEX] = swe_d2l(x[0] * DEG);	/* true ecliptic */
    out[EP_NUT_INDEX] = swe_d2l(x[2] * DEG);	/* nutation */
    out[EP_ECL_INDEX + EP_NP] = 0;
    out[EP_NUT_INDEX + EP_NP] = 0;
    return out;
  }
  return NULL;
}	/* ephread */

// same in double
double *dephread2(double jd, int plalist, int flag, char *errtext)
{
  static int jdbase =  INVALID_BASE;
  static int lastplalist = 0;
  static double lon[EP_NP][EPBS]; // buffer for 20 days unpacked ephe 
  static double out[2 * EP_NP];	 // buffer for return longitude and return speed 
  int p, pf;
  int ix, jdlong, iflagret;
  double lp;
  double jfract;
  double x[6];
  if (errtext != NULL)
    *errtext = '\0';
  if (plalist == 0)
    plalist = EP_ALL_BITS;		/* default: all logitudes, no speeds */
  /*
   * we must determine when to reload the lon buffer: if the contents do
   * not allow immediate interpolation or if the plalist selection has
   * changed since the last call.
   */
  if ((plalist & lastplalist) != plalist) {	/* new set is not contained in old */
    jdbase = INVALID_BASE;
  }
  lastplalist = plalist;
  jdlong = floor(jd - 0.5);
  ix = jdlong - jdbase;
  if (ix < EP_MIN_IX || ix >= EPBS) {	/* must reload full buffer */
    jdbase = ((jdlong - EP_MIN_IX) / NDB) * NDB;		/* new base */
    if (jdbase > jdlong - EP_MIN_IX) jdbase -= NDB;	/* fix bug for neg. */
    if (ephe4_unpack_d(jdbase, plalist, lon, 0, errtext) != OK)
      goto err_exit;
    if (ephe4_unpack_d(jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
      goto err_exit;
    ix = jdlong - jdbase;
  } else if (ix > EP_MAX_IX) {		/* must shift upper half down
					   and reload upper half of buffer */
    jdbase +=  NDB;		/* new base */
    for (p = 0; p < EP_NP; p++)
      memcpy(&lon[p][0], &lon[p][NDB], NDB * sizeof(double));
    if (ephe4_unpack_d(jdbase + NDB, plalist, lon, 0 + NDB, errtext) != OK)
      goto err_exit;
    ix = jdlong - jdbase;
  }
  jfract = jd - 0.5 - jdlong;        
  /*
   * we use the interpolator even for jfract = 0, because it delivers
   * the speed term. The computation overhead is unimportant
   * in any case.
   */
  for (p = 0, pf = 1; p < EP_NP; p++, pf = pf << 1)
    if ((plalist & pf) != 0) {
      inpolq((int) ix, qod[p], jfract, &(lon[p][0]), &(out[p]), &lp);
      if (p <= PLACALC_CHIRON) {	/* normalize all except ecl and nut */
	if (out[p] < 0)
	  out[p] += 360.0;
	else if (out[p] >= 360.0)
	  out[p] -= 360.0;
      }
#ifdef DEBUG
      fprintf(stderr,"ephread p=%d, lon=%.3lf\n", p, out[p]);
#endif
      if (flag & EP_BIT_SPEED)
	out[p+EP_NP] = lp;
    }
  return out;
err_exit:
  jdbase = INVALID_BASE;
  lastplalist = 0;
  if ((flag & EP_BIT_MUST_USE_EPHE) == 0) {	/* try using calc */
    int sweflag = 0;
    char serr[AS_MAXCH];
    if (flag & EP_BIT_SPEED)
      sweflag = SEFLG_SPEED;
    if (errtext != NULL)
      sprintf(errtext,"ephread failed for jd=%f; used swe_calc().", jd);
    for (p = 0, pf = 1; p < PLACALC_CALC_N; p++, pf = pf << 1) {
      if ((plalist & pf) != 0) {
	if ((iflagret = swe_calc(jd, ephe_plac2swe(p), sweflag, x, serr)) != ERR) {
	  out[p] = x[0];
	  if (flag & EP_BIT_SPEED)
	    out[p + EP_NP] = x[3];
	  if (out[p] < 0)
	    out[p] += 360.0;
	  else if (out[p] >= 360.0)
	    out[p] -= 360.0;
	} else {
	  swe_close();
          if (errtext != NULL)
	    strcat(errtext," calc failed too.");
          return NULL;
	}
      }
    }
    if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
      swe_close();
      sprintf(errtext, "error in swe_calc() %s\n", serr);
      return NULL;
    }
    out[EP_ECL_INDEX] = x[0];	/* true ecliptic */
    out[EP_NUT_INDEX] = x[2];	/* nutation */
    out[EP_ECL_INDEX + EP_NP] = 0;
    out[EP_NUT_INDEX + EP_NP] = 0;
    return out;
  }
  return NULL;
}	

/****************************************************
  unpack an ephe file block specified by jlong 
  and the planets specified by pflag into
  the array lon[p][EPBS], starting at index i0.
  jdl is (long) floor(full julian date);
  If the reading failes, ERR is returned and
  an error text of max. 79 char in errs, except when errs = NULL.
****************************************************/
static int ephe4_unpack(int jdl, int plalist, centisec lon[][EPBS], int i0,char *errs)
{
  int p, i, pf;
  centisec l_ret, d_ret;
  struct ep4 e;
  if (eph4_posit (jdl, FALSE, errs) != OK)
    return (ERR);
  if (fread (&e, sizeof(struct ep4), 1, ephfp) != 1) {
    if (errs != NULL)
      sprintf (errs, "ephe4_unpack: fread for jd=%d failed", jdl);
    return (ERR);
  }
#ifdef INTEL_BYTE_ORDER
  shortreorder((UCHAR *) &e, sizeof(struct ep4));
#endif
  for (p = PLACALC_SUN, pf = 1; p <= PLACALC_CHIRON; p++, pf = pf << 1) {
    if ((plalist & pf) == 0) continue; 
    l_ret = e.elo[p].p0m * 6000L + e.elo[p].p0s;	/* csec */
    d_ret = e.elo[p].pd1m * 6000L + e.elo[p].pd1s;	/* csec */
    lon[p][i0] = l_ret;
    l_ret +=  d_ret;
    if (l_ret < 0) {
      lon[p][i0+1] = l_ret + DEG360;
    } else if (l_ret >= DEG360) {
      lon[p][i0+1] = l_ret - DEG360;
    } else {
      lon[p][i0+1] = l_ret;
    }
    for (i = 2; i < NDB; i++) {
      if (p == PLACALC_MOON || p == PLACALC_MERCURY)
	d_ret += e.elo[p].pd2[i-2] * 10L;
      else
	d_ret += e.elo[p].pd2[i-2];
      l_ret += d_ret;
      if (l_ret < 0) {
	lon[p][i0+i] = l_ret + DEG360;
      } else if (l_ret >= DEG360) {
	lon[p][i0+i] = l_ret - DEG360;
      } else {
	lon[p][i0+i] = l_ret;
      }
    }
  }	/* for p */
  if ( plalist & EP_ECL_BIT) { /* unpack ecl */
    l_ret = e.ecl0m * 6000L + e.ecl0s;
    lon[EP_ECL_INDEX][i0] = l_ret;
    for (i = 1; i < NDB; i++)
      lon[EP_ECL_INDEX][i0+i] =  l_ret + e.ecld1[i-1]; 
  }
  if ( plalist & EP_NUT_BIT) { /* unpack nut */
    for (i = 0; i < NDB; i++) 
      lon[EP_NUT_INDEX][i0+i] = e.nuts[i] ;
  }
  return OK;
}	/* ephe4_unpack */

// same in double
static int ephe4_unpack_d(int jdl, int plalist, double lon[][EPBS], int i0,char *errs)
{
  int p, i, pf;
  double l_ret, d_ret;
  struct ep4 e;
  if (eph4_posit (jdl, FALSE, errs) != OK)
    return (ERR);
  if (fread (&e, sizeof(struct ep4), 1, ephfp) != 1) {
    if (errs != NULL)
      sprintf (errs, "ephe4_unpack: fread for jd=%d failed", jdl);
    return (ERR);
  }
#ifdef INTEL_BYTE_ORDER
  shortreorder((UCHAR *) &e, sizeof(struct ep4));
#endif
  for (p = PLACALC_SUN, pf = 1; p <= PLACALC_CHIRON; p++, pf = pf << 1) {
    if ((plalist & pf) == 0) continue; 
    l_ret = (e.elo[p].p0m * 6000 + e.elo[p].p0s) * CS2DEG;	
    d_ret = (e.elo[p].pd1m * 6000 + e.elo[p].pd1s) * CS2DEG;	
    lon[p][i0] = l_ret;
    l_ret += d_ret;
    if (l_ret < 0) {
      lon[p][i0+1] = l_ret + 360.0;
    } else if (l_ret >= 360.0) {
      lon[p][i0+1] = l_ret - 360.0;
    } else {
      lon[p][i0+1] = l_ret;
    }
    for (i = 2; i < NDB; i++) {
      if (p == PLACALC_MOON || p == PLACALC_MERCURY)
	d_ret += (e.elo[p].pd2[i-2] * 10 * CS2DEG);
      else
	d_ret += (e.elo[p].pd2[i-2] * CS2DEG);
      l_ret += d_ret;
      if (l_ret < 0) {
	lon[p][i0+i] = l_ret + 360.0;
      } else if (l_ret >= 360.0) {
	lon[p][i0+i] = l_ret - 360.0;
      } else {
	lon[p][i0+i] = l_ret;
      }
    }
  }	/* for p */
  if ( plalist & EP_ECL_BIT) { /* unpack ecl */
    l_ret = (e.ecl0m * 6000L + e.ecl0s) * CS2DEG;
    lon[EP_ECL_INDEX][i0] = l_ret;
    for (i = 1; i < NDB; i++)
      lon[EP_ECL_INDEX][i0+i] =  l_ret + e.ecld1[i-1] * CS2DEG; 
  }
  if ( plalist & EP_NUT_BIT) { /* unpack nut */
    for (i = 0; i < NDB; i++) 
      lon[EP_NUT_INDEX][i0+i] = e.nuts[i] * CS2DEG ;
  }
  return OK;
}	

/****************************************************
  position ephe file at proper position for julian 
  date jd; if writeflag = TRUE (write mode), create file
  if required. Return OK or ERR.
  globals used: ephfp.
*****************************************************/
int eph4_posit (int jlong, AS_BOOL writeflag, char *errtext)
{
  int filenr;
  long posit;
  static int open_filenr = -10000;
  char fname[AS_MAXCH], s[80], *sp;
  filenr = jlong / EP4_NDAYS;   
  if (jlong < 0 && filenr * EP4_NDAYS != jlong) filenr--;
  posit = jlong - filenr * EP4_NDAYS;
  posit = posit / NDB * sizeof(struct ep4);
  if (open_filenr  != filenr) {
    if (ephfp != NULL) {
      fclose(ephfp);
      open_filenr = -10000;
    }
    if (filenr >= 0)
      sprintf (s, "%s%s%d", EP4_PATH, EP4_FILE, filenr);
    else
      sprintf (s, "%s%sM%d", EP4_PATH, EP4_FILE, -filenr);
    my_makepath(fname, s);
    if (writeflag)
      sp = BFILE_W_CREATE;
    else
      sp = BFILE_R_ACCESS;
    ephfp = fopen (fname, sp);
    if (ephfp == NULL) {
      if (errtext != NULL) {
	if (! writeflag)  {
	  sprintf (errtext,"eph4_posit: file %s does not exist\n",  fname);
	} else {
	  sprintf (errtext,"eph4_posit: could not create file %s\n",  fname);
	}
      }
      return (ERR);	
    }
    open_filenr = filenr;
  }
  if  (fseek (ephfp, posit, 0) == 0 && ftell(ephfp) == posit) {
    return (OK);
  } else {
    if (errtext != NULL)
      sprintf (errtext,"eph4_posit: fseek(%ld) of file nr %d failed\n",  
	      posit, open_filenr);
    return (ERR);
  }
}	/* end eph4_posit */

/*****************************************************
quicker Everett interpolation, after Pottenger
version  for long, 17.7.91 by Alois Treindl
*****************************************************/
static void inpolq_l(int n, int o, double p, centisec *x, centisec *axu, centisec *adxu)
/* 
 * interpolate between x[n] and x[n-1], at argument n+p 
 * o = order of interpolation, maximum 5 
 * p = argument in [0..1] 
 * x[] array of function values, x[n-2]..x[n+3] must exist 
 * axu pointer for storage of result 
 * adxu pointer for storage of dx/dt  
 */
{
  static double	q,q2,q3,q4,q5,
	  p2,p3,p4,p5,
	  u,u0,u1,u2;
  static double lastp = 9999;
  double rl, rlp;
  centisec	dm2,dm1,d0,dp1,dp2,
	  d2m1,d20,d2p1,d2p2,
	  d30,d3p1,d3p2,
	  d4p1,d4p2;
  centisec offset = 0;
  if (lastp != p) {	/* recompute the interpolator factors */
    q=1.0-p;
    q2 = q*q;
    q3 = (q+1.0)*q*(q-1.0)/6.0;	/* q - 1  over 3; u5 */
    p2 = p*p;
    p3 = (p+1.0)*p*(p-1.0)/6.0;	/* p - 1  over 3; u8 */
    u = (3.0*p2-1.0)/6.0;
    u0 = (3.0*q2-1.0)/6.0;
    q4 = q2*q2;		/* f5 */
    p4 = p2*p2;		/* f4 */
    u1 = (5.0*p4-15.0*p2+4.0)/120.0;	/* u1 */
    u2 = (5.0*q4-15.0*q2+4.0)/120.0;	/* u2 */
    q5 = q3*(q+2.0)*(q-2.0)/20.0;   /* q - 2  over 5; u6 */
    p5 = (p+2.0)*p3*(p-2.0)/20.0;	/* p - 2  over 5; u9 */
    lastp = p;
  }
  dm1 = x[n]   - x[n-1];
  if (dm1 >= DEG180)
    dm1 -= DEG360;
  else if (dm1 < -DEG180)
    dm1 += DEG360;
  d0  = x[n+1] - x[n];
  if (d0 >= DEG180) {
    d0 -= DEG360;
    offset = DEG360;
  } else if (d0 < -DEG180) {
    d0 += DEG360;
    offset = -DEG360;
  }
  dp1 = x[n+2] - x[n+1];
  if (dp1 >= DEG180)
    dp1 -= DEG360;
  else if (dp1 < -DEG180)
    dp1 += DEG360;
  d20  = d0 - dm1;	/* f8 */
  d2p1 = dp1 - d0;	/* f9 */
  /*
   * Everett interpolation 3rd order
   */
  rl = q*(x[n] + offset)   + q3*d20 + p*x[n+1]  + p3*d2p1;
  rlp = d0 + u*d2p1  - u0*d20;
  if ( o > 3 ) {	/* 5th order */
    dm2 = x[n-1] - x[n-2];
    if (dm2 >= DEG180)
      dm2 -= DEG360;
    else if (dm2 < -DEG180)
      dm2 += DEG360;
    dp2 = x[n+3] - x[n+2];
    if (dp2 >= DEG180)
      dp2 -= DEG360;
    else if (dp2 < -DEG180)
      dp2 += DEG360;
    d2m1 = dm1 - dm2;
    d2p2 = dp2 - dp1;
    d30  = d20 - d2m1;
    d3p1 = d2p1 - d20;
    d3p2 = d2p2 - d2p1;
    d4p1 = d3p1 - d30;	/* f7 */
    d4p2 = d3p2 - d3p1;	/* f */
    rl  += p5*d4p2 + q5*d4p1;
    rlp += u1*d4p2 - u2*d4p1;
  }
  *axu = swe_d2l (rl);
  *adxu = swe_d2l (rlp);
}	/* end inpolq_l() */

/*****************************************************
quicker Everett interpolation, after Pottenger
version for double  9 Jul 1988    by Alois Treindl
return OK, no error checking
Was used in Placalc to interpolate 80-day stored ephe for outer planets.
*****************************************************/
static int inpolq(int n, int o, double p, double *x, double *axu, double *adxu)
// n	interpolate between x[n] and x[n-1], at argument n+p 
// o	order of interpolation, maximum 5 
// p,	argument , intervall [0..1] 
// x[]	array of function values, x[n-o]..x[n+o] must exist 
// *axu	pointer for storage of result 
// *adxu pointer for storage of dx/dt  
{
  static double	q,q2,q3,q4,q5,p2,p3,p4,p5, u,u0,u1,u2;
  static double lastp = 9999.0;
  double	dm2,dm1,d0,dp1,dp2,
	  d2m1,d20,d2p1,d2p2,
	  d30,d3p1,d3p2,
	  d4p1,d4p2;
  double offset = 0.0;
  if (lastp != p) {
    q=1.0-p;
    q2 = q*q;
    q3 = (q+1.0)*q*(q-1.0)/6.0;	/* q - 1  over 3; u5 */
    p2 = p*p;
    p3 = (p+1.0)*p*(p-1.0)/6.0;	/* p - 1  over 3; u8 */
    u = (3.0*p2-1.0)/6.0;
    u0 = (3.0*q2-1.0)/6.0;
    q4 = q2*q2;		/* f5 */
    p4 = p2*p2;		/* f4 */
    u1 = (5.0*p4-15.0*p2+4.0)/120.0;	/* u1 */
    u2 = (5.0*q4-15.0*q2+4.0)/120.0;	/* u2 */
    q5 = q3*(q+2.0)*(q-2.0)/20.0;   /* q - 2  over 5; u6 */
    p5 = (p+2.0)*p3*(p-2.0)/20.0;	/* p - 2  over 5; u9 */
    lastp = p;
  }
  dm1 = x[n]   - x[n-1];
  if (dm1 > 180.0) dm1 -= 360.0;
  if (dm1 < -180.0) dm1 += 360.0;
  d0  = x[n+1] - x[n];
  if (d0 > 180.0) {
    d0 -= 360.0;
    offset = 360.0;
  }
  if (d0 < -180.0) {
    d0 += 360.0;
    offset = -360.0;
  }
  dp1 = x[n+2] - x[n+1];
  if (dp1 > 180.0) dp1 -= 360.0;
  if (dp1 < -180.0) dp1 += 360.0;
  d20  = d0 - dm1;	/* f8 */
  d2p1 = dp1 - d0;	/* f9 */
  /* Everett interpolation 3rd order */
  *axu = q*(x[n] + offset)   + q3*d20
       + p*x[n+1]  + p3*d2p1;
  *adxu  =  d0 + u*d2p1  - u0*d20;
  if ( o > 3 ) {	/* 5th order */
    dm2 = x[n-1] - x[n-2];
    if (dm2 > 180.0) dm2 -= 360.0;
    if (dm2 < -180.0) dm2 += 360.0;
    dp2 = x[n+3] - x[n+2];
    if (dp2 > 180.0) dp2 -= 360.0;
    if (dp2 < -180.0) dp2 += 360.0;
    d2m1 = dm1 - dm2;
    d2p2 = dp2 - dp1;
    d30  = d20 - d2m1;
    d3p1 = d2p1 - d20;
    d3p2 = d2p2 - d2p1;
    d4p1 = d3p1 - d30;	/* f7 */
    d4p2 = d3p2 - d3p1;	/* f */
    *axu  += p5*d4p2 + q5*d4p1;
    *adxu += u1*d4p2 - u2*d4p1;
  }
  return (OK);
}	/* end inpolq() */


static char *my_makepath(char *d, char *s)
{
  char *getenv();
  if (*s == *DIR_GLUE || *s == '/' || strchr (s, ':') != NULL) {
    strcpy (d, s);	/* s is absolute path name */
  }
# if MSDOS
  while ((p = strchr(d, '/')) != NULL) *p = '\\';
# endif
  return (d);
}

int ephe_plac2swe(int p)
{
  if (p >= PLACALC_SUN && p <= PLACALC_TRUE_NODE) return p;
  if (p == PLACALC_CHIRON) return SE_CHIRON;
  if (p == PLACALC_LILITH) return SE_MEAN_APOG;
  if (p == PLACALC_CERES) return SE_CERES;
  if (p == PLACALC_PALLAS) return SE_PALLAS;
  if (p == PLACALC_JUNO) return SE_JUNO;
  if (p == PLACALC_VESTA) return SE_VESTA;
  if (p == PLACALC_EARTHHEL) return SE_EARTH;
  return -1;
}

int old_d2l (double x)		
{
  if (x >=0)
    return ((int) (x + 0.5));
  else
    return (- (int) (0.5 - x));
}