placalc.c

占星术4.0源码

	2362,	2386,	2449,	2434,	2408,
	2402,	2400,	2387,	2395,	2386,
	2393,	2373,	2392,	2396,	2402,
	2433,	2483,	2530,	2570,	2624,
	2677,	2728,	2778,	2825,	2871,
	2915,	2957,	2997,	3036,	3072,	/* 1950 - 1954 */
	3107,	3135,	3168,	3218,	3268,
	3315,	3359,	3400,	3447,	3503,
	3573,	3654,	3743,	3829,	3920,
	4018,	4117,	4223,	4337,	4449,
	4548,	4646,	4752,	4853,	4959,
	5054,	5138,	5217,	5296,	5379,	/* 1980 - 1984 */
	5434,	5487,	5532,	5582,	5630,	/* 1985 - 89 from AE 1991 */
	5686,	5757,	5900,	5900,	6000,	/* AE 1993 and extrapol */
	6050,	6100,	6150,	6200,	6250,	/* 1995 - 1999 */
	6300};					/* 2000 */
  double yr, cy, delta;
  long iyr, i;
  yr = (jd_ad + 18262) / 365.25 + 100.0;    /*  year  relative 1800 */
  cy = yr / 100;
  iyr =  (long) (floor (yr) + 1800);	 /* truncated to integer, rel 0 */
# if TIDAL_26			/* Stephenson formula only when 26" tidal
				   term in lunar motion */
  if ( iyr >= 1637  && iyr < 2000 ) {
    i = iyr - 1637;
    delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - floor (yr)) * 0.01;
  } else if (iyr >= 2000 ) {	/* parabola, fitted at value[2000] */
    delta = 25.5 * cy * cy  - 25.5 * 4 + 63.00;	
  } else if (iyr >= 948) {	/* from 948 - 1637 use parabola */
    delta = 25.5 * cy * cy;	
  } else {	/* before 984 use other parabola */	
    delta = 1361.7	+ 320 * cy + 44.3 * cy * cy;	/* fits at 948 */
  }
# else	/* use Clemence value + 5 sec before 1690, new dt afterwards */
  cy -= 1;	/* epoch 1900 */
  if ( iyr >= 1690  && iyr < 2000 ) {
    i = iyr - 1637;
    delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - floor (yr)) * 0.01;
  } else if (iyr >= 2000 ) {	/* parabola, fitted at value[2000] */
    delta = 29.949 * cy * cy  - 29.949 * 4 + 63.0;	
  } else {	
    delta = 5 + 24.349 + 72.3165 * cy + 29.949 * cy * cy; /* fits at 1690 */
  }
# endif
  return (delta / 86400.0);
}
#ifdef ASTROLOG
/* End contents of deltat.c */
#endif
#ifdef ASTROLOG
/* Begin contents of d2l.c */
#endif
/*******************************************
$Header: d2l.c,v 1.9 91/11/16 16:24:20 alois Exp $
********************************************/

/*************************************
double to long with rounding, no overflow check
*************************************/ 
long d2l (double x)		
{
  if (x >=0)
    return ((long) (x + 0.5));
  else
    return (- (long) (0.5 - x));
}
#ifdef ASTROLOG
/* End contents of d2l.c */
#endif
#ifdef ASTROLOG
/* Begin contents of csec.c */
#endif
/*
 * $Header$
 *
 * A collection of useful functions for centisec calculations.

  ---------------------------------------------------------------
  | Copyright Astrodienst Zurich AG and Alois Treindl, 1991.	|
  | The use of this source code is subject to regulations made	|
  | by Astrodienst Zurich. The code is NOT in the public domain.|
  |								|
  | This copyright notice must not be changed or removed	|
  | by any user of this program.				|
  ---------------------------------------------------------------
*******************************************************/
#ifndef ASTROLOG
#include "ourdef.h"
#include "astrolib.h"
#include "housasp.h"
#endif

