formulas.c

占星术4.0源码

        ret[i] = DTOR(1.0);          /* if -v0 is in effect   */
    return;
  }
#endif /* MATRIX */

  /* A second loop is needed for geocentric charts or central bodies other */
  /* than the Sun. For example, we can't find the position of Mercury in   */
  /* relation to Pluto until we know the position of Pluto in relation to  */
  /* the Sun, and since Mercury is calculated first, another pass needed.  */

  ind = centerplanet;
  for (i = 0; i <= BASE; i++) {
    helio[i]    = planet[i];
    helioalt[i] = planetalt[i];
    helioret[i] = ret[i];
    if (i != _MOO && i != ind) {
      spacex[i] -= spacex[ind];
      spacey[i] -= spacey[ind];
      spacez[i] -= spacez[ind];
    }
  }
  spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  SwapReal(&spacex[0], &spacex[ind]);
  SwapReal(&spacey[0], &spacey[ind]);    /* Do some swapping - we want   */
  SwapReal(&spacez[0], &spacez[ind]);    /* the central body to be in    */
  SwapReal(&spacex[1], &spacex[ind]);    /* object position number zero. */
  SwapReal(&spacey[1], &spacey[ind]);
  SwapReal(&spacez[1], &spacez[ind]);
  for (i = 1; i <= (operation & DASHu ? BASE : PLANETS+1);
      i += (i == 1 ? 2 : (i != PLANETS+1 ? 1 : 10))) {
    XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
    if (ind != _SUN || i != _SUN)
      ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
        (XS*XS+YS*YS);
    if (ind == _SUN)
      aber = 0.0057756*sqrt(XS*XS+YS*YS+ZS*ZS)*RTOD(ret[i]);  /* Aberration */
    ProcessPlanet(i, aber);
    if (exdisplay & DASHv0)                 /* Use relative velocity */
      ret[i] = DTOR(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  }
}


#ifdef MATRIX
/*
******************************************************************************
** Lunar Position Calculations
******************************************************************************
*/

/* Calculate the position and declination of the Moon, and the Moon's North  */
/* Node. This has to be done separately from the other planets, because they */
/* all orbit the Sun, while the Moon orbits the Earth.                       */

void ComputeLunar(moonlo, moonla, nodelo, nodela)
real *moonlo, *moonla, *nodelo, *nodela;
{
  real LL, G, N, G1, D, L, ML, L1, MB, T1, M = 3600.0;

  LL = 973563.0+1732564379.0*T-4.0*T*T;  /* Compute mean lunar longitude     */
  G = 1012395.0+6189.0*T;                /* Sun's mean longitude of perigee  */
  N = 933060.0-6962911.0*T+7.5*T*T;      /* Compute mean lunar node          */
  G1 = 1203586.0+14648523.0*T-37.0*T*T;  /* Mean longitude of lunar perigee  */
  D = 1262655.0+1602961611.0*T-5.0*T*T;  /* Mean elongation of Moon from Sun */
  L = (LL-G1)/M; L1 = ((LL-D)-G)/M;      /* Some auxiliary angles            */
  T1 = (LL-N)/M; D = D/M; Y = 2.0*D;

  /* Compute Moon's perturbations. */

  ML = 22639.6*SIND(L)-4586.4*SIND(L-Y)+2369.9*SIND(Y)+769.0*SIND(2.0*L)-
    669.0*SIND(L1)-411.6*SIND(2.0*T1)-212.0*SIND(2.0*L-Y)-206.0*SIND(L+L1-Y);
  ML += 192.0*SIND(L+Y)-165.0*SIND(L1-Y)+148.0*SIND(L-L1)-125.0*SIND(D)-110.0*
    SIND(L+L1)-55.0*SIND(2.0*T1-Y)-45.0*SIND(L+2.0*T1)+40.0*SIND(L-2.0*T1);

  *moonlo = G = Mod((LL+ML)/M+SD);       /* Lunar longitude */

  /* Compute lunar latitude. */

  MB = 18461.5*SIND(T1)+1010.0*SIND(L+T1)-999.0*SIND(T1-L)-624.0*SIND(T1-Y)+
    199.0*SIND(T1+Y-L)-167.0*SIND(L+T1-Y);
  MB += 117.0*SIND(T1+Y)+62.0*SIND(2.0*L+T1)-
    33.0*SIND(T1-Y-L)-32.0*SIND(T1-2.0*L)-30.0*SIND(L1+T1-Y);
  *moonla = MB =
    Sgn(MB)*((dabs(MB)/M)/DEGREES-floor((dabs(MB)/M)/DEGREES))*DEGREES;

  /* Compute position of the north node. One can compute the true North */
  /* Node here if they like, but the Mean Node seems to be the one used */
  /* in astrology, so that's the one Astrolog returns by default.       */

#ifdef TRUENODE
  N = N+5392.0*SIND(2.0*T1-Y)-541.0*SIND(L1)-442.0*SIND(Y)+423.0*SIND(2.0*T1)-
    291.0*SIND(2.0*L-2.0*T1);
#endif
  *nodelo = Mod(N/M+SD);
  *nodela = 0.0;
}
#endif /* MATRIX */


