placalc.c

占星术4.0源码

     * Our current use of the axis a is wrong.
     */
    *arad = pd[MOON].axis;	
    *alat = 0.0;
    *alngspeed = diff8360( knn, knv ) / NODE_INTERVAL;
    }
    break;
  case LILITH: {
    /*
     * Added 22-Jun-93
     * Lilith or Dark Moon is the empty focal point of the mean lunar ellipse.
     * This is 180 degrees from the perihel.
     * Because the lunar orbit is not in the ecliptic, it must be
     * projected onto the ecliptic in the same way as the planetary orbits
     * are (see for example Montenbruck, Grundlagen der Ephemeridenrechnung).
     *
     * We compute the MEAN Lilith, not the TRUE one which would have to be
     * derived in a similar way as the true node.
     * For the radius vector of Lilith we use a simple formula; 
     * to get a precise value, the fact that the focal point of the ellipse
     * is not at the center of the earth but at the barycenter moon-earth
     * would have to be accounted for.
     * For the speed we always return a constant: the T term from the
     * lunar perihel.
     * Joelle de Gravelaine publishes in her book "Lilith der schwarze Mond"
     * (Astrodata, 1990) an ephemeris which gives noon (12.00) positions
     * but does not project onto the ecliptic.
     * This creates deviations 
     */
    double arg_lat, lon, cosi;
    struct elements *e = &el[MOON];
    arg_lat = degnorm(e->pe - e->kn + 180.0);
    cosi = COSDEG(e->in);
    if (e->in == 0 || ABS8( arg_lat -  90.0 ) < TANERRLIMIT 
    || ABS8( arg_lat - 270.0 ) < TANERRLIMIT ) {
      lon = arg_lat;
    } else {
      lon = ATAN8( TANDEG( arg_lat ) * cosi );
      lon = RADTODEG * lon;
      if ( arg_lat > 90.0 && arg_lat < 270.0 )  lon += 180.0;
    }
    lon = degnorm(lon + e->kn);
    *alng = lon;
    *alngspeed = 0.111404;	/* 6'41.05" per day */
    *arad = 2 * pd[MOON].axis * e->ex;
    /*
     *  To test Gravalaines error, return unprojected long in alat.
     *  the correct latitude would be:
     *   *alat = RADTODEG * ASIN8(SINDEG(arg_lat) * SINDEG(e->in));
     */
    *alat = degnorm(arg_lat + e->kn);	/* unprojected longitude, no nut */
    }
    break;
  default:
    fprintf(stderr, "calc() called with illegal planet %d\n", planet);
    return ERR;
  } /* end switch */ 
  if (calc_nut)
    *alng += nut;
  *alng = smod8360( *alng);	/* normalize to circle */
  return (OK);
} /* end calc */

int rel_geo(int planet, double rau)
{
  /*
   * get relative earth distance in range 0..1000:
   * To compute the relative distance we use fixed values of 
   * mimimum and maximum distance measured empirically between
   * 1300 AD and 2300 AD (see helconst.c). 
   * This approach is certainly fine for the
   * outer planets, but not the best for Sun and Moon, where it
   * would be better to look at the mean anomaly, i.e. the progress
   * the planet makes on it's Kepler orbit.
   * Considering the low importance astrologers give to the relative
   * distance, we consider the effort not worthwile.
   * Now we compare real radius with longtime-averaged distances.
   */
  int rgeo;
  if (planet == MEAN_NODE || planet == TRUE_NODE || planet == LILITH) {
    return 0;
  } else {
#ifndef ASTROLOG
    rgeo = 1000 * (1.0 - (rau - rmima[planet][0]) / (rmima[planet][1] - rmima[planet][0]));
#else
    rgeo = 1000 * (int)(1.0 - (rau - rmima[planet][0]) / (rmima[planet][1] - rmima[planet][0]));
#endif
  }
  if (rgeo < 0)
    rgeo = 0;
  else if (rgeo > 999)
    rgeo = 999;
  return rgeo;
}

/******************************************************************
helio to geocentric conversion
******************************************************************/
void togeo(REAL8 lngearth,
	 REAL8 radearth,
	 REAL8 lng,
	 REAL8 rad,
	 REAL8 zet,
	 REAL8 *alnggeo,
	 REAL8 *aradgeo )
{
  REAL8 r1, x, y;
  r1 = sqrt( rad * rad - zet * zet );
  x = r1 * COS8( DEGTORAD * lng ) - radearth * COS8( DEGTORAD * lngearth );
  y = r1 * SIN8( DEGTORAD * lng ) - radearth * SIN8( DEGTORAD * lngearth );
  *aradgeo = sqrt( x * x + y * y + zet * zet );
  x = test_near_zero( x );
  *alnggeo =  smod8360( RADTODEG * ATAN28( y, x ) );
} /* end togeo */

