placalc.c

占星术4.0源码

   * otherwise as universal time.
   */
  if ((flag & CALC_BIT_EPHE) == 0) {
    jd_ad += deltat (jd_ad);
  }
  for (p = SUN; p < CALC_N; p++) {
    if (! check_bit(plalist, p)) continue;
    if (calc (p, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
      lcs [p] = d2l (rlng * DEG);
      lpcs [p] = d2l (rspeed * DEG);
      betcs [p] = d2l (rlat * DEG);
      rau [p] = rrad;
    } else {
      sprintf(so,"error at planet %d", p);
      return ( ERR);
    }
  }
  /*
   * format comma separated list: id,teph,flag,plalist,ekl,nut
   * REAL8 is given with 8 digits precision after decimal point, 
   * all angles are given in centiseconds.
   * then for each requested planet: longitude (csec)
   * then for each requested planet, if wanted: speed (csec/day)
   * then for each requested planet, if wanted: latitude (csec)
   * then for each requested planet, if wanted: rgeo (units 0..999)
   * then for each requested planet, if wanted: rau  (A.U.)
   */
  sprintf (so, "%d,%.8lf,%d,%d,%ld,%ld", id, jd_ad, flag, plalist, 
	       d2l(ekl * DEG), d2l (nut * DEG) );
  so_len = strlen (so);
  for (planet = SUN; planet < CALC_N; planet++) {
    if (! check_bit(plalist, planet)) continue;
    sprintf (s ,",%ld", lcs[planet]);
    strcat (so + so_len, s);
    so_len += strlen (s);
  }
  if (flag & CALC_BIT_SPEED) {
    for (planet = SUN; planet < CALC_N; planet++)  {
      if (! check_bit(plalist, planet)) continue;
      sprintf (s ,",%ld", lpcs[planet]);
      strcat (so + so_len, s);
      so_len += strlen (s);
    }
  }
  if (flag & CALC_BIT_BETA) {
    for (planet = SUN; planet < CALC_N; planet++)  {
      if (! check_bit(plalist, planet)) continue;
      sprintf (s ,",%ld", betcs[planet]);
      strcat (so + so_len, s);
      so_len += strlen (s);
    }
  }
  if (flag & CALC_BIT_RGEO) {
    for (planet = SUN; planet < CALC_N; planet++)  {
      if (! check_bit(plalist, planet)) continue;
      sprintf (s ,",%ld", rel_geo(planet,rau[planet]));
      strcat (so + so_len, s);
      so_len += strlen (s);
    }
  }
  if (flag & CALC_BIT_RAU) {
    for (planet = SUN; planet < CALC_N; planet++)  {
      if (! check_bit(plalist, planet)) continue;
      sprintf (s ,",%.8lf", rau[planet]);
      strcat (so + so_len, s);
      so_len += strlen (s);
    }
  }
  return (OK);
}	/* end calcserv */
#endif /* !ASTROLOG */