/************************************
normalize argument into interval [0..DEG360]
*************************************/
centisec csnorm(centisec p)
{
  if (p < 0) 
    do { p += DEG360; } while (p < 0);
  else if (p >= DEG360)
    do { p -= DEG360; } while (p >= DEG360);
  return (p);
}

double degnorm(double p)
{
  if (p < 0) 
    do { p += 360.0; } while (p < 0);
  else if (p >= 360.0)
    do { p -= 360.0; } while (p >= 360.0);
  return (p);
}

/************************************
distance in centisecs p1 - p2
normalized to [0..360[
**************************************/
centisec difcsn (centisec p1, centisec p2)
{ 
  return (csnorm(p1 - p2));
}

double difdegn (double p1, double p2)
{ 
  return (degnorm(p1 - p2));
}

/************************************
distance in centisecs p1 - p2
normalized to [-180..180[
**************************************/
centisec difcs2n (centisec p1, centisec p2)
{ centisec dif;
  dif = csnorm(p1 - p2);
  if (dif  >= DEG180) return (dif - DEG360);
  return (dif);
}

double difdeg2n (double p1, double p2)
{ double dif;
  dif = degnorm(p1 - p2);
  if (dif  >= 180.0) return (dif - 360.0);
  return (dif);
}

/*************************************
round second, but at 29.5959 always down
*************************************/ 
centisec roundsec(centisec x)	
{
  centisec t;
  t = (x + 50) / 100 *100L;	/* round to seconds */
  if (t > x && t % DEG30 == 0)  /* was rounded up to next sign */
    t = x / 100 * 100L;		/* round last second of sign downwards */
  return (t);
}


/******************************/
double dcos(centisec x)
{
  return (COS8 (CSTORAD * x));
}

/******************************/
double dsin(centisec x)
{
  return (SIN8 (CSTORAD * x));
}

/******************************/
double dtan(centisec x)
{
  return (TAN8 (CSTORAD * x));
}

/******************************/
centisec datan(double x)
{
  return (d2l (RADTOCS * ATAN8 (x)) );	
}

/******************************/
centisec dasin(double x)
{
  return (d2l (RADTOCS * ASIN8 (x)));	
}
#ifdef ASTROLOG
/* End contents of csec.c */
#endif
#endif /* ASTROLOG */
#include <string.h>

#ifndef ASTROLOG
#ifndef DIR_GLUE
# if MSDOS
# define DIR_GLUE "\\"
# else
# define DIR_GLUE "/"
# endif
#endif
#else
#define DIR_GLUE ""
#endif /* ASTROLOG */

 
/************************************************************
externally accessible globals, defined as extern in placalc.h 
************************************************************/


REAL8 meanekl, ekl, nut;
struct elements el [MARS + 1];

/*
 * The global variable ephe_path indicates where the ephemeris files
 * LRZ5_xx and CHI_xx are stored.
 * By default it is set to the #defined EPHE_PATH, but the user of the
 * placalc module can change it to any other location before he
 * starts calling calc().
 */
char *ephe_path = EPHE_PATH;
/*
 * If there occurs an internal error in placalc, a message is 
 * written to stderr or into the string variable placalc_err_text.
 * By default it is written to stderr, but if placalc_err_text is
 * not a NULL pointer, it is written to the string variable.
 * The user must set this pointer to a string of at least 160 bytes long,
 * if he/she wants to use it.
 */
char *placalc_err_text = NULL;


