formulas.c

占星术4.0源码

    if (XS < 0.0)
      XS += PI;
    if (neg)
      R2 = RA+PI-(XS/FF);
    else
      R2 = RA+(XS/FF);
    R1 = R2;
  }
  LO = ATAN(tan(R1)/cos(OB));
  if (LO < 0.0)
    LO += PI;
  if (sin(R1) < 0.0)
    LO += PI;
  return RTOD(LO);
}

void HousePlacidus()
{
  int i;

  house[1] = Mod(Asc-SD);
  house[4] = Mod(MC+DEGHALF-SD);
  house[5] = CuspPlacidus(30.0, 3.0, FALSE) + DEGHALF;
  house[6] = CuspPlacidus(60.0, 1.5, FALSE) + DEGHALF;
  house[2] = CuspPlacidus(120.0, 1.5, TRUE);
  house[3] = CuspPlacidus(150.0, 3.0, TRUE);
  for (i = 1; i <= SIGNS; i++) {
    if (i <= 6)
      house[i] = Mod(house[i]+SD);
    else
      house[i] = Mod(house[i-6]+DEGHALF);
  }
}

void HouseKoch()
{
  real A1, A2, A3, KN;
  int i;

  A1 = sin(RA)*tan(AA)*tan(OB);
  A1 = ASIN(A1);
  for (i = 1; i <= SIGNS; i++) {
    D = Mod(60.0+30.0*(real)i);
    A2 = D/DEGQUAD-1.0; KN = 1.0;
    if (D >= DEGHALF) {
      KN = -1.0;
      A2 = D/DEGQUAD-3.0;
    }
    A3 = DTOR(Mod(RTOD(RA)+D+A2*RTOD(A1)));
    X = Angle(cos(A3)*cos(OB)-KN*tan(AA)*sin(OB), sin(A3));
    house[i] = Mod(RTOD(X)+SD);
  }
}

void HouseEqual()
{
  int i;

  for (i = 1; i <= SIGNS; i++)
    house[i] = Mod(Asc-30.0+30.0*(real)i);
}

void HouseCampanus()
{
  real KO, DN;
  int i;

  for (i = 1; i <= SIGNS; i++) {
    KO = DTOR(60.000001+30.0*(real)i);
    DN = ATAN(tan(KO)*cos(AA));
    if (DN < 0.0)
      DN += PI;
    if (sin(KO) < 0.0)
      DN += PI;
    X = Angle(cos(RA+DN)*cos(OB)-sin(DN)*tan(AA)*sin(OB), sin(RA+DN));
    house[i] = Mod(RTOD(X)+SD);
  }
}

void HouseMeridian()
{
  int i;

  for (i = 1; i <= SIGNS; i++) {
    D = DTOR(60.0+30.0*(real)i);
    X = Angle(cos(RA+D)*cos(OB), sin(RA+D));
    house[i] = Mod(RTOD(X)+SD);
  }
}

void HouseRegiomontanus()
{
  int i;

  for (i = 1; i <= SIGNS; i++) {
    D = DTOR(60.0+30.0*(real)i);
    X = Angle(cos(RA+D)*cos(OB)-sin(D)*tan(AA)*sin(OB), sin(RA+D));
    house[i] = Mod(RTOD(X)+SD);
  }
}

void HousePorphyry()
{
  int i;

  X = Asc-MC;
  if (X < 0.0)
    X += DEGREES;
  Y = X/3.0;
  for (i = 1; i <= 2; i++)
    house[i+4] = Mod(DEGHALF+MC+i*Y);
  X = Mod(DEGHALF+MC)-Asc;
  if (X < 0.0)
    X += DEGREES;
  house[1]=Asc;
  Y = X/3.0;
  for (i = 1; i <= 3; i++)
    house[i+1] = Mod(Asc+i*Y);
  for (i = 1; i <= 6; i++)
    house[i+6] = Mod(house[i]+DEGHALF);
}