/*
******************************************************************************
** Star Position Calculations
******************************************************************************
*/

/* This is used by the chart calculation routine to calculate the positions */
/* of the fixed stars. Since the stars don't move in the sky over time,     */
/* getting their positions is mostly just reading info from an array and    */
/* converting it to the correct reference frame. However, we have to add    */
/* in the correct precession for the tropical zodiac, and sort the final    */
/* index list based on what order the stars are supposed to be printed in.  */

void ComputeStars(SD)
real SD;
{
  int i, j;
  real x, y, z;

  /* Read in star positions. */

  for (i = 1; i <= STARS; i++) {
    x = stardata[i*6-6]; y = stardata[i*6-5]; z = stardata[i*6-4];
    planet[BASE+i] = DTOR(x*DEGREES/24.0+y*15.0/60.0+z*0.25/60.0);
    x = stardata[i*6-3]; y = stardata[i*6-2]; z = stardata[i*6-1];
    planetalt[BASE+i] = DTOR(x+y/60.0+z/60.0/60.0);
    EquToEcl(&planet[BASE+i], &planetalt[BASE+i]);           /* Convert to */
    planet[BASE+i] = Mod(RTOD(planet[BASE+i])+SD2000+SD);    /* ecliptic.  */
    planetalt[BASE+i] = RTOD(planetalt[BASE+i]);
    starname[i] = i;
  }

  /* Sort the index list if -Uz, -Ul, -Un, or -Ub switch in effect. */

  if (universe > 1) for (i = 2; i <= STARS; i++) {
    j = i-1;

    /* Compare star names for -Un switch. */

    if (universe == 'n') while (j > 0 && StringCmp(
      objectname[BASE+starname[j]], objectname[BASE+starname[j+1]]) > 0) {
      SWAP(starname[j], starname[j+1]);
      j--;

    /* Compare star brightnesses for -Ub switch. */

    } else if (universe == 'b') while (j > 0 &&
      starbright[starname[j]] > starbright[starname[j+1]]) {
      SWAP(starname[j], starname[j+1]);
      j--;

    /* Compare star zodiac locations for -Uz switch. */

    } else if (universe == 'z') while (j > 0 &&
      planet[BASE+starname[j]] > planet[BASE+starname[j+1]]) {
      SWAP(starname[j], starname[j+1]);
      j--;

    /* Compare star declinations for -Ul switch. */

    } else if (universe == 'l') while (j > 0 &&
      planetalt[BASE+starname[j]] < planetalt[BASE+starname[j+1]]) {
      SWAP(starname[j], starname[j+1]);
      j--;
    }
  }
}


/*
******************************************************************************
** Chart Calculation.
******************************************************************************
*/

#ifdef PLACALC
/* Compute the positions of the planets at a certain time using the Placalc */
/* accurate formulas and ephemeris. This will superseed the Matrix routine  */
/* values and is only called with the -b switch is in effect. Not all       */
/* objects or modes are available using this, but some additional values    */
/* such as Moon and Node velocities not available without -b are. (This is  */
/* the one place in Astrolog which calls the Placalc package functions.)    */

void ComputePlacalc(t)
real t;
{
  int i;
  real r1, r2, r3, r4, r;

  /* We can compute the positions of Sun through Pluto, Chiron, and the */
  /* North Node using Placalc. The other object must be done elsewhere. */

  for (i = _SUN; i <= _NOD; i++) if (i <= _CHI || i >= _NOD) {
    if (PlacalcPlanet(i, t*36525.0+2415020.0, centerplanet != _SUN,
      &r1, &r2, &r3, &r4)) {

      /* Note that this can't compute charts with central planets other */
      /* than the Sun or Earth or relative velocities in current state. */

      planet[i]    = Mod(r1 + SD);
      planetalt[i] = r2;
      ret[i]       = DTOR(r3);

      /* Compute x,y,z coordinates from azimuth, altitude, and distance. */

      spacez[i] = r4*SIND(planetalt[i]);
      r         = r4*COSD(planetalt[i]);
      spacex[i] = r*COSD(planet[i]);
      spacey[i] = r*SIND(planet[i]);
    }
  }
}
#endif


/* This is probably the main routine in all of Astrolog. It generates a   */
/* chart, calculating the positions of all the celestial bodies and house */
/* cusps, based on the current chart information, and saves them for use  */
/* by any of the display routines.                                        */