/******************************************************************
helup()
prepares the orbital elements and the disturbation arguments for the
inner planets and the moon. helup(t) is called by hel() and by calc().
helup() returns its results in global variables.
helup() remembers the t it has been called with before and does
not recalculate its results when it is called more than once with
the same t.
******************************************************************/
void helup (REAL8  jd_ad ) /* relative julian date, ephemeris time */
{
  int i;
  static REAL8 thelup = HUGE;	/* is initialized only once at load time */
  struct elements *e = el;      /* pointer to el[i] */  
  struct elements *ee = el;     /* pointer to el[EARTH] */  
  struct eledata  *d = pd;      /* pointer to pd[i] */
  REAL8 td, ti, ti2, tj1, tj2, tj3;

  if ( thelup == jd_ad ) return;	/* if already calculated then return */
  for ( i = SUN; i <= MARS; i++, d++, e++ ) 
  {
    td = jd_ad - d->epoch;
    ti = e->tj = td / 36525.0;  /* julian centuries from epoch */
    ti2 = ti * ti;
    tj1 = ti / 3600.0;	/* used when coefficients are in seconds of arc */
    tj2 = ti * tj1;
    tj3 = ti * tj2;
    e->lg = mod8360( d->lg0 + d->lg1 * td  + d->lg2 * tj2 + d->lg3 * tj3 );
      /* also with moon lg1 *td is exact to 10e-8 degrees within 5000 years */
    e->pe = mod8360( d->pe0 + d->pe1 * tj1 + d->pe2 * tj2 + d->pe3 * tj3 );
    e->ex = d->ex0 + d->ex1 * ti + d->ex2 * ti2;
    e->kn = mod8360( d->kn0 + d->kn1 * tj1 + d->kn2 * tj2 + d->kn3 * tj3 );
    e->in = d->in0 + d->in1 * tj1 + d->in2 * tj2;
    e->ma = smod8360( e->lg - e->pe );
    if ( i == MOON ) {  /* calculate ekliptic according Newcomb, APAE VI,
			and nutation according Exp.Suppl. 1961, identical
			with Mark Potttenger elemnut()
			all terms >= 0.01" only .
			The 1984 IAU Theory of Nutation, as published in
			AE 1984 suppl. has not yet been implemented
			because it would mean to use other elements of
			moon and sun */
      REAL8 mnode, mlong2, slong2, mg, sg, d2;
      mnode  = DEGTORAD * e->kn;	/* moon's mean node */
      mlong2 = DEGTORAD * 2.0 * e->lg; 	/* 2 x moon's mean longitude */
      mg     = DEGTORAD * e->ma;	/* moon's mean anomaly (g1) */
      slong2 = DEGTORAD * 2.0 * ee->lg; /* 2 x sun's mean longitude (L), with
					the phase 180 deg earth-sun irrelevant
					because 2 x 180 = 360 deg */
      sg     = DEGTORAD * ee->ma;	/* sun's mean anomaly = earth's */
      d2     = mlong2 - slong2;		/* 2 x elongation of moon from sun */
      meanekl = ekld[0] + ekld[1] * tj1 + ekld[2] * tj2 + ekld[3] * tj3;
      ekl = meanekl + 
	    ( 9.2100 * COS8( mnode ) 
	    - 0.0904 * COS8( 2.0 * mnode ) 
	    + 0.0183 * COS8( mlong2 - mnode )
	    + 0.0884 * COS8( mlong2 ) 
	    + 0.0113 * COS8( mlong2 + mg )
	    + 0.5522 * COS8( slong2 )
	    + 0.0216 * COS8( slong2 + sg ) ) / 3600.0;
      nut = ( ( -17.2327 - 0.01737 * ti ) * SIN8( mnode ) 
     	    + 0.2088 * SIN8( 2.0 * mnode )
	    + 0.0675 * SIN8( mg )
	    - 0.0149 * SIN8( mg - d2 )
	    - 0.0342 * SIN8( mlong2 - mnode)
	    + 0.0114 * SIN8( mlong2 - mg)
            - 0.2037 * SIN8( mlong2 ) 
	    - 0.0261 * SIN8( mlong2 + mg )
	    + 0.0124 * SIN8( slong2 - mnode)
	    + 0.0214 * SIN8( slong2 - sg)
	    - 1.2729 * SIN8( slong2 ) 
	    - 0.0497 * SIN8( slong2 + sg)
            + 0.1261 * SIN8( sg ) ) / 3600.0; 
    }      
  }  /* for i */
  
  /* calculate the arguments sa[] for the disturbation terms */ 
  ti = (jd_ad - EPOCH1850) / 365.25;	/* julian years from 1850 */
  for ( i = 0; i < SDNUM; i++ ) 
    sa [i] = mod8360 (sd [i].sd0 + sd [i].sd1 * ti);
  /*
    sa[2] += 0.3315 * SIN8 (DEGTORAD *(133.9099 + 38.39365 * el[SUN].tj));
  */
  /* correction of jupiter perturbation argument for sun from Pottenger;
  creates only .03" and 1e-7 rad, not applied because origin unclear */
  thelup = jd_ad;               /* note the last helup time */
}    /* end helup() */