void HouseMorinus()
{
  int i;

  for (i = 1; i <= SIGNS; i++) {
    D = DTOR(60.0+30.0*(real)i);
    X = Angle(cos(RA+D), sin(RA+D)*cos(OB));
    house[i] = Mod(RTOD(X)+SD);
  }
}

real CuspTopocentric(deg)
real deg;
{
  real OA, X, LO;

  OA = Mod(RA+DTOR(deg));
  X = ATAN(tan(AA)/cos(OA));
  LO = ATAN(cos(X)*tan(OA)/cos(X+OB));
  if (LO < 0.0)
    LO += PI;
  if (sin(OA) < 0.0)
    LO += PI;
  return LO;
}

void HouseTopocentric()
{
  real TL, P1, P2, LT;
  int i;

  modulus = 2.0*PI;
  house[4] = Mod(DTOR(MC+DEGHALF-SD));
  TL = tan(AA); P1 = ATAN(TL/3.0); P2 = ATAN(TL/1.5); LT = AA;
  AA = P1; house[5] = CuspTopocentric(30.0) + PI;
  AA = P2; house[6] = CuspTopocentric(60.0) + PI;
  AA = LT; house[1] = CuspTopocentric(90.0);
  AA = P2; house[2] = CuspTopocentric(120.0);
  AA = P1; house[3] = CuspTopocentric(150.0);
  AA = LT; modulus = DEGREES;
  for (i = 1; i <= 6; i++) {
    house[i] = Mod(RTOD(house[i])+SD);
    house[i+6] = Mod(house[i]+DEGHALF);
  }
}
#endif /* MATRIX */

/* In "null" houses, the cusps are always fixed to start at their            */
/* corresponding sign, i.e. the 1st house is always at 0 degrees Aries, etc. */

void HouseNull()
{
  int i;

  for (i = 1; i <= SIGNS; i++)
    house[i] = Mod(STOZ(i)+SD);
}


/* Calculate the house cusp positions, using the specified algorithm. */

void ComputeHouses(housesystem)
int housesystem;
{
  char string[STRING];

  if (dabs(AA) > DTOR(DEGQUAD-TROPIC) && housesystem == 0) {
    sprintf(string,
      "The %s system of houses is not defined at extreme latitudes.",
      systemname[housesystem]);
    PrintError(string);
    Terminate(_FATAL);
  }
  switch (housesystem) {
  case  1: HouseKoch();          break;
  case  2: HouseEqual();         break;
  case  3: HouseCampanus();      break;
  case  4: HouseMeridian();      break;
  case  5: HouseRegiomontanus(); break;
  case  6: HousePorphyry();      break;
  case  7: HouseMorinus();       break;
  case  8: HouseTopocentric();   break;
  case  9: HouseNull();          break;
  default: HousePlacidus();
  }
}


#ifdef MATRIX
/*
******************************************************************************
** Planetary Position Calculations.
******************************************************************************
*/

/* Read the next three values from the planet data stream, and return them */
/* combined as the coefficients of a quadratic equation in the chart time. */

real ReadThree()
{
  real S1, S2;

  S0 = ReadPlanetData(); S1 = ReadPlanetData();
  S2 = ReadPlanetData();
  return S0 = DTOR(S0+S1*T+S2*T*T);
}


/* Another coordinate transformation. This is used by the ComputePlanets() */
/* procedure to rotate rectangular coordinates by a certain amount.        */

void RecToSph2(AP, AN, IN)
real AP, AN, IN;
{
  RecToPol(X, Y, &A, &R); A += AP; PolToRec(A, R, &X, &Y);
  D = X; X = Y; Y = 0.0; RecToPol(X, Y, &A, &R);
  A += IN; PolToRec(A, R, &X, &Y);
  G = Y; Y = X; X = D; RecToPol(X, Y, &A, &R); A += AN;
  if (A < 0.0)
    A += 2.0*PI;
  PolToRec(A, R, &X, &Y);
}


/* Calculate some harmonic delta error correction factors to add onto the */
/* coordinates of Jupiter through Pluto, for better accuracy.             */