#ifndef ASTROLOG
/**********************************************************
function nacalc	()
calculates an array of planet longitudes and speeds,
as needed for complete nathan data records.
The function knows itself how many planets and in which mode
they have to be calculated for Nathan. 

return OK or ERR

The returned positions are in centiseconds, our standard
coordinate format for fast mathematics with planetary positions.

This function is just a template of how the calc() package
can be used. 
**********************************************************/
int nacalc (REAL8	jd_ad,	/* universal time relative  julian date */
	    centisec *plon,	/* returned longitudes */
	    centisec *pspe	/* returned speeds, if not NULL pointer */
	   )
{
  int planet, flag;
  REAL8 rlng, rrad, rlat, rspeed;
  int result = OK;
  flag = CALC_BIT_SPEED;	/* same, with speed */
  jd_ad += deltat( jd_ad );	/* ET = UT + Delta_T */
  for (planet = SUN; planet <= TRUE_NODE; planet++) {
    if (calc (planet, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
      plon [planet] = d2l (rlng * DEG);
      if (pspe != NULL) pspe [planet] = d2l (rspeed * DEG);
    } else {
      plon [planet] = -1;
      if (pspe != NULL) pspe [planet] = 0;
      result = ERR;
    }
  }
  planet = TRUE_NODE + 1;	/* CHIRON may not be TRUE_NODE + 1 */
  if (calc (CHIRON, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
    plon [planet] = d2l (rlng * DEG);
    if (pspe != NULL) pspe [planet] = d2l (rspeed * DEG);
  } else {
    plon [planet] = -1;
    if (pspe != NULL) pspe [planet] = 0;
    result = ERR;
  }
  return result;
}	/* end nacalc */
#endif /* !ASTROLOG */

#ifdef ASTROLOG
/* Given an object index and a Julian Day time, get its zodiac and        */
/* declination position (planetary longitude and latitude) of the object  */
/* and its velocity and distance from the Earth or Sun. This basically    */
/* just calls the Placalc calculation function to actually do it, but as  */
/* this is the one routine called from Astrolog, this is the one routine  */
/* which has knowledge of and uses both Astrolog and Placalc definitions, */
/* and does things such as translation to Placalc indices and formats.    */

int PlacalcPlanet(ind, jd, helio, planet, planetalt, ret, space)
int ind, helio;
double jd;
real *planet, *planetalt, *ret, *space;
{
  int iplanet, flag;
  REAL8 jd_ad, rlng, rrad, rlat, rspeed;

  if (ind <= _PLU)    /* Convert Astrolog object index to Placalc index. */
    iplanet = ind-1;
  else if (ind == _CHI)
    iplanet = CHIRON;
  else if (ind == _NOD)
#ifdef TRUENODE
    iplanet = TRUE_NODE;
#else
    iplanet = MEAN_NODE;
#endif
  else
    return FALSE;

  jd_ad = jd - JUL_OFFSET;
  flag = helio ? CALC_BIT_SPEED | CALC_BIT_HELIO : CALC_BIT_SPEED;
  jd_ad += deltat(jd_ad);
  if (calc(iplanet, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
    *planet    = rlng;
    *planetalt = rlat;
    *ret       = rspeed;
    *space     = rrad;
    return TRUE;
  }
  return FALSE;
}
#endif /* ASTROLOG */

#ifndef ASTROLOG
/******************************************************************
 * calculation server
 * delivers positions in string format which can be sent easily
 * over a communication line to the calculation client.
 ******************************************************************/
int calcserv(int id,	/* request id, random number to prevent phase err */
	    REAL8 jd_ad,	/* time as relative Astrodienst julian date */
	    int	flag,	/* a set of CALC_BIT_ bitflags */
	    int  plalist,/* bit list of planets to be computed, 0 = all */
	    char *so)	/* output string, MUST BE LONG ENOUGH (800 bytes)*/
{
  int p, planet, so_len;
  REAL8 rlng, rrad, rlat, rspeed, rau[CALC_N];
  centisec lcs[CALC_N], lpcs[CALC_N], betcs[CALC_N]; 
#ifndef ASTROLOG
  int rgeo[CALC_N];
#endif
  char s[MAXCHAR];
  if (plalist == 0) plalist = CALC_ALL_PLANET_BITS;
  /*
   * flag determines whether deltat is added to t;
   * if CALC_BIT_EPHE is set, jd_ad is considered as ephemeris time,