/******************************************************************
   function calc(): 
   This is the main routine for computing a planets position.
   The function has several modes, which are controlled by bits in
   the parameter 'flag'. The normal mode (flag == 0) computes
   a planets apparent geocentric position in ecliptic coordinates relative to
   the true equinox of date, without speed

   Explanation of the arguments: see the functions header.

   Returns OK or ERR (if some planet out of time range). OK and ERR are
   defined in ourdef.h and must not be confused with TRUE and FALSE.
   OK and ERR are of type int, not of type BOOLEAN.

   Bits used in flag:
   CALC_BIT_HELIO  		0 = geocentric, 1 = heliocentric
   CALC_BIT_NOAPP 	   	0 = apparent positions, 1 = true positions
   CALC_BIT_NONUT 		0 = do nutation (true equinox of date)
				1 = don't do nutation (mean equinox of date).

   CALC_BIT_SPEED 		0 = don't calc speed,
				1 = calc speed, takes quite long for moon
				(is observed only for moon, with other
				planets speed is cheap)

   Side effects and local memory:
   For doing heliocentric positions the fucntion must know the
   earth's position for the desired time t. It remembers the earth
   position so it does not have to recompute it each time a planet
   position is wanted for the same time t.
   It calls helup(t), which leaves as a side effect the global
   variables meanekl, ekl and nut for the time t.

   Functions called by calc():
   helup(t)
   hel(t)
   moon(t)
   togeo()

   Time range:
   The function can be used savely in the time range 5000 BC to
   3000 AD. The stored ephemeris is available only for this time
   range, so Jupiter ... Pluto cannot be computed outside. The
   function will return results for the other planets also outside
   of this time range, but they become meaningless pretty soon
   before 5000 BC, because Newcombs time series expansions for the
   elements will not work anymore.

******************************************************************/
int calc(int  planet,  	/* planet index as defined in placalc.h,
			 SUN = 0, MOON = 1 etc.
			 planet == -1 calc calculates only nut and ecl */
	 REAL8 jd_ad,	/* relative Astrodienst Juldate, ephemeris time.
			 Astrodienst Juldate is relative 31 Dec 1949, noon. */
	 int  flag,	/* See definition of flag bits above */
	 REAL8 *alng,
	 REAL8 *arad,
	 REAL8 *alat,
	 REAL8 *alngspeed)
      /* pointers to the return variables:
	 alng = ecliptic longitude in degrees
	 arad = radius vector in AU (astronomic units)
	 alat = ecliptic latitude in degrees
	 alngspeed = speed of planet in degrees per day
	 
	 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	 The precision of the speed is quite limited.
	 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	 For Sun, Mercury, Venus and Mars we take only the speed from
	 the undisturbed Kepler orbit. For the Moon there is no
	 reasonable undisturbed orbit and we derive the speed from
	 its position at t + dt and t - dt. We need these 
	 moon positions anyway for the true node calculation.
	 For the outer planets and Chiron we derive the precise
	 speed from the stored ephemeris by high order inter-
	 polation; the precision is limited for the geocentric
	 case due to the limited precision of the earth's/sun's
	 speed.
	 Applications who need precise speeds should
	 get them by calling calc() with slightly different times.
	 */
  /*
   * Comment 7 May 1991 by Alois Treindl:
   * Center of Earth versus Barycenter Earth-Moon:
   * Brown's theory of the moon gives the moon's coordinates relative
   * to the center of the earth. Newcomb's theory of the Sun gives the
   * coordinates of the earth's center relative to the center of the Sun.
   * This is what we need.
   *
   * How about the Mean Lunar Node?
   * The orbital elements of the Sun in Newcomb's theory are given
   * relative to the barycenter Earth-Moon; the reduction to geocentric
   * is only applied after doing the Kepler ellipse calculation.
   * Are the Lunar elements also relative to the barycenter??
   * If yes:
   * When we use the moon's mean node out of the elements, it is still
   * as seen from the barycenter. Because the node is close to the
   * earth, we would have to apply a considerable correction, which is of
   * the order of 4000/384000 km or 35' (minutes of arc).
   * Nobody has ever applied such a correction to the mean node.
   *
   * And the True Node?
   * When we calculate the osculating orbital elements of the Moon (true node),
   * are they relative to the barycenter or to the Earth's center?
   * Our derivation of true node from the actual Moon positions considers
   * the earth's center as the focal point of the osculating lunar ellipse.
   * A more correct approach would first reduce the lunar position from
   * geocentric to barycentric, then compute the orbital elements from
   * the reduced positions, and then reduce the desired items
   * (node, apogaeum, 'dark moon') to geocentric positions.
   * No known astrological ephemeris has ever used such a correction, which is
   * of the same order of magnitude as the correction to the meannode above.
   * When the moon is going through the ecliptic, the geocenter, barycenter
   * moon (and the node identical to the moon itself) line up; this is why
   * the error does not show up in normal considerations.
   */
{
  struct rememberdat  /* time for which the datas are calculated */
    { REAL8  calculation_time, lng, rad, zet, lngspeed, radspeed, zetspeed; };
  static struct rememberdat earthrem = 
    { HUGE, HUGE, HUGE, HUGE, HUGE, HUGE, HUGE };
  static struct rememberdat moonrem  = 
    { HUGE, HUGE, HUGE, HUGE, HUGE, HUGE, HUGE };
  REAL8 c, s, x, knn, knv; 
  REAL8 rp, zp; /* needed to call hel! */
  REAL8 *azet = alat; 
  BOOLEAN calc_geo, calc_helio, calc_apparent, calc_speed,
	  calc_nut;

  helup (jd_ad); /* helup checks whether it was already called with same time*/
  if ( planet == CALC_ONLY_ECL_NUT ) return (OK);  

  calc_helio =  flag & CALC_BIT_HELIO;
  calc_geo = ! calc_helio;
  calc_apparent = ! (flag & CALC_BIT_NOAPP);
  calc_nut = ! (flag & CALC_BIT_NONUT);
  calc_speed = flag & CALC_BIT_SPEED;
  /*
   * it is necessary to compute EARTH in the following cases:
   * heliocentric MOON or EARTH
   * geocentric any planet except MOON or nodes or LILITH
   */
  if (calc_helio && (planet == MOON || planet == EARTH)
  || calc_geo && planet != MOON
	      && planet != MEAN_NODE
	      && planet != TRUE_NODE
              && planet != LILITH) {
    if ( earthrem.calculation_time != jd_ad ) {	 
      hel (  EARTH, jd_ad, alng, arad, azet, alngspeed, &rp, &zp );
      /* store earthdata for geocentric calculation: */
      earthrem.lng = *alng;  
      earthrem.rad = *arad;
      earthrem.zet = *azet;
      earthrem.lngspeed = *alngspeed;
      earthrem.radspeed = rp;
      earthrem.zetspeed = zp;
      earthrem.calculation_time = jd_ad;
    }
  }
  switch(planet) {
  case EARTH:	/* has been already computed */
    *alng = earthrem.lng;
    *arad = earthrem.rad;
    *azet = earthrem.zet;
    *alngspeed = earthrem.lngspeed;
    rp = earthrem.radspeed;
    zp = earthrem.zetspeed;
    if ( calc_geo ) {	/* SUN seen from earth */
      *alng = smod8360( *alng + 180.0 );	
      *azet = - *azet;
    }
    if (calc_apparent)
      *alng = *alng - 0.0057683 * (*arad) * (*alngspeed);
    break; 
  case MOON:
    moon( alng, arad, azet );  
    moonrem.lng = *alng;  /* moonrem will be used for TRUE_NODE */
    moonrem.rad = *arad;
    moonrem.zet = *azet;
    *alngspeed = 12;
    moonrem.calculation_time = jd_ad;
    if ( calc_helio || calc_speed ) {/* get a second moon position */
      REAL8 lng2, rad2, zet2;
      helup( jd_ad + MOON_SPEED_INTERVAL );		
      moon( &lng2, &rad2, &zet2 );  
      helup( jd_ad );
      if ( calc_helio ) { /* moon as seen from sun */
	togeo( earthrem.lng, -earthrem.rad, moonrem.lng, moonrem.rad,
	     moonrem.zet, alng, arad );
	togeo( earthrem.lng + MOON_SPEED_INTERVAL * earthrem.lngspeed, 
	    -( earthrem.rad + MOON_SPEED_INTERVAL * earthrem.radspeed ), 
	    lng2, rad2, zet2, &lng2, &rad2 );
      }
      *alngspeed =  diff8360( lng2, *alng ) / MOON_SPEED_INTERVAL;
      /* rp = ( rad2 - *arad ) / MOON_SPEED_INTERVAL;
	 zp = ( zet2 - moonrem.zet ) / MOON_SPEED_INTERVAL; */
    }
    *alat = RADTODEG * ASIN8( *azet / *arad );
    /*
     * light time correction, not applied for moon or nodes;
     * moon would have only term of ca. 0.02", see Expl.Sup.1961 p.109
     */
    break;
  case MERCURY:
  case VENUS:
  case MARS:
  case JUPITER:
  case SATURN:
  case URANUS:
  case NEPTUNE:
  case PLUTO:
  case CHIRON:
    if (hel ( planet, jd_ad, alng, arad, azet, alngspeed, &rp, &zp ) != OK)
      return (ERR);	/* outer planets can fail if out of ephemeris range */
    if ( calc_geo ) {       /* geocentric */
      REAL8 lng1, rad1, lng2, rad2;
      togeo( earthrem.lng, earthrem.rad, *alng, *arad, *azet, &lng1, &rad1 );
      togeo( earthrem.lng + earthrem.lngspeed, 
	   earthrem.rad + earthrem.radspeed, 
	   *alng + *alngspeed, *arad + rp, *azet + zp, &lng2, &rad2 ); 
      *alng = lng1;
      *arad = rad1; 
      *alngspeed = diff8360( lng2, lng1 );
      /* rp = rad2 - rad1; */
    }
    *alat = RADTODEG * ASIN8( *azet / *arad );
    if (calc_apparent)
      *alng = *alng - 0.0057683 * (*arad) * (*alngspeed);
    break;
  case MEAN_NODE:
    *alng = smod8360( el[MOON].kn);
    /*
     * the distance of the node is the 'orbital parameter' p = a (1-e^2);
     * Our current use of the axis a is wrong, but is never used.
     */
    *arad = pd[MOON].axis;
    *alat = 0.0;
    *alngspeed = -0.053;
    break;
  case TRUE_NODE: {
    /* see comment 'Note 7 May 1991' above */
    REAL8  ln, rn, zn, 
	   lv, rv, zv, 
	   l1, r1, z1, 
	   xn, yn, xv, yv, r0, x0, y0;
    helup( jd_ad + NODE_INTERVAL );
    moon( &ln, &rn, &zn );
    helup( jd_ad - NODE_INTERVAL );
    moon( &lv, &rv, &zv );
    helup( jd_ad );
    if ( moonrem.calculation_time != jd_ad ) 
      moon( &l1, &r1, &z1 );  
    else { 	/* moon is already calculated */
      l1 = moonrem.lng;
      r1 = moonrem.rad;
      z1 = moonrem.zet;
    }
    rn = sqrt( rn * rn - zn * zn );
    rv = sqrt( rv * rv - zv * zv );
    r0 = sqrt( r1 * r1 - z1 * z1 );
    xn = rn * COS8( DEGTORAD * ln );
    yn = rn * SIN8( DEGTORAD * ln );
    xv = rv * COS8( DEGTORAD * lv );
    yv = rv * SIN8( DEGTORAD * lv );
    x0 = r0 * COS8( DEGTORAD * l1 );   
    y0 = r0 * SIN8( DEGTORAD * l1 );   
    x = test_near_zero( x0 * yn - xn * y0 );
    s = ( y0 * zn - z1 * yn ) / x;
    c = test_near_zero( ( x0 * zn - z1 * xn ) / x );
    knn =  smod8360( RADTODEG * ATAN28( s, c )); /* = ATAN8( s / c ) */
    x = test_near_zero( y0 * xv - x0 * yv );
    s = ( yv * z1 - zv * y0 ) / x;
    c = test_near_zero( ( xv * z1 - zv * x0 ) / x );
    knv =  smod8360( RADTODEG * ATAN28( s, c ));        
    *alng = smod8360( ( knv + knn ) / 2 );
    /*
     * the distance of the node is the 'orbital parameter' p = a (1-e^2);