void ErrorCorrect(ind, x, y, z)
int ind;
real *x, *y, *z;
{
  real U, V, W, T0[4];
  int IK, IJ, errorindex;

  errorindex = errorcount[ind];
  for (IK = 1; IK <= 3; IK++) {
    if (ind == 6 && IK == 3) {
      T0[3] = 0.0;
      break;
    }
    if (IK == 3)
      errorindex--;
    ReadThree(); A = 0.0;
    for (IJ = 1; IJ <= errorindex; IJ++) {
      U = ReadPlanetData(); V = ReadPlanetData();
      W = ReadPlanetData();
      A = A+DTOR(U)*cos((V*T+W)*PI/DEGHALF);
    }
    T0[IK] = RTOD(S0+A);
  }
  *x += T0[2]; *y += T0[1]; *z += T0[3];
}


/* Another subprocedure of the ComputePlanets() routine. Convert the final */
/* rectangular coordinates of a planet to zodiac position and declination. */

void ProcessPlanet(ind, aber)
int ind;
real aber;
{
  real ang, rad;

  RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  planet[ind] = Mod(RTOD(ang)+NU-aber+SD);
  RecToPol(rad, spacez[ind], &ang, &rad);
  if (centerplanet == _SUN && ind == _SUN)
    ang = 0.0;
  planetalt[ind] = RTOD(ang);
}


/* This is probably the heart of the whole program of Astrolog. Calculate  */
/* the position of each body that orbits the Sun. A heliocentric chart is  */
/* most natural; extra calculation is needed to have other central bodies. */

void ComputePlanets()
{
  real helio[BASE+1], helioalt[BASE+1], helioret[BASE+1],
    heliox[BASE+1], helioy[BASE+1], helioz[BASE+1];
  real aber = 0.0, AU, E, EA, E1, XW, YW, AP, AN, IN, XS, YS, ZS;
  int ind = 1, i;

  datapointer = planetdata;
  while (ind <= (operation & DASHu ? BASE : PLANETS+1)) {
    modulus = 2.0*PI;
    EA = M = Mod(ReadThree());      /* Calculate mean anomaly */
    E = RTOD(ReadThree());          /* Calculate eccentricity */
    for (i = 1; i <= 5; i++)
      EA = M+E*sin(EA);             /* Solve Kepler's equation */
    AU = ReadPlanetData();          /* Semi-major axis         */
    E1 = 0.01720209/(pow(AU, 1.5)*
      (1.0-E*cos(EA)));                     /* Begin velocity coordinates */
    XW = -AU*E1*sin(EA);                    /* Perifocal coordinates      */
    YW = AU*E1*pow(1.0-E*E,0.5)*cos(EA);
    AP = ReadThree(); AN = ReadThree();
    IN = ReadThree();                       /* Calculate inclination       */
    X = XW; Y = YW; RecToSph2(AP, AN, IN);  /* Rotate velocity coordinates */
    heliox[ind] = X; helioy[ind] = Y;
    helioz[ind] = G;                        /* Helio ecliptic rectangtular */
    modulus = DEGREES;
    X = AU*(cos(EA)-E);                 /* Perifocal coordinates for        */
    Y = AU*sin(EA)*pow(1.0-E*E,0.5);    /* rectangular position coordinates */
    RecToSph2(AP, AN, IN);              /* Rotate for rectangular */
    XS = X; YS = Y; ZS = G;             /* position coordinates   */
    if (ind >= _JUP && ind <= _PLU)
      ErrorCorrect(ind, &XS, &YS, &ZS);
    ret[ind] =                                        /* Helio daily motion */
      (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
    spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
    ProcessPlanet(ind, 0.0);
    ind += (ind == 1 ? 2 : (ind != PLANETS+1 ? 1 : 10));
  }
  spacex[0] = spacey[0] = spacez[0] = 0.0;

  if (!centerplanet) {
    if (exdisplay & DASHv0)
      for (i = 0; i <= BASE; i++)    /* Use relative velocity */