formulas.c

占星术4.0源码

/*
** Astrolog (Version 4.00) File: formulas.c
**
** IMPORTANT NOTICE: the graphics database and chart display routines
** used in this program are Copyright (C) 1991-1993 by Walter D. Pullen
** ([e-mail]). Permission is granted to freely
** use and distribute these routines provided one doesn't sell,
** restrict, or profit from them in any way. Modification is allowed
** provided these notices remain with any altered or edited versions of
** the program.
**
** The main planetary calculation routines used in this program have
** been Copyrighted and the core of this program is basically a
** conversion to C of the routines created by James Neely as listed in
** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
** available from Matrix Software. The copyright gives us permission to
** use the routines for personal use but not to sell them or profit from
** them in any way.
**
** The PostScript code within the core graphics routines are programmed
** and Copyright (C) 1992-1993 by Brian D. Willoughby
** ([e-mail]). Conditions are identical to those above.
**
** The extended accurate ephemeris databases and formulas are from the
** calculation routines in the program "Placalc" and are programmed and
** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
** ([e-mail]). The use of that source code is subject to
** regulations made by Astrodienst Zurich, and the code is not in the
** public domain. This copyright notice must not be changed or removed
** by any user of this program.
**
** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
** X Window graphics initially programmed 10/23-29/1991.
** PostScript graphics initially programmed 11/29-30/1992.
** Last code change made 12/31/1993.
*/

#include "astrolog.h"

real MC, Asc, Vtx, OB,
  X, Y, A, R, M, D, B, Q, L, G, O, RA, S0, NU;

real planet1[TOTAL+1], planet2[TOTAL+1],
  planetalt1[TOTAL+1], planetalt2[TOTAL+1],
  house1[SIGNS+1], house2[SIGNS+1], ret1[TOTAL+1], ret2[TOTAL+1];

real FAR *datapointer;


/*
******************************************************************************
** Specific Calculations.
******************************************************************************
*/


/* Integer division - like the "/" operator but always rounds result down. */

long Dvd(x, y)
long x, y;
{
  long z;

  if (y == 0)
    return x;
  z = x / y;
  if ((x >= 0 == y >= 0) || x-z*y == 0)
    return z;
  return z - 1;
}


#ifdef MATRIX
/* Given a month, day, and year, convert it into a single Julian day value, */
/* i.e. the number of days passed since a fixed reference date.             */

long MdyToJulian(mon, day, yea)
int mon, day, yea;
{
  long im, j;

  im = 12*((long)yea+4800)+(long)mon-3;
  j = (2*(im%12) + 7 + 365*im)/12;
  j += (long)day + im/48 - 32083;
  if (j > 2299171)                   /* Take care of dates in */
    j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  return j;
}


/* Take a Julian day value, and convert it back into the corresponding */
/* month, day, and year.                                               */

void JulianToMdy(JD, mon, day, yea)
real JD;
int *mon, *day, *yea;
{
  long L, N, IT, JT, K, IK;

  L  = (long)floor(JD+ROUND)+68569L;
  N  = Dvd(4L*L, 146097L);
  L  -= Dvd(146097L*N + 3L, 4L);
  IT = Dvd(4000L*(L+1L), 1461001L);
  L  -= Dvd(1461L*IT, 4L) - 31L;
  JT = Dvd(80L*L, 2447L);
  K  = L-Dvd(2447L*JT, 80L);
  L  = Dvd(JT, 11L);
  JT += 2L - 12L*L;
  IK = 100L*(N-49L) + IT + L;
  *mon = (int)JT; *day = (int)K; *yea = (int)IK;
}


/* This is a subprocedure of CastChart(). Once we have the chart parameters, */
/* calculate a few important things related to the date, i.e. the Greenwich  */
/* time, the Julian day and fractional part of the day, the offset to the    */
/* sidereal, and a couple of other things.                                   */

real ProcessInput(var)
int var;
{
  real Off, Ln;

  TT = Sgn(TT)*floor(dabs(TT))+FRACT(dabs(TT))*100.0/60.0+DecToDeg(ZZ);
  OO = DecToDeg(OO);
  AA = DTOR(DecToDeg(AA));
  AA = MIN(AA, 89.9999);      /* Make sure the chart isn't being cast */
  AA = MAX(AA, -89.9999);     /* on the precise north or south pole.  */

  /* if parameter 'var' isn't set, then we can assume that the true time   */
  /* has already been determined (as in a -rm switch time midpoint chart). */

  if (var) {
    JD = (real)MdyToJulian(MM, DD, YY);
    if (!progress || (operation & DASHp0) > 0)
      T = ((JD-2415020.0)+TT/24.0-0.5)/36525.0;
    else

      /* Determine actual time that a progressed chart is to be cast for. */

      T = (((Jdp-JD)/progday+JD)-2415020.0+TT/24.0-0.5)/36525.0;
  }

  /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  OB = DTOR(23.452294-0.0130125*T);

  Ln = Mod((933060-6962911*T+7.5*T*T)/3600.0);    /* Mean lunar node */
  Off = (259205536.0*T+2013816.0)/3600.0;         /* Mean Sun        */
  Off = 17.23*sin(DTOR(Ln))+1.27*sin(DTOR(Off))-(5025.64+1.11*T)*T;
  Off = (Off-84038.27)/3600.0;
  SD = operation & DASHs ? Off : 0.0;
  return Off;
}