/******************************************************************
hel()
Computes the heliocentric positions for all planets except the moon.
The outer planets from Jupiter onwards, including Chiron, are 
actually done by a subsequent call to outer_hel() which takes
exactly the same parameters.
hel() does true position relative to the mean ecliptic and equinox
of date. Nutation is not added and must be done so by the caller.
The latitude of the Sun (max. 0.5") is neglected and always returned
as zero.

return: OK or ERR
******************************************************************/
int hel( int	planet,		/* planet index as defined by placalc.h */
	REAL8	t,		/* relative juliand date, ephemeris time */
				/* Now come 6 pointers to return values. */
	REAL8	*al,		/* longitude in degrees */
	REAL8   *ar,		/* radius in AU */
	REAL8   *az,		/* distance from ecliptic in AU */
	REAL8   *alp, 		/* speed in longitude, degrees per day */
	REAL8   *arp,		/* speed in radius, AU per day */
	REAL8   *azp)   	/* speed in z, AU per day */
{
  void  disturb();
  REAL8 fnu();
  register struct elements *e;
  register struct eledata  *d;
  REAL8 lk = 0.0;
  REAL8 rk = 0.0;
  REAL8 b,  h1,  sini, sinv, cosi, cosu, cosv, man, truanom, esquare, 
	k8,  u, up, v, vp;
      
  if (planet >= JUPITER ) 
    return ( outer_hel( planet, t, al, ar, az, alp, arp, azp ));
  if (planet < SUN || planet == MOON)
    return (ERR);

  e = &el[planet];
  d = &pd[planet];
  sini = SIN8( DEGTORAD * e->in );
  cosi = COS8( DEGTORAD * e->in );
  esquare = sqrt( ( 1.0 + e->ex ) / ( 1.0 - e->ex ) ); /* H6 in old version */
  man = e->ma;
  if ( planet == EARTH ) /* some longperiodic terms in mean longitude */
    man += (	0.266 * SIN8 ( DEGTORAD * ( 31.8 + 119.0 * e->tj ) )
	    + 6.40  * SIN8 ( DEGTORAD * ( 231.19 + 20.2 * e->tj ) )
	    + (1.882-0.016*e->tj) * SIN8( DEGTORAD * (57.24 + 150.27 * e->tj))
	   ) / 3600.0;
  if ( planet == MARS )   /* some longperiodic terms */
    man += ( 0.606 * SIN8( DEGTORAD * (212.87 + e->tj * 119.051) )
	  + 52.490 * SIN8( DEGTORAD * (47.48 + e->tj * 19.771) )
	  +  0.319 * SIN8( DEGTORAD * (116.88 + e->tj * 773.444) )
	  +  0.130 * SIN8( DEGTORAD * (74 + e->tj * 163) )
	  +  0.280 * SIN8( DEGTORAD * (300 + e->tj * 40.8) )
	  -  ( 37.05 +13.5 * e->tj )
	  ) / 3600.0;
  u = fnu ( man, e->ex, 0.0000003 ); /* error 0.001" returns radians */  
  cosu = COS8( u );
  h1 = 1 - e->ex * cosu;
  *ar = d->axis * h1;
  if ( ABS8( M_PI - u ) < TANERRLIMIT ) 
    truanom = u; /* very close to aphel */
  else
    truanom = 2.0 * ATAN8( esquare * TAN8( u * 0.5 ) ); /* true anomaly, rad*/
  v = smod8360( truanom * RADTODEG + e->pe - e->kn ); /* argument of latitude */
  if ( sini == 0.0 || ABS8( v -  90.0 ) < TANERRLIMIT 
                   || ABS8( v - 270.0 ) < TANERRLIMIT ) {
    *al = v;
  } else {
    *al = RADTODEG * ATAN8( TAN8( v * DEGTORAD ) * cosi );
    if ( v > 90.0 && v < 270.0 )  *al += 180.0;
  }
  *al = smod8360( *al + e->kn );
  sinv = SIN8( v * DEGTORAD );
  cosv = COS8( v * DEGTORAD );
  *az = *ar * sinv * sini;
  b = ASIN8( sinv * sini );   	/* latitude in radians */
  k8 = cosv / COS8( b ) * sini;
  up = 360.0 / d->period / h1;    /* du/dt degrees/day */
  if ( ABS8 ( M_PI - u ) < TANERRLIMIT ) 
    vp = up / esquare;	/* speed at aphel */
  else
    vp = up * esquare * ( 1 + COS8 ( truanom ) ) / ( 1 + cosu ); 
    /* dv/dt degrees/day */
  *arp = d->axis * up * DEGTORAD * SIN8( u ) * e->ex; 
	/* dr/dt  AU/day */
  *azp = *arp * sinv * sini + *ar * vp * DEGTORAD * cosv * sini;    /* dz/dt */
  *alp = vp / cosi * ( 1 - k8 * k8 );
  /* now come the disturbations */
  switch ( planet ) {
  REAL8 am, mma, ema, u2;
  case EARTH:  
  /* 
   * earth has some special moon values and a disturbation series due to the
   * planets. The moon stuff is due to the fact, that the mean elements
   * give the coordinates of the earth-moon barycenter. By adding the
   * corrections we effectively reduce to the center of the earth.
   * We neglect the correction in latitude, which is about 0.5", because
   * for astrological purposes we want the Sun to have latitude zero.
   */
      am = DEGTORAD * smod8360( el[MOON].lg - e->lg + 180.0 ); /* degrees */
      mma = DEGTORAD * el[MOON].ma;
      ema = DEGTORAD * e->ma;
      u2 = 2.0 * DEGTORAD * (e->lg - 180.0 - el[MOON].kn); /* 2u' */
      lk =   6.454 * SIN8( am ) 
	   + 0.013 * SIN8( 3.0 * am )
	   + 0.177 * SIN8( am + mma )
	   - 0.424 * SIN8( am - mma )
	   + 0.039 * SIN8( 3.0 * am - mma )
	   - 0.064 * SIN8( am + ema )
	   + 0.172 * SIN8( am - ema )
	   - 0.013 * SIN8( am - mma - ema)
	   - 0.013 * SIN8( u2 );
      rk =   13360 * COS8( am ) 
	   + 30    * COS8( 3.0 * am )
	   + 370   * COS8( am + mma )
	   - 1330  * COS8( am - mma )
	   + 80    * COS8( 3.0 * am - mma )
	   - 140   * COS8( am + ema )
	   + 360   * COS8( am - ema )
	   - 30    * COS8( am - mma - ema)
	   + 30    * COS8( u2 );
      /* long periodic term from mars 15g''' - 8g'', Vol 6 p19, p24 */
      lk +=  0.202 * SIN8( DEGTORAD * (315.6 + 893.3 * e->tj));
      disturb( earthkor, al, ar, lk, rk, man );
      break;
  case MERCURY:    /* only normal disturbation series */
      disturb( mercurykor, al, ar, 0.0, 0.0, man );
      break;
  case VENUS:  /* some longperiod terms and normal series */
      lk = (2.761 - 0.22*e->tj) * SIN8( DEGTORAD * (237.24 + 150.27 * e->tj))
              + 0.269 * SIN8( DEGTORAD * (212.2  + 119.05 * e->tj))
	      - 0.208 * SIN8( DEGTORAD * (175.8  + 1223.5 * e->tj));
	      /* make seconds */
      disturb( venuskor, al, ar, lk, 0.0, man );
      break;
  case MARS:   /* only normal disturbation series */
      disturb( marskor, al, ar, 0.0, 0.0, man );
      break;
  }    /* switch planet */
  return (OK);
}    /* hel */

/**************