home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource1 / ast40dos / formulas.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-01-04  |  30.4 KB  |  1,018 lines
  1. /*
  2. ** Astrolog (Version 4.00) File: formulas.c
  3. **
  4. ** IMPORTANT NOTICE: the graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1993 by Walter D. Pullen
  6. ** (cruiser1@stein.u.washington.edu). Permission is granted to freely
  7. ** use and distribute these routines provided one doesn't sell,
  8. ** restrict, or profit from them in any way. Modification is allowed
  9. ** provided these notices remain with any altered or edited versions of
  10. ** the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 12/31/1993.
  36. */
  37.  
  38. #include "astrolog.h"
  39.  
  40. real MC, Asc, Vtx, OB,
  41.   X, Y, A, R, M, D, B, Q, L, G, O, RA, S0, NU;
  42.  
  43. real planet1[TOTAL+1], planet2[TOTAL+1],
  44.   planetalt1[TOTAL+1], planetalt2[TOTAL+1],
  45.   house1[SIGNS+1], house2[SIGNS+1], ret1[TOTAL+1], ret2[TOTAL+1];
  46.  
  47. real FAR *datapointer;
  48.  
  49.  
  50. /*
  51. ******************************************************************************
  52. ** Specific Calculations.
  53. ******************************************************************************
  54. */
  55.  
  56.  
  57. /* Integer division - like the "/" operator but always rounds result down. */
  58.  
  59. long Dvd(x, y)
  60. long x, y;
  61. {
  62.   long z;
  63.  
  64.   if (y == 0)
  65.     return x;
  66.   z = x / y;
  67.   if ((x >= 0 == y >= 0) || x-z*y == 0)
  68.     return z;
  69.   return z - 1;
  70. }
  71.  
  72.  
  73. #ifdef MATRIX
  74. /* Given a month, day, and year, convert it into a single Julian day value, */
  75. /* i.e. the number of days passed since a fixed reference date.             */
  76.  
  77. long MdyToJulian(mon, day, yea)
  78. int mon, day, yea;
  79. {
  80.   long im, j;
  81.  
  82.   im = 12*((long)yea+4800)+(long)mon-3;
  83.   j = (2*(im%12) + 7 + 365*im)/12;
  84.   j += (long)day + im/48 - 32083;
  85.   if (j > 2299171)                   /* Take care of dates in */
  86.     j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  87.   return j;
  88. }
  89.  
  90.  
  91. /* Take a Julian day value, and convert it back into the corresponding */
  92. /* month, day, and year.                                               */
  93.  
  94. void JulianToMdy(JD, mon, day, yea)
  95. real JD;
  96. int *mon, *day, *yea;
  97. {
  98.   long L, N, IT, JT, K, IK;
  99.  
  100.   L  = (long)floor(JD+ROUND)+68569L;
  101.   N  = Dvd(4L*L, 146097L);
  102.   L  -= Dvd(146097L*N + 3L, 4L);
  103.   IT = Dvd(4000L*(L+1L), 1461001L);
  104.   L  -= Dvd(1461L*IT, 4L) - 31L;
  105.   JT = Dvd(80L*L, 2447L);
  106.   K  = L-Dvd(2447L*JT, 80L);
  107.   L  = Dvd(JT, 11L);
  108.   JT += 2L - 12L*L;
  109.   IK = 100L*(N-49L) + IT + L;
  110.   *mon = (int)JT; *day = (int)K; *yea = (int)IK;
  111. }
  112.  
  113.  
  114. /* This is a subprocedure of CastChart(). Once we have the chart parameters, */
  115. /* calculate a few important things related to the date, i.e. the Greenwich  */
  116. /* time, the Julian day and fractional part of the day, the offset to the    */
  117. /* sidereal, and a couple of other things.                                   */
  118.  
  119. real ProcessInput(var)
  120. int var;
  121. {
  122.   real Off, Ln;
  123.  
  124.   TT = Sgn(TT)*floor(dabs(TT))+FRACT(dabs(TT))*100.0/60.0+DecToDeg(ZZ);
  125.   OO = DecToDeg(OO);
  126.   AA = DTOR(DecToDeg(AA));
  127.   AA = MIN(AA, 89.9999);      /* Make sure the chart isn't being cast */
  128.   AA = MAX(AA, -89.9999);     /* on the precise north or south pole.  */
  129.  
  130.   /* if parameter 'var' isn't set, then we can assume that the true time   */
  131.   /* has already been determined (as in a -rm switch time midpoint chart). */
  132.  
  133.   if (var) {
  134.     JD = (real)MdyToJulian(MM, DD, YY);
  135.     if (!progress || (operation & DASHp0) > 0)
  136.       T = ((JD-2415020.0)+TT/24.0-0.5)/36525.0;
  137.     else
  138.  
  139.       /* Determine actual time that a progressed chart is to be cast for. */
  140.  
  141.       T = (((Jdp-JD)/progday+JD)-2415020.0+TT/24.0-0.5)/36525.0;
  142.   }
  143.  
  144.   /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  145.   OB = DTOR(23.452294-0.0130125*T);
  146.  
  147.   Ln = Mod((933060-6962911*T+7.5*T*T)/3600.0);    /* Mean lunar node */
  148.   Off = (259205536.0*T+2013816.0)/3600.0;         /* Mean Sun        */
  149.   Off = 17.23*sin(DTOR(Ln))+1.27*sin(DTOR(Off))-(5025.64+1.11*T)*T;
  150.   Off = (Off-84038.27)/3600.0;
  151.   SD = operation & DASHs ? Off : 0.0;
  152.   return Off;
  153. }
  154.  
  155.  
  156. /* Convert polar to rectangular coordinates. */
  157.  
  158. void PolToRec(A, R, X, Y)
  159. real A, R, *X, *Y;
  160. {
  161.   if (A == 0.0)
  162.     A = SMALL;
  163.   *X = R*cos(A);
  164.   *Y = R*sin(A);
  165. }
  166.  
  167.  
  168. /* Convert rectangular to polar coordinates. */
  169.  
  170. void RecToPol(X, Y, A, R)
  171. real X, Y, *A, *R;
  172. {
  173.   if (Y == 0.0)
  174.     Y = SMALL;
  175.   *R = sqrt(X*X+Y*Y);
  176.   *A = Angle(X, Y);
  177. }
  178.  
  179.  
  180. /* Convert rectangular to spherical coordinates. */
  181.  
  182. void RecToSph()
  183. {
  184.   A = B; R = 1.0;
  185.   PolToRec(A, R, &X, &Y);
  186.   Q = Y; R = X; A = L;
  187.   PolToRec(A, R, &X, &Y);
  188.   G = X; X = Y; Y = Q;
  189.   RecToPol(X, Y, &A, &R);
  190.   A += O;
  191.   PolToRec(A, R, &X, &Y);
  192.   Q = ASIN(Y);
  193.   Y = X; X = G;
  194.   RecToPol(X, Y, &A, &R);
  195.   if (A < 0.0)
  196.     A += 2*PI;
  197.   G = A;
  198. }
  199.  
  200.  
  201. /* Do a coordinate transformation: Given a longitude and latitude value,    */
  202. /* return the new longitude and latitude values that the same location      */
  203. /* would have, were the equator tilted by a specified number of degrees.    */
  204. /* In other words, do a pole shift! This is used to convert among ecliptic, */
  205. /* equatorial, and local coordinates, each of which have zero declination   */
  206. /* in different planes. In other words, take into account the Earth's axis. */
  207.  
  208. void CoorXform(azi, alt, tilt)
  209. real *azi, *alt, tilt;
  210. {
  211.   real x, y, a1, l1;
  212.   real sinalt, cosalt, sinazi, sintilt, costilt;
  213.  
  214.   sinalt = sin(*alt); cosalt = cos(*alt); sinazi = sin(*azi);
  215.   sintilt = sin(tilt); costilt = cos(tilt);
  216.  
  217.   x = cosalt * sinazi * costilt;
  218.   y = sinalt * sintilt;
  219.   x -= y;
  220.   a1 = cosalt;
  221.   y = cosalt * cos(*azi);
  222.   l1 = Angle(y, x);
  223.   a1 = a1 * sinazi * sintilt + sinalt * costilt;
  224.   a1 = ASIN(a1);
  225.   *azi = l1; *alt = a1;
  226. }
  227.  
  228.  
  229. /* This is another subprocedure of CastChart(). Calculate a few variables */
  230. /* corresponding to the chart parameters that are used later on. The      */
  231. /* astrological vertex (object number twenty) is also calculated here.    */
  232.  
  233. void ComputeVariables()
  234. {
  235.   real R2;
  236.  
  237.   RA = DTOR(Mod((6.6460656+2400.0513*T+2.58E-5*T*T+TT)*15.0-OO));
  238.   R2 = RA; O = -OB; B = AA; A = R2; R = 1.0;
  239.   PolToRec(A, R, &X, &Y);
  240.   X *= cos(O);
  241.   RecToPol(X, Y, &A, &R);
  242.   MC = Mod(SD+RTOD(A));            /* Midheaven */
  243.   L = R2;
  244.   RecToSph();
  245. #if FALSE
  246.   AZ = Mod(SD+Mod(G+PI/2.0));      /* Ascendant */
  247. #endif
  248.   L= R2+PI; B = PI/2.0-dabs(B);
  249.   if (AA < 0.0)
  250.     B = -B;
  251.   RecToSph();
  252.   Vtx = Mod(SD+RTOD(G+PI/2.0));    /* Vertex */
  253. }
  254. #endif /* MATRIX */
  255.  
  256.  
  257. /*
  258. ******************************************************************************
  259. ** House Cusp Calculations.
  260. ******************************************************************************
  261. */
  262.  
  263. /* This is a subprocedure of HousePlace(). Given a zodiac position, return */
  264. /* which of the twelve houses it falls in. Remember that a special check   */
  265. /* has to be done for the house that spans 0 degrees Aries.                */
  266.  
  267. int HousePlaceIn(point)
  268. real point;
  269. {
  270.   int i = 0;
  271.  
  272.   point = Mod(point + 0.5/60.0/60.0);
  273.   do {
  274.     i++;
  275.   } while (!(i >= SIGNS ||
  276.       (point >= house[i] && point < house[Mod12(i+1)]) ||
  277.       (house[i] > house[Mod12(i+1)] &&
  278.       (point >= house[i] || point < house[Mod12(i+1)]))));
  279.   return i;
  280. }
  281.  
  282.  
  283. /* For each object in the chart, determine what house it belongs in. */
  284.  
  285. void HousePlace()
  286. {
  287.   int i;
  288.  
  289.   for (i = 1; i <= total; i++)
  290.     inhouse[i] = HousePlaceIn(planet[i]);
  291. }
  292.  
  293.  
  294. #ifdef MATRIX
  295. /* The following two functions calculate the midheaven and ascendant of  */
  296. /* the chart in question, based on time and location. They are also used */
  297. /* in some of the house cusp calculation routines as a quick way to get  */
  298. /* the 10th and 1st house cusps.                                         */
  299.  
  300. real CuspMidheaven()
  301. {
  302.   real MC;
  303.  
  304.   MC = ATAN(tan(RA)/cos(OB));
  305.   if (MC < 0.0)
  306.     MC += PI;
  307.   if (RA > PI)
  308.     MC += PI;
  309.   return Mod(RTOD(MC)+SD);
  310. }
  311.  
  312. real CuspAscendant()
  313. {
  314.   real Asc;
  315.  
  316.   Asc = Angle(-sin(RA)*cos(OB)-tan(AA)*sin(OB), cos(RA));
  317.   return Mod(RTOD(Asc)+SD);
  318. }
  319.  
  320.  
  321. /* These are various different algorithms for calculating the house cusps: */
  322.  
  323. real CuspPlacidus(deg, FF, neg)
  324. real deg, FF;
  325. bool neg;
  326. {
  327.   int i;
  328.   real LO, R1, R2, XS;
  329.  
  330.   R1 = RA+DTOR(deg);
  331.   if (neg)
  332.     X = 1.0;
  333.   else
  334.     X = -1.0;
  335.   for (i = 1; i <= 10; i++) {
  336.  
  337.     /* This formula works except at 0 latitude (AA == 0.0). */
  338.  
  339.     XS = X*sin(R1)*tan(OB)*tan(AA == 0.0 ? 0.0001 : AA);
  340.     XS = ACOS(XS);
  341.     if (XS < 0.0)
  342.       XS += PI;
  343.     if (neg)
  344.       R2 = RA+PI-(XS/FF);
  345.     else
  346.       R2 = RA+(XS/FF);
  347.     R1 = R2;
  348.   }
  349.   LO = ATAN(tan(R1)/cos(OB));
  350.   if (LO < 0.0)
  351.     LO += PI;
  352.   if (sin(R1) < 0.0)
  353.     LO += PI;
  354.   return RTOD(LO);
  355. }
  356.  
  357. void HousePlacidus()
  358. {
  359.   int i;
  360.  
  361.   house[1] = Mod(Asc-SD);
  362.   house[4] = Mod(MC+DEGHALF-SD);
  363.   house[5] = CuspPlacidus(30.0, 3.0, FALSE) + DEGHALF;
  364.   house[6] = CuspPlacidus(60.0, 1.5, FALSE) + DEGHALF;
  365.   house[2] = CuspPlacidus(120.0, 1.5, TRUE);
  366.   house[3] = CuspPlacidus(150.0, 3.0, TRUE);
  367.   for (i = 1; i <= SIGNS; i++) {
  368.     if (i <= 6)
  369.       house[i] = Mod(house[i]+SD);
  370.     else
  371.       house[i] = Mod(house[i-6]+DEGHALF);
  372.   }
  373. }
  374.  
  375. void HouseKoch()
  376. {
  377.   real A1, A2, A3, KN;
  378.   int i;
  379.  
  380.   A1 = sin(RA)*tan(AA)*tan(OB);
  381.   A1 = ASIN(A1);
  382.   for (i = 1; i <= SIGNS; i++) {
  383.     D = Mod(60.0+30.0*(real)i);
  384.     A2 = D/DEGQUAD-1.0; KN = 1.0;
  385.     if (D >= DEGHALF) {
  386.       KN = -1.0;
  387.       A2 = D/DEGQUAD-3.0;
  388.     }
  389.     A3 = DTOR(Mod(RTOD(RA)+D+A2*RTOD(A1)));
  390.     X = Angle(cos(A3)*cos(OB)-KN*tan(AA)*sin(OB), sin(A3));
  391.     house[i] = Mod(RTOD(X)+SD);
  392.   }
  393. }
  394.  
  395. void HouseEqual()
  396. {
  397.   int i;
  398.  
  399.   for (i = 1; i <= SIGNS; i++)
  400.     house[i] = Mod(Asc-30.0+30.0*(real)i);
  401. }
  402.  
  403. void HouseCampanus()
  404. {
  405.   real KO, DN;
  406.   int i;
  407.  
  408.   for (i = 1; i <= SIGNS; i++) {
  409.     KO = DTOR(60.000001+30.0*(real)i);
  410.     DN = ATAN(tan(KO)*cos(AA));
  411.     if (DN < 0.0)
  412.       DN += PI;
  413.     if (sin(KO) < 0.0)
  414.       DN += PI;
  415.     X = Angle(cos(RA+DN)*cos(OB)-sin(DN)*tan(AA)*sin(OB), sin(RA+DN));
  416.     house[i] = Mod(RTOD(X)+SD);
  417.   }
  418. }
  419.  
  420. void HouseMeridian()
  421. {
  422.   int i;
  423.  
  424.   for (i = 1; i <= SIGNS; i++) {
  425.     D = DTOR(60.0+30.0*(real)i);
  426.     X = Angle(cos(RA+D)*cos(OB), sin(RA+D));
  427.     house[i] = Mod(RTOD(X)+SD);
  428.   }
  429. }
  430.  
  431. void HouseRegiomontanus()
  432. {
  433.   int i;
  434.  
  435.   for (i = 1; i <= SIGNS; i++) {
  436.     D = DTOR(60.0+30.0*(real)i);
  437.     X = Angle(cos(RA+D)*cos(OB)-sin(D)*tan(AA)*sin(OB), sin(RA+D));
  438.     house[i] = Mod(RTOD(X)+SD);
  439.   }
  440. }
  441.  
  442. void HousePorphyry()
  443. {
  444.   int i;
  445.  
  446.   X = Asc-MC;
  447.   if (X < 0.0)
  448.     X += DEGREES;
  449.   Y = X/3.0;
  450.   for (i = 1; i <= 2; i++)
  451.     house[i+4] = Mod(DEGHALF+MC+i*Y);
  452.   X = Mod(DEGHALF+MC)-Asc;
  453.   if (X < 0.0)
  454.     X += DEGREES;
  455.   house[1]=Asc;
  456.   Y = X/3.0;
  457.   for (i = 1; i <= 3; i++)
  458.     house[i+1] = Mod(Asc+i*Y);
  459.   for (i = 1; i <= 6; i++)
  460.     house[i+6] = Mod(house[i]+DEGHALF);
  461. }
  462.  
  463. void HouseMorinus()
  464. {
  465.   int i;
  466.  
  467.   for (i = 1; i <= SIGNS; i++) {
  468.     D = DTOR(60.0+30.0*(real)i);
  469.     X = Angle(cos(RA+D), sin(RA+D)*cos(OB));
  470.     house[i] = Mod(RTOD(X)+SD);
  471.   }
  472. }
  473.  
  474. real CuspTopocentric(deg)
  475. real deg;
  476. {
  477.   real OA, X, LO;
  478.  
  479.   OA = Mod(RA+DTOR(deg));
  480.   X = ATAN(tan(AA)/cos(OA));
  481.   LO = ATAN(cos(X)*tan(OA)/cos(X+OB));
  482.   if (LO < 0.0)
  483.     LO += PI;
  484.   if (sin(OA) < 0.0)
  485.     LO += PI;
  486.   return LO;
  487. }
  488.  
  489. void HouseTopocentric()
  490. {
  491.   real TL, P1, P2, LT;
  492.   int i;
  493.  
  494.   modulus = 2.0*PI;
  495.   house[4] = Mod(DTOR(MC+DEGHALF-SD));
  496.   TL = tan(AA); P1 = ATAN(TL/3.0); P2 = ATAN(TL/1.5); LT = AA;
  497.   AA = P1; house[5] = CuspTopocentric(30.0) + PI;
  498.   AA = P2; house[6] = CuspTopocentric(60.0) + PI;
  499.   AA = LT; house[1] = CuspTopocentric(90.0);
  500.   AA = P2; house[2] = CuspTopocentric(120.0);
  501.   AA = P1; house[3] = CuspTopocentric(150.0);
  502.   AA = LT; modulus = DEGREES;
  503.   for (i = 1; i <= 6; i++) {
  504.     house[i] = Mod(RTOD(house[i])+SD);
  505.     house[i+6] = Mod(house[i]+DEGHALF);
  506.   }
  507. }
  508. #endif /* MATRIX */
  509.  
  510. /* In "null" houses, the cusps are always fixed to start at their            */
  511. /* corresponding sign, i.e. the 1st house is always at 0 degrees Aries, etc. */
  512.  
  513. void HouseNull()
  514. {
  515.   int i;
  516.  
  517.   for (i = 1; i <= SIGNS; i++)
  518.     house[i] = Mod(STOZ(i)+SD);
  519. }
  520.  
  521.  
  522. /* Calculate the house cusp positions, using the specified algorithm. */
  523.  
  524. void ComputeHouses(housesystem)
  525. int housesystem;
  526. {
  527.   char string[STRING];
  528.  
  529.   if (dabs(AA) > DTOR(DEGQUAD-TROPIC) && housesystem == 0) {
  530.     sprintf(string,
  531.       "The %s system of houses is not defined at extreme latitudes.",
  532.       systemname[housesystem]);
  533.     PrintError(string);
  534.     Terminate(_FATAL);
  535.   }
  536.   switch (housesystem) {
  537.   case  1: HouseKoch();          break;
  538.   case  2: HouseEqual();         break;
  539.   case  3: HouseCampanus();      break;
  540.   case  4: HouseMeridian();      break;
  541.   case  5: HouseRegiomontanus(); break;
  542.   case  6: HousePorphyry();      break;
  543.   case  7: HouseMorinus();       break;
  544.   case  8: HouseTopocentric();   break;
  545.   case  9: HouseNull();          break;
  546.   default: HousePlacidus();
  547.   }
  548. }
  549.  
  550.  
  551. #ifdef MATRIX
  552. /*
  553. ******************************************************************************
  554. ** Planetary Position Calculations.
  555. ******************************************************************************
  556. */
  557.  
  558. /* Read the next three values from the planet data stream, and return them */
  559. /* combined as the coefficients of a quadratic equation in the chart time. */
  560.  
  561. real ReadThree()
  562. {
  563.   real S1, S2;
  564.  
  565.   S0 = ReadPlanetData(); S1 = ReadPlanetData();
  566.   S2 = ReadPlanetData();
  567.   return S0 = DTOR(S0+S1*T+S2*T*T);
  568. }
  569.  
  570.  
  571. /* Another coordinate transformation. This is used by the ComputePlanets() */
  572. /* procedure to rotate rectangular coordinates by a certain amount.        */
  573.  
  574. void RecToSph2(AP, AN, IN)
  575. real AP, AN, IN;
  576. {
  577.   RecToPol(X, Y, &A, &R); A += AP; PolToRec(A, R, &X, &Y);
  578.   D = X; X = Y; Y = 0.0; RecToPol(X, Y, &A, &R);
  579.   A += IN; PolToRec(A, R, &X, &Y);
  580.   G = Y; Y = X; X = D; RecToPol(X, Y, &A, &R); A += AN;
  581.   if (A < 0.0)
  582.     A += 2.0*PI;
  583.   PolToRec(A, R, &X, &Y);
  584. }
  585.  
  586.  
  587. /* Calculate some harmonic delta error correction factors to add onto the */
  588. /* coordinates of Jupiter through Pluto, for better accuracy.             */
  589.  
  590. void ErrorCorrect(ind, x, y, z)
  591. int ind;
  592. real *x, *y, *z;
  593. {
  594.   real U, V, W, T0[4];
  595.   int IK, IJ, errorindex;
  596.  
  597.   errorindex = errorcount[ind];
  598.   for (IK = 1; IK <= 3; IK++) {
  599.     if (ind == 6 && IK == 3) {
  600.       T0[3] = 0.0;
  601.       break;
  602.     }
  603.     if (IK == 3)
  604.       errorindex--;
  605.     ReadThree(); A = 0.0;
  606.     for (IJ = 1; IJ <= errorindex; IJ++) {
  607.       U = ReadPlanetData(); V = ReadPlanetData();
  608.       W = ReadPlanetData();
  609.       A = A+DTOR(U)*cos((V*T+W)*PI/DEGHALF);
  610.     }
  611.     T0[IK] = RTOD(S0+A);
  612.   }
  613.   *x += T0[2]; *y += T0[1]; *z += T0[3];
  614. }
  615.  
  616.  
  617. /* Another subprocedure of the ComputePlanets() routine. Convert the final */
  618. /* rectangular coordinates of a planet to zodiac position and declination. */
  619.  
  620. void ProcessPlanet(ind, aber)
  621. int ind;
  622. real aber;
  623. {
  624.   real ang, rad;
  625.  
  626.   RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  627.   planet[ind] = Mod(RTOD(ang)+NU-aber+SD);
  628.   RecToPol(rad, spacez[ind], &ang, &rad);
  629.   if (centerplanet == _SUN && ind == _SUN)
  630.     ang = 0.0;
  631.   planetalt[ind] = RTOD(ang);
  632. }
  633.  
  634.  
  635. /* This is probably the heart of the whole program of Astrolog. Calculate  */
  636. /* the position of each body that orbits the Sun. A heliocentric chart is  */
  637. /* most natural; extra calculation is needed to have other central bodies. */
  638.  
  639. void ComputePlanets()
  640. {
  641.   real helio[BASE+1], helioalt[BASE+1], helioret[BASE+1],
  642.     heliox[BASE+1], helioy[BASE+1], helioz[BASE+1];
  643.   real aber = 0.0, AU, E, EA, E1, XW, YW, AP, AN, IN, XS, YS, ZS;
  644.   int ind = 1, i;
  645.  
  646.   datapointer = planetdata;
  647.   while (ind <= (operation & DASHu ? BASE : PLANETS+1)) {
  648.     modulus = 2.0*PI;
  649.     EA = M = Mod(ReadThree());      /* Calculate mean anomaly */
  650.     E = RTOD(ReadThree());          /* Calculate eccentricity */
  651.     for (i = 1; i <= 5; i++)
  652.       EA = M+E*sin(EA);             /* Solve Kepler's equation */
  653.     AU = ReadPlanetData();          /* Semi-major axis         */
  654.     E1 = 0.01720209/(pow(AU, 1.5)*
  655.       (1.0-E*cos(EA)));                     /* Begin velocity coordinates */
  656.     XW = -AU*E1*sin(EA);                    /* Perifocal coordinates      */
  657.     YW = AU*E1*pow(1.0-E*E,0.5)*cos(EA);
  658.     AP = ReadThree(); AN = ReadThree();
  659.     IN = ReadThree();                       /* Calculate inclination       */
  660.     X = XW; Y = YW; RecToSph2(AP, AN, IN);  /* Rotate velocity coordinates */
  661.     heliox[ind] = X; helioy[ind] = Y;
  662.     helioz[ind] = G;                        /* Helio ecliptic rectangtular */
  663.     modulus = DEGREES;
  664.     X = AU*(cos(EA)-E);                 /* Perifocal coordinates for        */
  665.     Y = AU*sin(EA)*pow(1.0-E*E,0.5);    /* rectangular position coordinates */
  666.     RecToSph2(AP, AN, IN);              /* Rotate for rectangular */
  667.     XS = X; YS = Y; ZS = G;             /* position coordinates   */
  668.     if (ind >= _JUP && ind <= _PLU)
  669.       ErrorCorrect(ind, &XS, &YS, &ZS);
  670.     ret[ind] =                                        /* Helio daily motion */
  671.       (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
  672.     spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
  673.     ProcessPlanet(ind, 0.0);
  674.     ind += (ind == 1 ? 2 : (ind != PLANETS+1 ? 1 : 10));
  675.   }
  676.   spacex[0] = spacey[0] = spacez[0] = 0.0;
  677.  
  678.   if (!centerplanet) {
  679.     if (exdisplay & DASHv0)
  680.       for (i = 0; i <= BASE; i++)    /* Use relative velocity */
  681.         ret[i] = DTOR(1.0);          /* if -v0 is in effect   */
  682.     return;
  683.   }
  684. #endif /* MATRIX */
  685.  
  686.   /* A second loop is needed for geocentric charts or central bodies other */
  687.   /* than the Sun. For example, we can't find the position of Mercury in   */
  688.   /* relation to Pluto until we know the position of Pluto in relation to  */
  689.   /* the Sun, and since Mercury is calculated first, another pass needed.  */
  690.  
  691.   ind = centerplanet;
  692.   for (i = 0; i <= BASE; i++) {
  693.     helio[i]    = planet[i];
  694.     helioalt[i] = planetalt[i];
  695.     helioret[i] = ret[i];
  696.     if (i != _MOO && i != ind) {
  697.       spacex[i] -= spacex[ind];
  698.       spacey[i] -= spacey[ind];
  699.       spacez[i] -= spacez[ind];
  700.     }
  701.   }
  702.   spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  703.   SwapReal(&spacex[0], &spacex[ind]);
  704.   SwapReal(&spacey[0], &spacey[ind]);    /* Do some swapping - we want   */
  705.   SwapReal(&spacez[0], &spacez[ind]);    /* the central body to be in    */
  706.   SwapReal(&spacex[1], &spacex[ind]);    /* object position number zero. */
  707.   SwapReal(&spacey[1], &spacey[ind]);
  708.   SwapReal(&spacez[1], &spacez[ind]);
  709.   for (i = 1; i <= (operation & DASHu ? BASE : PLANETS+1);
  710.       i += (i == 1 ? 2 : (i != PLANETS+1 ? 1 : 10))) {
  711.     XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
  712.     if (ind != _SUN || i != _SUN)
  713.       ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
  714.         (XS*XS+YS*YS);
  715.     if (ind == _SUN)
  716.       aber = 0.0057756*sqrt(XS*XS+YS*YS+ZS*ZS)*RTOD(ret[i]);  /* Aberration */
  717.     ProcessPlanet(i, aber);
  718.     if (exdisplay & DASHv0)                 /* Use relative velocity */
  719.       ret[i] = DTOR(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  720.   }
  721. }
  722.  
  723.  
  724. #ifdef MATRIX
  725. /*
  726. ******************************************************************************
  727. ** Lunar Position Calculations
  728. ******************************************************************************
  729. */
  730.  
  731. /* Calculate the position and declination of the Moon, and the Moon's North  */
  732. /* Node. This has to be done separately from the other planets, because they */
  733. /* all orbit the Sun, while the Moon orbits the Earth.                       */
  734.  
  735. void ComputeLunar(moonlo, moonla, nodelo, nodela)
  736. real *moonlo, *moonla, *nodelo, *nodela;
  737. {
  738.   real LL, G, N, G1, D, L, ML, L1, MB, T1, M = 3600.0;
  739.  
  740.   LL = 973563.0+1732564379.0*T-4.0*T*T;  /* Compute mean lunar longitude     */
  741.   G = 1012395.0+6189.0*T;                /* Sun's mean longitude of perigee  */
  742.   N = 933060.0-6962911.0*T+7.5*T*T;      /* Compute mean lunar node          */
  743.   G1 = 1203586.0+14648523.0*T-37.0*T*T;  /* Mean longitude of lunar perigee  */
  744.   D = 1262655.0+1602961611.0*T-5.0*T*T;  /* Mean elongation of Moon from Sun */
  745.   L = (LL-G1)/M; L1 = ((LL-D)-G)/M;      /* Some auxiliary angles            */
  746.   T1 = (LL-N)/M; D = D/M; Y = 2.0*D;
  747.  
  748.   /* Compute Moon's perturbations. */
  749.  
  750.   ML = 22639.6*SIND(L)-4586.4*SIND(L-Y)+2369.9*SIND(Y)+769.0*SIND(2.0*L)-
  751.     669.0*SIND(L1)-411.6*SIND(2.0*T1)-212.0*SIND(2.0*L-Y)-206.0*SIND(L+L1-Y);
  752.   ML += 192.0*SIND(L+Y)-165.0*SIND(L1-Y)+148.0*SIND(L-L1)-125.0*SIND(D)-110.0*
  753.     SIND(L+L1)-55.0*SIND(2.0*T1-Y)-45.0*SIND(L+2.0*T1)+40.0*SIND(L-2.0*T1);
  754.  
  755.   *moonlo = G = Mod((LL+ML)/M+SD);       /* Lunar longitude */
  756.  
  757.   /* Compute lunar latitude. */
  758.  
  759.   MB = 18461.5*SIND(T1)+1010.0*SIND(L+T1)-999.0*SIND(T1-L)-624.0*SIND(T1-Y)+
  760.     199.0*SIND(T1+Y-L)-167.0*SIND(L+T1-Y);
  761.   MB += 117.0*SIND(T1+Y)+62.0*SIND(2.0*L+T1)-
  762.     33.0*SIND(T1-Y-L)-32.0*SIND(T1-2.0*L)-30.0*SIND(L1+T1-Y);
  763.   *moonla = MB =
  764.     Sgn(MB)*((dabs(MB)/M)/DEGREES-floor((dabs(MB)/M)/DEGREES))*DEGREES;
  765.  
  766.   /* Compute position of the north node. One can compute the true North */
  767.   /* Node here if they like, but the Mean Node seems to be the one used */
  768.   /* in astrology, so that's the one Astrolog returns by default.       */
  769.  
  770. #ifdef TRUENODE
  771.   N = N+5392.0*SIND(2.0*T1-Y)-541.0*SIND(L1)-442.0*SIND(Y)+423.0*SIND(2.0*T1)-
  772.     291.0*SIND(2.0*L-2.0*T1);
  773. #endif
  774.   *nodelo = Mod(N/M+SD);
  775.   *nodela = 0.0;
  776. }
  777. #endif /* MATRIX */
  778.  
  779.  
  780. /*
  781. ******************************************************************************
  782. ** Star Position Calculations
  783. ******************************************************************************
  784. */
  785.  
  786. /* This is used by the chart calculation routine to calculate the positions */
  787. /* of the fixed stars. Since the stars don't move in the sky over time,     */
  788. /* getting their positions is mostly just reading info from an array and    */
  789. /* converting it to the correct reference frame. However, we have to add    */
  790. /* in the correct precession for the tropical zodiac, and sort the final    */
  791. /* index list based on what order the stars are supposed to be printed in.  */
  792.  
  793. void ComputeStars(SD)
  794. real SD;
  795. {
  796.   int i, j;
  797.   real x, y, z;
  798.  
  799.   /* Read in star positions. */
  800.  
  801.   for (i = 1; i <= STARS; i++) {
  802.     x = stardata[i*6-6]; y = stardata[i*6-5]; z = stardata[i*6-4];
  803.     planet[BASE+i] = DTOR(x*DEGREES/24.0+y*15.0/60.0+z*0.25/60.0);
  804.     x = stardata[i*6-3]; y = stardata[i*6-2]; z = stardata[i*6-1];
  805.     planetalt[BASE+i] = DTOR(x+y/60.0+z/60.0/60.0);
  806.     EquToEcl(&planet[BASE+i], &planetalt[BASE+i]);           /* Convert to */
  807.     planet[BASE+i] = Mod(RTOD(planet[BASE+i])+SD2000+SD);    /* ecliptic.  */
  808.     planetalt[BASE+i] = RTOD(planetalt[BASE+i]);
  809.     starname[i] = i;
  810.   }
  811.  
  812.   /* Sort the index list if -Uz, -Ul, -Un, or -Ub switch in effect. */
  813.  
  814.   if (universe > 1) for (i = 2; i <= STARS; i++) {
  815.     j = i-1;
  816.  
  817.     /* Compare star names for -Un switch. */
  818.  
  819.     if (universe == 'n') while (j > 0 && StringCmp(
  820.       objectname[BASE+starname[j]], objectname[BASE+starname[j+1]]) > 0) {
  821.       SWAP(starname[j], starname[j+1]);
  822.       j--;
  823.  
  824.     /* Compare star brightnesses for -Ub switch. */
  825.  
  826.     } else if (universe == 'b') while (j > 0 &&
  827.       starbright[starname[j]] > starbright[starname[j+1]]) {
  828.       SWAP(starname[j], starname[j+1]);
  829.       j--;
  830.  
  831.     /* Compare star zodiac locations for -Uz switch. */
  832.  
  833.     } else if (universe == 'z') while (j > 0 &&
  834.       planet[BASE+starname[j]] > planet[BASE+starname[j+1]]) {
  835.       SWAP(starname[j], starname[j+1]);
  836.       j--;
  837.  
  838.     /* Compare star declinations for -Ul switch. */
  839.  
  840.     } else if (universe == 'l') while (j > 0 &&
  841.       planetalt[BASE+starname[j]] < planetalt[BASE+starname[j+1]]) {
  842.       SWAP(starname[j], starname[j+1]);
  843.       j--;
  844.     }
  845.   }
  846. }
  847.  
  848.  
  849. /*
  850. ******************************************************************************
  851. ** Chart Calculation.
  852. ******************************************************************************
  853. */
  854.  
  855. #ifdef PLACALC
  856. /* Compute the positions of the planets at a certain time using the Placalc */
  857. /* accurate formulas and ephemeris. This will superseed the Matrix routine  */
  858. /* values and is only called with the -b switch is in effect. Not all       */
  859. /* objects or modes are available using this, but some additional values    */
  860. /* such as Moon and Node velocities not available without -b are. (This is  */
  861. /* the one place in Astrolog which calls the Placalc package functions.)    */
  862.  
  863. void ComputePlacalc(t)
  864. real t;
  865. {
  866.   int i;
  867.   real r1, r2, r3, r4, r;
  868.  
  869.   /* We can compute the positions of Sun through Pluto, Chiron, and the */
  870.   /* North Node using Placalc. The other object must be done elsewhere. */
  871.  
  872.   for (i = _SUN; i <= _NOD; i++) if (i <= _CHI || i >= _NOD) {
  873.     if (PlacalcPlanet(i, t*36525.0+2415020.0, centerplanet != _SUN,
  874.       &r1, &r2, &r3, &r4)) {
  875.  
  876.       /* Note that this can't compute charts with central planets other */
  877.       /* than the Sun or Earth or relative velocities in current state. */
  878.  
  879.       planet[i]    = Mod(r1 + SD);
  880.       planetalt[i] = r2;
  881.       ret[i]       = DTOR(r3);
  882.  
  883.       /* Compute x,y,z coordinates from azimuth, altitude, and distance. */
  884.  
  885.       spacez[i] = r4*SIND(planetalt[i]);
  886.       r         = r4*COSD(planetalt[i]);
  887.       spacex[i] = r*COSD(planet[i]);
  888.       spacey[i] = r*SIND(planet[i]);
  889.     }
  890.   }
  891. }
  892. #endif
  893.  
  894.  
  895. /* This is probably the main routine in all of Astrolog. It generates a   */
  896. /* chart, calculating the positions of all the celestial bodies and house */
  897. /* cusps, based on the current chart information, and saves them for use  */
  898. /* by any of the display routines.                                        */
  899.  
  900. real CastChart(var)
  901. int var;
  902. {
  903.   int i, k;
  904.   real housetemp[SIGNS+1], Off = 0.0, j;
  905.  
  906.   if (MM == -1) {
  907.  
  908.     /* Hack: If month is negative, then we know chart was read in through a  */
  909.     /* -o0 position file, so the planet positions are already in the arrays. */
  910.  
  911.     MC = planet[_MC]; Asc = planet[_ASC]; Vtx = planet[_VTX];
  912.   } else {
  913.     Off = ProcessInput(var);
  914.     ComputeVariables();
  915.     if (operation & DASHG)          /* Check for -G geodetic chart. */
  916.       RA = DTOR(Mod(-OO));
  917.     MC  = CuspMidheaven();          /* Calculate our Ascendant & Midheaven. */
  918.     Asc = CuspAscendant();
  919.     ComputeHouses(housesystem);     /* Go calculate house cusps. */
  920.     for (i = 1; i <= total; i++)
  921.       ret[i] = 1.0;                 /* Assume direct until we say otherwise. */
  922.  
  923.     /* Go calculate planet, Moon, and North Node positions. */
  924.  
  925.     ComputePlanets();
  926.     ComputeLunar(&planet[_MOO], &planetalt[_MOO],
  927.       &planet[_NOD], &planetalt[_NOD]);
  928.     ret[_NOD] = -1.0;
  929.  
  930.     /* Compute more accurate ephemeris positions for certain objects. */
  931.  
  932. #ifdef PLACALC
  933.     if (placalc)
  934.       ComputePlacalc(T);
  935. #endif
  936.  
  937.     /* Calculate position of Part of Fortune. */
  938.  
  939.     j = planet[_MOO]-planet[_SUN];
  940.     j = dabs(j) < DEGQUAD ? j : j - Sgn(j)*DEGREES;
  941.     planet[_FOR] = Mod(j+Asc);
  942.  
  943.     /* Fill in "planet" positions corresponding to house cusps. */
  944.  
  945.     planet[_MC] = MC; planet[_ASC] = Asc; planet[_VTX] = Vtx;
  946.     planet[C_LO]   = house[11]; planet[C_LO+1] = house[12];
  947.     planet[C_LO+2] = house[2];  planet[C_LO+3] = house[3];
  948.   }
  949.  
  950.   /* Go calculate star positions if -U switch in effect. */
  951.  
  952.   if (universe)
  953.     ComputeStars(operation & DASHs ? 0.0 : -Off);
  954.  
  955.   /* Now, we may have to modify the base positions we calculated above based */
  956.   /* on what type of chart we are generating.                                */
  957.  
  958.   if (operation & DASHp0) {         /* Are we doing a -p0 solar arc chart?   */
  959.     for (i = 1; i <= total; i++)
  960.       planet[i] = Mod(planet[i] + (Jdp - JD) / progday);
  961.     for (i = 1; i <= SIGNS; i++)
  962.       house[i]  = Mod(house[i]  + (Jdp - JD) / progday);
  963.     }
  964.   if (multiplyfactor > 1)           /* Are we doing a -x harmonic chart?     */
  965.     for (i = 1; i <= total; i++)
  966.       planet[i] = Mod(planet[i] * (real)multiplyfactor);
  967.   if (onasc) {
  968.     if (onasc > 0)                  /* Is -1 put on Ascendant in effect?     */
  969.       j = planet[onasc]-Asc;
  970.     else                            /* Or -2 put object on Midheaven switch? */
  971.       j = planet[-onasc]-MC;
  972.     for (i = 1; i <= SIGNS; i++)    /* If so, rotate the houses accordingly. */
  973.       house[i] = Mod(house[i]+j);
  974.   }
  975.  
  976.   /* Check to see if we are -F forcing any objects to be particular values. */
  977.  
  978.   for (i = 1; i <= total; i++)
  979.     if (force[i] != 0.0) {
  980.       planet[i] = force[i]-DEGREES;
  981.       planetalt[i] = ret[i] = 0.0;
  982.     }
  983.   HousePlace();               /* Figure out what house everything falls in. */
  984.  
  985.   /* If -f domal chart switch in effect, switch planet and house positions. */
  986.  
  987.   if (operation & DASHf) {
  988.     for (i = 1; i <= total; i++) {
  989.       k = inhouse[i];
  990.       inhouse[i] = ZTOS(planet[i]);
  991.       planet[i] = STOZ(k)+MinDistance(house[k], planet[i]) /
  992.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  993.     }
  994.     for (i = 1; i <= SIGNS; i++) {
  995.       k = HousePlaceIn(STOZ(i));
  996.       housetemp[i] = STOZ(k)+MinDistance(house[k], STOZ(i)) /
  997.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  998.     }
  999.     for (i = 1; i <= SIGNS; i++)
  1000.       house[i] = housetemp[i];
  1001.   }
  1002.  
  1003.   /* If -3 decan chart switch in effect, edit planet positions accordingly. */
  1004.  
  1005.   if (operation & DASH3)
  1006.     for (i = 1; i <= total; i++) {
  1007.       k = ZTOS(planet[i]);
  1008.       j = planet[i] - STOZ(k);
  1009.       k = Mod12(k + 4*((int)floor(j/10.0)));
  1010.       j = (j - floor(j/10.0)*10.0)*3.0;
  1011.       planet[i] = STOZ(k)+j;
  1012.       HousePlace();
  1013.     }
  1014.   return T;
  1015. }
  1016.  
  1017. /* formulas.c */
  1018.