real CastChart(var)
int var;
{
  int i, k;
  real housetemp[SIGNS+1], Off = 0.0, j;

  if (MM == -1) {

    /* Hack: If month is negative, then we know chart was read in through a  */
    /* -o0 position file, so the planet positions are already in the arrays. */

    MC = planet[_MC]; Asc = planet[_ASC]; Vtx = planet[_VTX];
  } else {
    Off = ProcessInput(var);
    ComputeVariables();
    if (operation & DASHG)          /* Check for -G geodetic chart. */
      RA = DTOR(Mod(-OO));
    MC  = CuspMidheaven();          /* Calculate our Ascendant & Midheaven. */
    Asc = CuspAscendant();
    ComputeHouses(housesystem);     /* Go calculate house cusps. */
    for (i = 1; i <= total; i++)
      ret[i] = 1.0;                 /* Assume direct until we say otherwise. */

    /* Go calculate planet, Moon, and North Node positions. */

    ComputePlanets();
    ComputeLunar(&planet[_MOO], &planetalt[_MOO],
      &planet[_NOD], &planetalt[_NOD]);
    ret[_NOD] = -1.0;

    /* Compute more accurate ephemeris positions for certain objects. */

#ifdef PLACALC
    if (placalc)
      ComputePlacalc(T);
#endif

    /* Calculate position of Part of Fortune. */

    j = planet[_MOO]-planet[_SUN];
    j = dabs(j) < DEGQUAD ? j : j - Sgn(j)*DEGREES;
    planet[_FOR] = Mod(j+Asc);

    /* Fill in "planet" positions corresponding to house cusps. */

    planet[_MC] = MC; planet[_ASC] = Asc; planet[_VTX] = Vtx;
    planet[C_LO]   = house[11]; planet[C_LO+1] = house[12];
    planet[C_LO+2] = house[2];  planet[C_LO+3] = house[3];
  }

  /* Go calculate star positions if -U switch in effect. */

  if (universe)
    ComputeStars(operation & DASHs ? 0.0 : -Off);

  /* Now, we may have to modify the base positions we calculated above based */
  /* on what type of chart we are generating.                                */

  if (operation & DASHp0) {         /* Are we doing a -p0 solar arc chart?   */
    for (i = 1; i <= total; i++)
      planet[i] = Mod(planet[i] + (Jdp - JD) / progday);
    for (i = 1; i <= SIGNS; i++)
      house[i]  = Mod(house[i]  + (Jdp - JD) / progday);
    }
  if (multiplyfactor > 1)           /* Are we doing a -x harmonic chart?     */
    for (i = 1; i <= total; i++)
      planet[i] = Mod(planet[i] * (real)multiplyfactor);
  if (onasc) {
    if (onasc > 0)                  /* Is -1 put on Ascendant in effect?     */
      j = planet[onasc]-Asc;
    else                            /* Or -2 put object on Midheaven switch? */
      j = planet[-onasc]-MC;
    for (i = 1; i <= SIGNS; i++)    /* If so, rotate the houses accordingly. */
      house[i] = Mod(house[i]+j);
  }

  /* Check to see if we are -F forcing any objects to be particular values. */

  for (i = 1; i <= total; i++)
    if (force[i] != 0.0) {
      planet[i] = force[i]-DEGREES;
      planetalt[i] = ret[i] = 0.0;
    }
  HousePlace();               /* Figure out what house everything falls in. */

  /* If -f domal chart switch in effect, switch planet and house positions. */

  if (operation & DASHf) {
    for (i = 1; i <= total; i++) {
      k = inhouse[i];
      inhouse[i] = ZTOS(planet[i]);
      planet[i] = STOZ(k)+MinDistance(house[k], planet[i]) /
        MinDistance(house[k], house[Mod12(k+1)])*30.0;
    }
    for (i = 1; i <= SIGNS; i++) {
      k = HousePlaceIn(STOZ(i));
      housetemp[i] = STOZ(k)+MinDistance(house[k], STOZ(i)) /
        MinDistance(house[k], house[Mod12(k+1)])*30.0;
    }
    for (i = 1; i <= SIGNS; i++)
      house[i] = housetemp[i];
  }

  /* If -3 decan chart switch in effect, edit planet positions accordingly. */

  if (operation & DASH3)
    for (i = 1; i <= total; i++) {
      k = ZTOS(planet[i]);
      j = planet[i] - STOZ(k);
      k = Mod12(k + 4*((int)floor(j/10.0)));
      j = (j - floor(j/10.0)*10.0)*3.0;
      planet[i] = STOZ(k)+j;
      HousePlace();
    }
  return T;
}

/* formulas.c */