/* Convert polar to rectangular coordinates. */

void PolToRec(A, R, X, Y)
real A, R, *X, *Y;
{
  if (A == 0.0)
    A = SMALL;
  *X = R*cos(A);
  *Y = R*sin(A);
}


/* Convert rectangular to polar coordinates. */

void RecToPol(X, Y, A, R)
real X, Y, *A, *R;
{
  if (Y == 0.0)
    Y = SMALL;
  *R = sqrt(X*X+Y*Y);
  *A = Angle(X, Y);
}


/* Convert rectangular to spherical coordinates. */

void RecToSph()
{
  A = B; R = 1.0;
  PolToRec(A, R, &X, &Y);
  Q = Y; R = X; A = L;
  PolToRec(A, R, &X, &Y);
  G = X; X = Y; Y = Q;
  RecToPol(X, Y, &A, &R);
  A += O;
  PolToRec(A, R, &X, &Y);
  Q = ASIN(Y);
  Y = X; X = G;
  RecToPol(X, Y, &A, &R);
  if (A < 0.0)
    A += 2*PI;
  G = A;
}


/* Do a coordinate transformation: Given a longitude and latitude value,    */
/* return the new longitude and latitude values that the same location      */
/* would have, were the equator tilted by a specified number of degrees.    */
/* In other words, do a pole shift! This is used to convert among ecliptic, */
/* equatorial, and local coordinates, each of which have zero declination   */
/* in different planes. In other words, take into account the Earth's axis. */

void CoorXform(azi, alt, tilt)
real *azi, *alt, tilt;
{
  real x, y, a1, l1;
  real sinalt, cosalt, sinazi, sintilt, costilt;

  sinalt = sin(*alt); cosalt = cos(*alt); sinazi = sin(*azi);
  sintilt = sin(tilt); costilt = cos(tilt);

  x = cosalt * sinazi * costilt;
  y = sinalt * sintilt;
  x -= y;
  a1 = cosalt;
  y = cosalt * cos(*azi);
  l1 = Angle(y, x);
  a1 = a1 * sinazi * sintilt + sinalt * costilt;
  a1 = ASIN(a1);
  *azi = l1; *alt = a1;
}


/* This is another subprocedure of CastChart(). Calculate a few variables */
/* corresponding to the chart parameters that are used later on. The      */
/* astrological vertex (object number twenty) is also calculated here.    */

void ComputeVariables()
{
  real R2;

  RA = DTOR(Mod((6.6460656+2400.0513*T+2.58E-5*T*T+TT)*15.0-OO));
  R2 = RA; O = -OB; B = AA; A = R2; R = 1.0;
  PolToRec(A, R, &X, &Y);
  X *= cos(O);
  RecToPol(X, Y, &A, &R);
  MC = Mod(SD+RTOD(A));            /* Midheaven */
  L = R2;
  RecToSph();
#if FALSE
  AZ = Mod(SD+Mod(G+PI/2.0));      /* Ascendant */
#endif
  L= R2+PI; B = PI/2.0-dabs(B);
  if (AA < 0.0)
    B = -B;
  RecToSph();
  Vtx = Mod(SD+RTOD(G+PI/2.0));    /* Vertex */
}
#endif /* MATRIX */


/*
******************************************************************************
** House Cusp Calculations.
******************************************************************************
*/

/* This is a subprocedure of HousePlace(). Given a zodiac position, return */
/* which of the twelve houses it falls in. Remember that a special check   */
/* has to be done for the house that spans 0 degrees Aries.                */

int HousePlaceIn(point)
real point;
{
  int i = 0;

  point = Mod(point + 0.5/60.0/60.0);
  do {
    i++;
  } while (!(i >= SIGNS ||
      (point >= house[i] && point < house[Mod12(i+1)]) ||
      (house[i] > house[Mod12(i+1)] &&
      (point >= house[i] || point < house[Mod12(i+1)]))));
  return i;
}


/* For each object in the chart, determine what house it belongs in. */

void HousePlace()
{
  int i;

  for (i = 1; i <= total; i++)
    inhouse[i] = HousePlaceIn(planet[i]);
}


#ifdef MATRIX
/* The following two functions calculate the midheaven and ascendant of  */
/* the chart in question, based on time and location. They are also used */
/* in some of the house cusp calculation routines as a quick way to get  */
/* the 10th and 1st house cusps.                                         */

real CuspMidheaven()
{
  real MC;

  MC = ATAN(tan(RA)/cos(OB));
  if (MC < 0.0)
    MC += PI;
  if (RA > PI)
    MC += PI;
  return Mod(RTOD(MC)+SD);
}

real CuspAscendant()
{
  real Asc;

  Asc = Angle(-sin(RA)*cos(OB)-tan(AA)*sin(OB), cos(RA));
  return Mod(RTOD(Asc)+SD);
}


/* These are various different algorithms for calculating the house cusps: */

real CuspPlacidus(deg, FF, neg)
real deg, FF;
bool neg;
{
  int i;
  real LO, R1, R2, XS;

  R1 = RA+DTOR(deg);
  if (neg)
    X = 1.0;
  else
    X = -1.0;
  for (i = 1; i <= 10; i++) {

    /* This formula works except at 0 latitude (AA == 0.0). */

    XS = X*sin(R1)*tan(OB)*tan(AA == 0.0 ? 0.0001 : AA);
    XS = ACOS(XS);