1
2 /*******************************************************
3 module swehouse.c
4 house and (simple) aspect calculation
5
6 ************************************************************/
7 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
8
9 License conditions
10 ------------------
11
12 This file is part of Swiss Ephemeris.
13
14 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
15 or distributor accepts any responsibility for the consequences of using it,
16 or for whether it serves any particular purpose or works at all, unless he
17 or she says so in writing.
18
19 Swiss Ephemeris is made available by its authors under a dual licensing
20 system. The software developer, who uses any part of Swiss Ephemeris
21 in his or her software, must choose between one of the two license models,
22 which are
23 a) GNU Affero General Public License (AGPL)
24 b) Swiss Ephemeris Professional License
25
26 The choice must be made before the software developer distributes software
27 containing parts of Swiss Ephemeris to others, and before any public
28 service using the developed software is activated.
29
30 If the developer choses the AGPL software license, he or she must fulfill
31 the conditions of that license, which includes the obligation to place his
32 or her whole software project under the AGPL or a compatible license.
33 See https://www.gnu.org/licenses/agpl-3.0.html
34
35 If the developer choses the Swiss Ephemeris Professional license,
36 he must follow the instructions as found in http://www.astro.com/swisseph/
37 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
38 and sign the corresponding license contract.
39
40 The License grants you the right to use, copy, modify and redistribute
41 Swiss Ephemeris, but only under certain conditions described in the License.
42 Among other things, the License requires that the copyright notices and
43 this notice be preserved on all copies.
44
45 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
46
47 The authors of Swiss Ephemeris have no control or influence over any of
48 the derived works, i.e. over software or services created by other
49 programmers which use Swiss Ephemeris functions.
50
51 The names of the authors or of the copyright holder (Astrodienst) must not
52 be used for promoting any software, product or service which uses or contains
53 the Swiss Ephemeris. This copyright notice is the ONLY place where the
54 names of the authors can legally appear, except in cases where they have
55 given special permission in writing.
56
57 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
58 for promoting such software, products or services.
59 */
60
61 //#include "sweodef.h"
62 #include "swephexp.h"
63 #include "sweph.h"
64 #include "swephlib.h"
65 #include "swehouse.h"
66 #include <string.h>
67
68 #define MILLIARCSEC (1.0 / 3600000.0)
69 #define SOLAR_YEAR 365.24219893
70 #define ARMCS ((SOLAR_YEAR+1) / SOLAR_YEAR * 360)
71
72 static double Asc1(double, double, double, double);
73 static double AscDash(double, double, double, double);
74 static double Asc2(double, double, double, double);
75 static int CalcH(double th, double fi, double ekl, char hsy, struct houses *hsp);
76 static int sidereal_houses_ecl_t0(double tjde,
77 double armc,
78 double eps,
79 double *nutlo,
80 double lat,
81 int hsys,
82 double *cusp,
83 double *ascmc,
84 double *cusp_speed,
85 double *ascmc_speed,
86 char *serr);
87 static int sidereal_houses_trad(double tjde,
88 int32 iflag,
89 double armc,
90 double eps,
91 double nutl,
92 double lat,
93 int hsys,
94 double *cusp,
95 double *ascmc,
96 double *cusp_speed,
97 double *ascmc_speed,
98 char *serr);
99 static int sidereal_houses_ssypl(double tjde,
100 double armc,
101 double eps,
102 double *nutlo,
103 double lat,
104 int hsys,
105 double *cusp,
106 double *ascmc,
107 double *cusp_speed,
108 double *ascmc_speed,
109 char *serr);
110 static int sunshine_solution_makransky(double ramc, double lat, double ecl, struct houses *hsp);
111 static int sunshine_solution_treindl(double ramc, double lat, double ecl, struct houses *hsp);
112 #if 0
113 static void test_Asc1();
114 #endif
115
116 /* housasp.c
117 * cusps are returned in double cusp[13],
118 * or cusp[37] with house system 'G'.
119 * cusp[1...12] houses 1 - 12
120 * additional points are returned in ascmc[10].
121 * ascmc[0] = ascendant
122 * ascmc[1] = mc
123 * ascmc[2] = armc
124 * ascmc[3] = vertex
125 * ascmc[4] = equasc * "equatorial ascendant" *
126 * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) *
127 * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) *
128 * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) *
129 */
swe_houses(double tjd_ut,double geolat,double geolon,int hsys,double * cusp,double * ascmc)130 int CALL_CONV swe_houses(double tjd_ut,
131 double geolat,
132 double geolon,
133 int hsys,
134 double *cusp,
135 double *ascmc)
136 {
137 int i, retc = 0;
138 double armc, eps, nutlo[2];
139 double tjde = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
140 eps = swi_epsiln(tjde, 0) * RADTODEG;
141 swi_nutation(tjde, 0, nutlo);
142 for (i = 0; i < 2; i++)
143 nutlo[i] *= RADTODEG;
144 armc = swe_degnorm(swe_sidtime0(tjd_ut, eps + nutlo[1], nutlo[0]) * 15 + geolon);
145 if (toupper(hsys) == 'I') { // compute sun declination for sunshine houses
146 int flags = SEFLG_SPEED| SEFLG_EQUATORIAL;
147 double xp[6];
148 int result = swe_calc_ut(tjd_ut, SE_SUN, flags, xp, NULL);
149 if (result < 0) {
150 // in case of failure, Porphyry houses
151 result = swe_houses_armc_ex2(armc, geolat, eps + nutlo[1], 'O', cusp, ascmc, NULL, NULL, NULL);
152 return ERR;
153 }
154 ascmc[9] = xp[1]; // declination in ascmc[9];
155 }
156 #ifdef TRACE
157 swi_open_trace(NULL);
158 if (swi_trace_count <= TRACE_COUNT_MAX) {
159 if (swi_fp_trace_c != NULL) {
160 fputs("\n/*SWE_HOUSES*/\n", swi_fp_trace_c);
161 fprintf(swi_fp_trace_c, "#if 0\n");
162 fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd_ut);
163 fprintf(swi_fp_trace_c, " geolon = %.9f;", geolon);
164 fprintf(swi_fp_trace_c, " geolat = %.9f;", geolat);
165 fprintf(swi_fp_trace_c, " hsys = %d;\n", hsys);
166 fprintf(swi_fp_trace_c, " retc = swe_houses(tjd, geolat, geolon, hsys, cusp, ascmc);\n");
167 fprintf(swi_fp_trace_c, " /* swe_houses calls swe_houses_armc as follows: */\n");
168 fprintf(swi_fp_trace_c, "#endif\n");
169 fflush(swi_fp_trace_c);
170 }
171 }
172 #endif
173 retc = swe_houses_armc_ex2(armc, geolat, eps + nutlo[1], hsys, cusp, ascmc, NULL, NULL, NULL);
174 return retc;
175 }
176
177 // For explanation see function swe_houses_ex2() below.
swe_houses_ex(double tjd_ut,int32 iflag,double geolat,double geolon,int hsys,double * cusp,double * ascmc)178 int CALL_CONV swe_houses_ex(double tjd_ut,
179 int32 iflag,
180 double geolat,
181 double geolon,
182 int hsys,
183 double *cusp,
184 double *ascmc)
185 {
186 return swe_houses_ex2(tjd_ut, iflag, geolat, geolon, hsys, cusp, ascmc, NULL, NULL, NULL);
187 }
188
189 /*
190 * Function returns OK or ERR.
191 * cusps are returned in double cusp[13],
192 * or cusp[37] with house system 'G'.
193 * cusp[1...12] houses 1 - 12
194 * ascmc[0...10] additional points:
195 * ascmc[0] = ascendant
196 * ascmc[1] = mc
197 * ascmc[2] = armc
198 * ascmc[3] = vertex
199 * ascmc[4] = equasc * "equatorial ascendant" *
200 * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) *
201 * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) *
202 * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) *
203 * cusp_speed[1...12] speeds (daily motions) of the cusps.
204 * ascmc_speed[0...10] speeds (daily motions) of the additional points.
205 * serr error message or warning
206 */
swe_houses_ex2(double tjd_ut,int32 iflag,double geolat,double geolon,int hsys,double * cusp,double * ascmc,double * cusp_speed,double * ascmc_speed,char * serr)207 int CALL_CONV swe_houses_ex2(double tjd_ut,
208 int32 iflag,
209 double geolat,
210 double geolon,
211 int hsys,
212 double *cusp,
213 double *ascmc,
214 double *cusp_speed,
215 double *ascmc_speed,
216 char *serr)
217 {
218 int i, retc = 0;
219 double armc, eps_mean, nutlo[2];
220 double tjde = tjd_ut + swe_deltat_ex(tjd_ut, iflag, NULL);
221 struct sid_data *sip = &swed.sidd;
222 double xp[6];
223 int retc_makr = 0;
224 int ito;
225 if (toupper(hsys) == 'G')
226 ito = 36;
227 else
228 ito = 12;
229 if ((iflag & SEFLG_SIDEREAL) && !swed.ayana_is_set)
230 swe_set_sid_mode(SE_SIDM_FAGAN_BRADLEY, 0, 0);
231 eps_mean = swi_epsiln(tjde, 0) * RADTODEG;
232 swi_nutation(tjde, 0, nutlo);
233 for (i = 0; i < 2; i++)
234 nutlo[i] *= RADTODEG;
235 if (iflag & SEFLG_NONUT) {
236 for (i = 0; i < 2; i++)
237 nutlo[i] = 0;
238 }
239 #ifdef TRACE
240 swi_open_trace(NULL);
241 if (swi_trace_count <= TRACE_COUNT_MAX) {
242 if (swi_fp_trace_c != NULL) {
243 fputs("\n/*SWE_HOUSES_EX*/\n", swi_fp_trace_c);
244 fprintf(swi_fp_trace_c, "#if 0\n");
245 fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd_ut);
246 fprintf(swi_fp_trace_c, " iflag = %d;\n", iflag);
247 fprintf(swi_fp_trace_c, " geolon = %.9f;", geolon);
248 fprintf(swi_fp_trace_c, " geolat = %.9f;", geolat);
249 fprintf(swi_fp_trace_c, " hsys = %d;\n", hsys);
250 fprintf(swi_fp_trace_c, " retc = swe_houses_ex(tjd, iflag, geolat, geolon, hsys, cusp, ascmc);\n");
251 fprintf(swi_fp_trace_c, " /* swe_houses calls swe_houses_armc as follows: */\n");
252 fprintf(swi_fp_trace_c, "#endif\n");
253 fflush(swi_fp_trace_c);
254 }
255 }
256 #endif
257 /*houses_to_sidereal(tjde, geolat, hsys, eps, cusp, ascmc, iflag);*/
258 armc = swe_degnorm(swe_sidtime0(tjd_ut, eps_mean + nutlo[1], nutlo[0]) * 15 + geolon);
259 //fprintf(stderr, "armc=%f, iflag=%d\n", armc, iflag);
260 if (toupper(hsys) == 'I') { // compute sun declination for sunshine houses
261 int flags = SEFLG_SPEED| SEFLG_EQUATORIAL;
262 retc_makr = swe_calc_ut(tjd_ut, SE_SUN, flags, xp, NULL);
263 if (retc_makr < 0) {
264 // in case of failure, provide Porphyry houses
265 hsys = (int) 'O';
266 }
267 ascmc[9] = xp[1]; // declination in ascmc[9];
268 }
269 if (iflag & SEFLG_SIDEREAL) {
270 if (sip->sid_mode & SE_SIDBIT_ECL_T0)
271 retc = sidereal_houses_ecl_t0(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);
272 else if (sip->sid_mode & SE_SIDBIT_SSY_PLANE)
273 retc = sidereal_houses_ssypl(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);
274 else
275 retc = sidereal_houses_trad(tjde, iflag, armc, eps_mean + nutlo[1], nutlo[0], geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);
276 } else {
277 retc = swe_houses_armc_ex2(armc, geolat, eps_mean + nutlo[1], hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);
278 if (toupper(hsys) == 'I')
279 ascmc[9] = xp[1]; // declination in ascmc[9];
280 }
281 if (iflag & SEFLG_RADIANS) {
282 for (i = 1; i <= ito; i++)
283 cusp[i] *= DEGTORAD;
284 for (i = 0; i < SE_NASCMC; i++)
285 ascmc[i] *= DEGTORAD;
286 }
287 if (retc_makr < 0)
288 return retc_makr;
289 return retc;
290 }
291
292 /*
293 * houses to sidereal
294 * ------------------
295 * there are two methods:
296 * a) the traditional one
297 * houses are computed tropically, then nutation and the ayanamsa
298 * are subtracted.
299 * b) the projection on the ecliptic of t0
300 * The house computation is then as follows:
301 *
302 * Be t the birth date and t0 the epoch at which ayanamsa = 0.
303 * 1. Compute the angle between the mean ecliptic at t0 and
304 * the true equator at t.
305 * The intersection point of these two circles we call the
306 * "auxiliary vernal point", and the angle between them the
307 * "auxiliary obliquity".
308 * 2. Compute the distance of the auxiliary vernal point from the
309 * vernal point at t. (this is a section on the equator)
310 * 3. subtract this value from the armc of t = aux. armc.
311 * 4. Compute the axes and houses for this aux. armc and aux. obliquity.
312 * 5. Compute the distance between the auxiliary vernal point and the
313 * vernal point at t0 (this is the ayanamsa at t, measured on the
314 * ecliptic of t0)
315 * 6. subtract this distance from all house cusps.
316 * 7. subtract ayanamsa_t0 from all house cusps.
317 */
sidereal_houses_ecl_t0(double tjde,double armc,double eps,double * nutlo,double lat,int hsys,double * cusp,double * ascmc,double * cusp_speed,double * ascmc_speed,char * serr)318 static int sidereal_houses_ecl_t0(double tjde,
319 double armc,
320 double eps,
321 double *nutlo,
322 double lat,
323 int hsys,
324 double *cusp,
325 double *ascmc,
326 double *cusp_speed,
327 double *ascmc_speed,
328 char *serr)
329 {
330 int i, j, retc = OK;
331 double x[6], xvpx[6], x2[6], epst0, xnorm[6];
332 double rxy, rxyz, c2, epsx, sgn, fac, dvpx, dvpxe;
333 double armcx;
334 struct sid_data *sip = &swed.sidd;
335 int ito;
336 if (toupper(hsys) == 'G')
337 ito = 36;
338 else
339 ito = 12;
340 /* epsilon at t0 */
341 epst0 = swi_epsiln(sip->t0, 0);
342 /* cartesian coordinates of an imaginary moving body on the
343 * the mean ecliptic of t0; we take the vernal point: */
344 x[0] = x[4] = 1;
345 x[1] = x[2] = x[3] = x[5] = 0;
346 /* to equator */
347 swi_coortrf(x, x, -epst0);
348 swi_coortrf(x+3, x+3, -epst0);
349 /* to tjd_et */
350 swi_precess(x, sip->t0, 0, J_TO_J2000);
351 swi_precess(x, tjde, 0, J2000_TO_J);
352 swi_precess(x+3, sip->t0, 0, J_TO_J2000);
353 swi_precess(x+3, tjde, 0, J2000_TO_J);
354 /* to true equator of tjd_et */
355 swi_coortrf(x, x, (eps - nutlo[1]) * DEGTORAD);
356 swi_coortrf(x+3, x+3, (eps - nutlo[1]) * DEGTORAD);
357 swi_cartpol_sp(x, x);
358 x[0] += nutlo[0] * DEGTORAD;
359 swi_polcart_sp(x, x);
360 swi_coortrf(x, x, -eps * DEGTORAD);
361 swi_coortrf(x+3, x+3, -eps * DEGTORAD);
362 /* now, we have the moving point precessed to tjd_et.
363 * next, we compute the auxiliary epsilon: */
364 swi_cross_prod(x, x+3, xnorm);
365 rxy = xnorm[0] * xnorm[0] + xnorm[1] * xnorm[1];
366 c2 = (rxy + xnorm[2] * xnorm[2]);
367 rxyz = sqrt(c2);
368 rxy = sqrt(rxy);
369 epsx = asin(rxy / rxyz) * RADTODEG; /* 1a */
370 /* auxiliary vernal point */
371 if (fabs(x[5]) < 1e-15)
372 x[5] = 1e-15;
373 fac = x[2] / x[5];
374 sgn = x[5] / fabs(x[5]);
375 for (j = 0; j <= 2; j++)
376 xvpx[j] = (x[j] - fac * x[j+3]) * sgn; /* 1b */
377 /* distance of the auxiliary vernal point from
378 * the zero point at tjd_et (a section on the equator): */
379 swi_cartpol(xvpx, x2);
380 dvpx = x2[0] * RADTODEG; /* 2 */
381 /* auxiliary armc */
382 armcx = swe_degnorm(armc - dvpx); /* 3 */
383 /* compute axes and houses: */
384 retc = swe_houses_armc_ex2(armcx, lat, epsx, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); /* 4 */
385 /* distance between auxiliary vernal point and
386 * vernal point of t0 (a section on the sidereal plane) */
387 dvpxe = acos(swi_dot_prod_unit(x, xvpx)) * RADTODEG; /* 5 */
388 if (tjde < sip->t0)
389 dvpxe = -dvpxe;
390 for (i = 1; i <= ito; i++) /* 6, 7 */
391 cusp[i] = swe_degnorm(cusp[i] - dvpxe - sip->ayan_t0);
392 for (i = 0; i <= SE_NASCMC; i++) {
393 if (i == 2) /* armc */
394 continue;
395 ascmc[i] = swe_degnorm(ascmc[i] - dvpxe - sip->ayan_t0);
396 }
397 if (hsys == 'N') { /* 1 = 0° Aries */
398 for (i = 1; i <= ito; i++) {
399 cusp[i] = (i - 1) * 30;
400 }
401 }
402 return retc;
403 }
404
405 /*
406 * Be t the birth date and t0 the epoch at which ayanamsa = 0.
407 * 1. Compute the angle between the solar system rotation plane and
408 * the true equator at t.
409 * The intersection point of these two circles we call the
410 * "auxiliary vernal point", and the angle between them the
411 * "auxiliary obliquity".
412 * 2. Compute the distance of the auxiliary vernal point from the
413 * zero point at t. (this is a section on the equator)
414 * 3. subtract this value from the armc of t = aux. armc.
415 * 4. Compute the axes and houses for this aux. armc and aux. obliquity.
416 * 5. Compute the distance between the auxiliary vernal point at t
417 * and the zero point of the solar system plane J2000
418 * (a section measured on the solar system plane)
419 * 6. subtract this distance from all house cusps.
420 * 7. compute the ayanamsa of J2000 on the solar system plane,
421 * referred to t0
422 * 8. subtract ayanamsa_t0 from all house cusps.
423 * 9. subtract ayanamsa_2000 from all house cusps.
424 */
sidereal_houses_ssypl(double tjde,double armc,double eps,double * nutlo,double lat,int hsys,double * cusp,double * ascmc,double * cusp_speed,double * ascmc_speed,char * serr)425 static int sidereal_houses_ssypl(double tjde,
426 double armc,
427 double eps,
428 double *nutlo,
429 double lat,
430 int hsys,
431 double *cusp,
432 double *ascmc,
433 double *cusp_speed,
434 double *ascmc_speed,
435 char *serr)
436 {
437 int i, j, retc = OK;
438 double x[6], x0[6], xvpx[6], x2[6], xnorm[6];
439 double rxy, rxyz, c2, epsx, eps2000, sgn, fac, dvpx, dvpxe, x00;
440 double armcx;
441 struct sid_data *sip = &swed.sidd;
442 int ito;
443 if (toupper(hsys) == 'G')
444 ito = 36;
445 else
446 ito = 12;
447 eps2000 = swi_epsiln(J2000, 0);
448 /* cartesian coordinates of the zero point on the
449 * the solar system rotation plane */
450 x[0] = x[4] = 1;
451 x[1] = x[2] = x[3] = x[5] = 0;
452 /* to ecliptic 2000 */
453 swi_coortrf(x, x, -SSY_PLANE_INCL);
454 swi_coortrf(x+3, x+3, -SSY_PLANE_INCL);
455 swi_cartpol_sp(x, x);
456 x[0] += SSY_PLANE_NODE_E2000;
457 swi_polcart_sp(x, x);
458 /* to equator 2000 */
459 swi_coortrf(x, x, -eps2000);
460 swi_coortrf(x+3, x+3, -eps2000);
461 /* to mean equator of t */
462 swi_precess(x, tjde, 0, J2000_TO_J);
463 swi_precess(x+3, tjde, 0, J2000_TO_J);
464 /* to true equator of t */
465 swi_coortrf(x, x, (eps - nutlo[1]) * DEGTORAD);
466 swi_coortrf(x+3, x+3, (eps - nutlo[1]) * DEGTORAD);
467 swi_cartpol_sp(x, x);
468 x[0] += nutlo[0] * DEGTORAD;
469 swi_polcart_sp(x, x);
470 swi_coortrf(x, x, -eps * DEGTORAD);
471 swi_coortrf(x+3, x+3, -eps * DEGTORAD);
472 /* now, we have the moving point precessed to tjd_et.
473 * next, we compute the auxiliary epsilon: */
474 swi_cross_prod(x, x+3, xnorm);
475 rxy = xnorm[0] * xnorm[0] + xnorm[1] * xnorm[1];
476 c2 = (rxy + xnorm[2] * xnorm[2]);
477 rxyz = sqrt(c2);
478 rxy = sqrt(rxy);
479 epsx = asin(rxy / rxyz) * RADTODEG; /* 1a */
480 /* auxiliary vernal point */
481 if (fabs(x[5]) < 1e-15)
482 x[5] = 1e-15;
483 fac = x[2] / x[5];
484 sgn = x[5] / fabs(x[5]);
485 for (j = 0; j <= 2; j++)
486 xvpx[j] = (x[j] - fac * x[j+3]) * sgn; /* 1b */
487 /* distance of the auxiliary vernal point from
488 * mean vernal point at tjd_et (a section on the equator): */
489 swi_cartpol(xvpx, x2);
490 dvpx = x2[0] * RADTODEG; /* 2 */
491 /* auxiliary armc */
492 armcx = swe_degnorm(armc - dvpx); /* 3 */
493 /* compute axes and houses: */
494 retc = swe_houses_armc_ex2(armcx, lat, epsx, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); /* 4 */
495 /* distance between the auxiliary vernal point at t and
496 * the sidereal zero point of 2000 at t
497 * (a section on the sidereal plane).
498 */
499 dvpxe = acos(swi_dot_prod_unit(x, xvpx)) * RADTODEG; /* 5 */
500 /* (always positive for dates after 5400 bc) */
501 dvpxe -= SSY_PLANE_NODE * RADTODEG;
502 /* ayanamsa between t0 and J2000, measured on solar system plane: */
503 /* position of zero point of t0 */
504 x0[0] = 1;
505 x0[1] = x0[2] = 0;
506 /* zero point of t0 in J2000 system */
507 if (sip->t0 != J2000)
508 swi_precess(x0, sip->t0, 0, J_TO_J2000);
509 /* zero point to ecliptic 2000 */
510 swi_coortrf(x0, x0, eps2000);
511 /* to solar system plane */
512 swi_cartpol(x0, x0);
513 x0[0] -= SSY_PLANE_NODE_E2000;
514 swi_polcart(x0, x0);
515 swi_coortrf(x0, x0, SSY_PLANE_INCL);
516 swi_cartpol(x0, x0);
517 x0[0] += SSY_PLANE_NODE;
518 x00 = x0[0] * RADTODEG; /* 7 */
519 for (i = 1; i <= ito; i++) /* 6, 8, 9 */
520 cusp[i] = swe_degnorm(cusp[i] - dvpxe - sip->ayan_t0 - x00);
521 for (i = 0; i <= SE_NASCMC; i++) {
522 if (i == 2) /* armc */
523 continue;
524 ascmc[i] = swe_degnorm(ascmc[i] - dvpxe - sip->ayan_t0 - x00);
525 }
526 if (hsys == 'N') { /* 1 = 0° Aries */
527 for (i = 1; i <= ito; i++) {
528 cusp[i] = (i - 1) * 30;
529 }
530 }
531 return retc;
532 }
533
534 /* common simplified procedure */
sidereal_houses_trad(double tjde,int32 iflag,double armc,double eps,double nutl,double lat,int hsys,double * cusp,double * ascmc,double * cusp_speed,double * ascmc_speed,char * serr)535 static int sidereal_houses_trad(double tjde,
536 int32 iflag,
537 double armc,
538 double eps,
539 double nutl,
540 double lat,
541 int hsys,
542 double *cusp,
543 double *ascmc,
544 double *cusp_speed,
545 double *ascmc_speed,
546 char *serr)
547 {
548 int i, retc = OK;
549 double ay;
550 int ito;
551 int ihs = toupper(hsys);
552 int ihs2 = ihs;
553 // ay = swe_get_ayanamsa(tjde);
554 //fprintf(stderr, "ay=%f\n", ay);
555 retc = swe_get_ayanamsa_ex(tjde, iflag, &ay, NULL);
556 //fprintf(stderr, "ay=%f\n", ay);
557 //fprintf(stderr, "nutl=%f\n", nutl);
558 if (ihs == 'G')
559 ito = 36;
560 else
561 ito = 12;
562 if (ihs == 'W') /* whole sign houses: treat as 'E' and fix later */
563 ihs2 = 'E';
564 //fprintf(stderr, "armc=%f\n", armc);
565 //if (hsys == 'P') fprintf(stderr, "ay=%f, t=%f %c", ay, tjde, (char) hsys);
566 retc = swe_houses_armc_ex2(armc, lat, eps, ihs2, cusp, ascmc, cusp_speed, ascmc_speed, serr);
567 //if (hsys == 'P') fprintf(stderr, " h1=%f", cusp[1]);
568 for (i = 1; i <= ito; i++) {
569 //cusp[i] = swe_degnorm(cusp[i] - ay - nutl);
570 cusp[i] = swe_degnorm(cusp[i] - ay);
571 if (ihs == 'W') /* whole sign houses */
572 cusp[i] -= fmod(cusp[i], 30);
573 }
574 if (ihs == 'N') { /* 1 = 0° Aries */
575 for (i = 1; i <= ito; i++) {
576 cusp[i] = (i - 1) * 30;
577 }
578 }
579 for (i = 0; i < SE_NASCMC; i++) {
580 if (i == 2) /* armc */
581 continue;
582 //ascmc[i] = swe_degnorm(ascmc[i] - ay - nutl);
583 ascmc[i] = swe_degnorm(ascmc[i] - ay);
584 }
585 //if (hsys == 'P') fprintf(stderr, " => %f\n", cusp[1]);
586 return retc;
587 }
588
589 // For explanation see function swe_houses_armc_ex2() below.
swe_houses_armc(double armc,double geolat,double eps,int hsys,double * cusp,double * ascmc)590 int CALL_CONV swe_houses_armc(
591 double armc,
592 double geolat,
593 double eps,
594 int hsys,
595 double *cusp,
596 double *ascmc)
597 {
598 return swe_houses_armc_ex2(armc, geolat, eps, hsys, cusp, ascmc, NULL, NULL, NULL);
599 }
600
601 /*
602 * Function returns OK or ERR.
603 * this function is required for very special computations
604 * where no date is given for house calculation,
605 * e.g. for composite charts or progressive charts.
606 * cusps are returned in double cusp[13],
607 * or cusp[37] with house system 'G'.
608 * cusp[1...12] houses 1 - 12
609 * ascmc[0...10] additional points:
610 * ascmc[0] = ascendant
611 * ascmc[1] = mc
612 * ascmc[2] = armc
613 * ascmc[3] = vertex
614 * ascmc[4] = equasc * "equatorial ascendant" *
615 * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) *
616 * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) *
617 * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) *
618 * cusp_speed[1...12] speeds (daily motions) of the cusps.
619 * ascmc_speed[0...10] speeds (daily motions) of the additional points.
620 * serr error message or warning
621 */
swe_houses_armc_ex2(double armc,double geolat,double eps,int hsys,double * cusp,double * ascmc,double * cusp_speed,double * ascmc_speed,char * serr)622 int CALL_CONV swe_houses_armc_ex2(
623 double armc,
624 double geolat,
625 double eps,
626 int hsys,
627 double *cusp,
628 double *ascmc,
629 double *cusp_speed,
630 double *ascmc_speed,
631 char *serr)
632 {
633 struct houses h, hm1, hp1;
634 int i, retc = 0, rm1, rp1;
635 int ito;
636 static double saved_sundec = 99;
637 if (toupper(hsys) == 'G')
638 ito = 36;
639 else
640 ito = 12;
641 armc = swe_degnorm(armc);
642 h.do_speed = FALSE;
643 h.do_hspeed = FALSE;
644 if (ascmc_speed != NULL || cusp_speed != NULL)
645 h.do_speed = TRUE; // is needed if cusp_speed wanted
646 if (cusp_speed != NULL)
647 h.do_hspeed = TRUE;
648 if (toupper(hsys) == 'I') { // declination for sunshine houses
649 if (ascmc[9] == 99) {
650 h.sundec = 0;
651 if (saved_sundec != 99) h.sundec = saved_sundec;
652 } else {
653 h.sundec = ascmc[9];
654 saved_sundec = h.sundec;
655 }
656 }
657 retc = CalcH(armc, geolat, eps, (char)hsys, &h);
658 cusp[0] = 0;
659 // on failure, we only have 12 Porphyry cusps
660 if (retc < 0) {
661 ito = 12;
662 if (serr != NULL) strcpy(serr, h.serr);
663 }
664 for (i = 1; i <= ito; i++) {
665 cusp[i] = h.cusp[i];
666 if (h.do_hspeed) cusp_speed[i] = h.cusp_speed[i];
667 }
668 ascmc[0] = h.ac; /* Asc */
669 ascmc[1] = h.mc; /* Mid */
670 ascmc[2] = armc;
671 ascmc[3] = h.vertex;
672 ascmc[4] = h.equasc;
673 ascmc[5] = h.coasc1; /* "co-ascendant" (W. Koch) */
674 ascmc[6] = h.coasc2; /* "co-ascendant" (M. Munkasey) */
675 ascmc[7] = h.polasc; /* "polar ascendant" (M. Munkasey) */
676 for (i = SE_NASCMC; i < 10; i++)
677 ascmc[i] = 0;
678 if (toupper(hsys) == 'I') // declination for sunshine houses
679 ascmc[9] = h.sundec ;
680 if (h.do_speed && ascmc_speed != NULL) {
681 ascmc_speed[0] = h.ac_speed; /* Asc */
682 ascmc_speed[1] = h.mc_speed; /* Mid */
683 ascmc_speed[2] = h.armc_speed;
684 ascmc_speed[3] = h.vertex_speed;
685 ascmc_speed[4] = h.equasc_speed;
686 ascmc_speed[5] = h.coasc1_speed; /* "co-ascendant" (W. Koch) */
687 ascmc_speed[6] = h.coasc2_speed; /* "co-ascendant" (M. Munkasey) */
688 ascmc_speed[7] = h.polasc_speed; /* "polar ascendant" (M. Munkasey) */
689 for (i = SE_NASCMC; i < 10; i++)
690 ascmc_speed[i] = 0;
691 }
692 if (h.do_interpol) { // must compute cusp_speed via interpolation
693 double dt = 1.0 / 86400;
694 double darmc = dt * ARMCS;
695 hm1.do_speed = FALSE;
696 hm1.do_hspeed = FALSE;
697 hp1.do_speed = FALSE;
698 hp1.do_hspeed = FALSE;
699 if (toupper(hsys) == 'I') {
700 hm1.sundec = h.sundec;
701 hp1.sundec = h.sundec;
702 }
703 rm1 = CalcH(armc - darmc, geolat, eps, (char)hsys, &hm1);
704 rp1 = CalcH(armc + darmc, geolat, eps, (char)hsys, &hp1);
705 if (rp1 >= 0 && rm1 >=0) {
706 if (fabs(swe_difdeg2n(hp1.ac, h.ac)) > 90) {
707 hp1 = h; // use only upper interval
708 dt = dt / 2;
709 } else if (fabs(swe_difdeg2n(hm1.ac, h.ac)) > 90) {
710 hm1 = h; // use only lower interval
711 dt = dt / 2;
712 }
713 for (i = 1; i <= 12; i++) {
714 double dx = swe_difdeg2n(hp1.cusp[i], hm1.cusp[i]);
715 cusp_speed[i] = dx / 2 / dt ;
716 }
717 }
718 }
719 #ifdef TRACE
720 swi_open_trace(NULL);
721 if (swi_trace_count <= TRACE_COUNT_MAX) {
722 if (swi_fp_trace_c != NULL) {
723 fputs("\n/*SWE_HOUSES_ARMC_EX2*/\n", swi_fp_trace_c);
724 fprintf(swi_fp_trace_c, " armc = %.9f;", armc);
725 fprintf(swi_fp_trace_c, " geolat = %.9f;", geolat);
726 fprintf(swi_fp_trace_c, " eps = %.9f;", eps);
727 fprintf(swi_fp_trace_c, " hsys = %d;\n", hsys);
728 fprintf(swi_fp_trace_c, " retc = swe_houses_armc_ex2(armc, geolat, eps, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);\n");
729 fputs(" printf(\"swe_houses_armc_ex2: %f\\t%f\\t%f\\t%c\\t\\n\", ", swi_fp_trace_c);
730 fputs(" armc, geolat, eps, hsys);\n", swi_fp_trace_c);
731 fputs(" printf(\"retc = %d\\n\", retc);\n", swi_fp_trace_c);
732 fputs(" printf(\"cusp:\\n\");\n", swi_fp_trace_c);
733 fputs(" for (i = 1; i <= 12; i++)\n", swi_fp_trace_c);
734 fputs(" printf(\" %d\\t%f\\n\", i, cusp[i]);\n", swi_fp_trace_c);
735 fputs(" printf(\"ascmc:\\n\");\n", swi_fp_trace_c);
736 fputs(" for (i = 0; i < 10; i++)\n", swi_fp_trace_c);
737 fputs(" printf(\" %d\\t%f\\n\", i, ascmc[i]);\n", swi_fp_trace_c);
738 fputs(" printf(\"cusp_speed:\\n\");\n", swi_fp_trace_c);
739 fputs(" for (i = 1; i <= 12; i++)\n", swi_fp_trace_c);
740 fputs(" printf(\" %d\\t%f\\n\", i, cusp_speed[i]);\n", swi_fp_trace_c);
741 fputs(" printf(\"ascmc_speed:\\n\");\n", swi_fp_trace_c);
742 fputs(" for (i = 0; i < 10; i++)\n", swi_fp_trace_c);
743 fputs(" printf(\" %d\\t%f\\n\", i, ascmc_speed[i]);\n", swi_fp_trace_c);
744 fflush(swi_fp_trace_c);
745 }
746 if (swi_fp_trace_out != NULL) {
747 fprintf(swi_fp_trace_out, "swe_houses_armc_ex2: %f\t%f\t%f\t%c\t\n", armc, geolat, eps, hsys);
748 fprintf(swi_fp_trace_out, "retc = %d\n", retc);
749 fputs("cusp:\n", swi_fp_trace_out);
750 for (i = 1; i <= 12; i++)
751 fprintf(swi_fp_trace_out, " %d\t%f\n", i, cusp[i]);
752 fputs("ascmc:\n", swi_fp_trace_out);
753 for (i = 0; i < 10; i++)
754 fprintf(swi_fp_trace_out, " %d\t%f\n", i, ascmc[i]);
755 fflush(swi_fp_trace_out);
756 }
757 }
758 #endif
759 #if 0
760 /* for test of swe_house_pos().
761 * 1st house will be 0, second 30, etc. */
762 for (i = 1; i <=12; i++) {
763 double x[6];
764 x[0] = cusp[i]; x[1] = 0; x[2] = 1;
765 cusp[i] = (swe_house_pos(armc, geolat, eps, hsys, x, NULL) - 1) * 30;
766 }
767 #endif
768 return retc;
769 }
770
771 /* for APC houses */
772 /* n number of house
773 * ph geographic latitude
774 * e ecliptic obliquity
775 * az armc
776 */
apc_sector(int n,double ph,double e,double az)777 static double apc_sector(int n, double ph, double e, double az)
778 {
779 int k, is_below_hor = 0;
780 double kv, a, dasc, dret;
781 /* kv: ascensional difference of the ascendant */
782 /* dasc: declination of the ascendant */
783 if (fabs(ph * RADTODEG) > 90 - VERY_SMALL) {
784 kv = 0;
785 dasc = 0;
786 } else {
787 kv = atan(tan(ph) * tan(e) * cos(az)/(1 + tan(ph) * tan(e) * sin(az)));
788 if (fabs(ph * RADTODEG) < VERY_SMALL) {
789 dasc = (90 - VERY_SMALL) * DEGTORAD;
790 if (ph < 0)
791 dasc = -dasc;
792 } else {
793 dasc = atan(sin(kv) / tan(ph));
794 }
795 }
796 /* note, at polar circles, when the mc sinks below the horizon,
797 * kv and dasc change sign in the above formulae.
798 * this is what we need, because the ascendand jumps by 180 deg */
799 /* printf("%f, %f\n", kv*RADTODEG, dasc*RADTODEG); */
800 if (n < 8) {
801 is_below_hor = 1; /* 1 and 7 are included here */
802 k = n - 1;
803 } else {
804 k = n - 13;
805 }
806 /* az + PI/2 + kv = armc + 90 + asc. diff. = right ascension of ascendant
807 * PI/2 +- kv = semi-diurnal or seminocturnal arc of ascendant
808 * a = right ascension of house cusp on apc circle (ascendant-parallel
809 * circle), with declination dasc */
810 if (is_below_hor) {
811 a = kv + az + PI/2 + k * (PI/2 - kv) / 3;
812 } else {
813 a = kv + az + PI/2 + k * (PI/2 + kv) / 3;
814 }
815 a = swe_radnorm(a);
816 dret = atan2(tan(dasc) * tan(ph) * sin(az) + sin(a),
817 cos(e) * (tan(dasc) * tan(ph) * cos(az) + cos(a)) + sin(e) * tan(ph) * sin(az - a));
818 dret = swe_degnorm(dret * RADTODEG);
819 return dret;
820 }
821
swe_house_name(int hsys)822 char *CALL_CONV swe_house_name(int hsys)
823 {
824 int h = hsys;
825 if (h != 'i') h = toupper(h);
826 switch (h) {
827 case 'A': return "equal";
828 case 'B': return "Alcabitius";
829 case 'C': return "Campanus";
830 case 'D': return "equal (MC)";
831 case 'E': return "equal";
832 case 'F': return "Carter poli-equ.";
833 case 'G': return "Gauquelin sectors";
834 case 'H': return "horizon/azimut";
835 case 'I': return "Sunshine";
836 case 'i': return "Sunshine/alt.";
837 case 'K': return "Koch";
838 case 'L': return "Pullen SD";
839 case 'M': return "Morinus";
840 case 'N': return "equal/1=Aries";
841 case 'O': return "Porphyry";
842 case 'Q': return "Pullen SR";
843 case 'R': return "Regiomontanus";
844 case 'S': return "Sripati";
845 case 'T': return "Polich/Page";
846 case 'U': return "Krusinski-Pisa-Goelzer";
847 case 'V': return "equal/Vehlow";
848 case 'W': return "equal/ whole sign";
849 case 'X': return "axial rotation system/Meridian houses";
850 case 'Y': return "APC houses";
851 default: return "Placidus";
852 }
853 }
854
855 // How to deal with Sunshine houses if the southern crossing point of Equator
856 // and Ecliptic is under the horizon:
857 // We follow the proposal by Dieter Koch, who wants to keep it in alalogy with
858 // Regiomontanus, where we keep the MC above the horozon, by switching it to the noth.
859 // This results in an clockwise sequence of house cusps in the chart.
860 //
861 // One can argue that the MC should be kept south, even when it is under the horizon.
862 // This would keep the sequence of houses in the chart counterclockwise as usual.
863 // To achieve it, the offsets on the diurnal arcs must be inverted.
864 #define SUNSHINE_KEEP_MC_SOUTH 0 // must be 0 or 1
865
swi_armc_to_mc(double armc,double eps)866 double swi_armc_to_mc(double armc, double eps)
867 {
868 double tant, mc;
869 if (fabs(armc - 90) > VERY_SMALL
870 && fabs(armc - 270) > VERY_SMALL) {
871 tant = tand(armc);
872 mc = atand(tant / cosd(eps));
873 if (armc > 90 && armc <= 270)
874 mc = swe_degnorm(mc + 180);
875 } else {
876 if (fabs(armc - 90) <= VERY_SMALL)
877 mc = 90;
878 else
879 mc = 270;
880 } /* if */
881 return mc;
882 }
883
884 //#define DEBUG_PLAC_ITER 1
885 #define VERY_SMALL_PLAC_ITER (1.0 / 360000.0 )
CalcH(double th,double fi,double ekl,char hsy,struct houses * hsp)886 static int CalcH(
887 double th, double fi, double ekl, char hsy, struct houses *hsp)
888 /* *********************************************************
889 * Arguments: th = sidereal time (angle 0..360 degrees
890 * hsy = letter code for house system;
891 * A equal
892 * E equal
893 * B Alcabitius
894 * C Campanus
895 * D equal (MC)
896 * F Carter "Poli-Equatorial"
897 * G 36 Gauquelin sectors
898 * H horizon / azimut
899 * I Sunshine solution Treindl
900 * i Sunshine solution Makransky
901 * K Koch
902 * L Pullen SD "sinusoidal delta", ex Neo-Porphyry
903 * M Morinus
904 * N equal/1=Aries
905 * O Porphyry
906 * P Placidus
907 * Q Pullen SR "sinusoidal ratio"
908 * R Regiomontanus
909 * S Sripati
910 * T Polich/Page ("topocentric")
911 * U Krusinski-Pisa-Goelzer
912 * V equal Vehlow
913 * W equal, whole sign
914 * X axial rotation system/ Meridian houses
915 * Y APC houses
916 * fi = geographic latitude
917 * ekl = obliquity of the ecliptic
918 * *********************************************************
919 * Koch and Placidus don't work in the polar circle.
920 * We swap MC/IC so that MC is always before AC in the zodiac
921 * We than divide the quadrants into 3 equal parts.
922 * *********************************************************
923 * All angles are expressed in degrees.
924 * Special trigonometric functions sind, cosd etc. are
925 * implemented for arguments in degrees.
926 ***********************************************************/
927 {
928 double tane, tanfi, cosfi, tant, sina, cosa, th2;
929 double a, c, f, fh1, fh2, xh1, xh2, rectasc, ad3, acmc, vemc;
930 int i, ih, ih2, retc = OK;
931 double sine, cose;
932 double x[3], krHorizonLon; /* BK 14.02.2006 */
933 int niter_max = 100; // maximum iterations allowed with Placidus
934 double cuspsv;
935 *hsp->serr = '\0';
936 hsp->do_interpol = 0;
937 cose = cosd(ekl);
938 sine = sind(ekl);
939 tane = tand(ekl);
940 /* north and south poles */
941 if (fabs(fabs(fi) - 90) < VERY_SMALL) {
942 if (fi < 0)
943 fi = -90 + VERY_SMALL;
944 else
945 fi = 90 - VERY_SMALL;
946 }
947 tanfi = tand(fi);
948 /* mc */
949 if (fabs(th - 90) > VERY_SMALL
950 && fabs(th - 270) > VERY_SMALL) {
951 tant = tand(th);
952 hsp->mc = atand(tant / cose);
953 if (th > 90 && th <= 270)
954 hsp->mc = swe_degnorm(hsp->mc + 180);
955 } else {
956 if (fabs(th - 90) <= VERY_SMALL)
957 hsp->mc = 90;
958 else
959 hsp->mc = 270;
960 } /* if */
961 hsp->mc = swe_degnorm(hsp->mc);
962 if (hsp->do_speed) hsp->mc_speed = AscDash(th, 0, sine, cose);
963 /* ascendant */
964 hsp->ac = Asc1(th + 90, fi, sine, cose);
965 if (hsp->do_speed)
966 hsp->ac_speed = AscDash(th + 90, fi, sine, cose);
967 if (hsp->do_hspeed) {
968 for (i = 0; i <= 12; i++)
969 hsp->cusp_speed[i] = 0;
970 }
971 hsp->armc_speed = ARMCS;
972 // these cusp[1] and cusp[10] values may be changed further down for some house systems
973 hsp->cusp[1] = hsp->ac;
974 hsp->cusp[10] = hsp->mc;
975 if (hsp->do_hspeed) {
976 hsp->cusp_speed[1] = hsp->ac_speed;
977 hsp->cusp_speed[10] = hsp->mc_speed;
978 }
979 /* we respect smaller case letter for i, otherwise they are deprecated */
980 if (hsy > 95 && hsy != 'i') {
981 sprintf(hsp->serr, "use of lower case letters like %c for house systems is deprecated", hsy);
982 hsy = (char) (hsy - 32);/* translate into capital letter */
983 }
984 switch (hsy) {
985 case 'A': /* equal houses */
986 case 'E':
987 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
988 if (acmc < 0) {
989 /* within polar circle we swap AC/DC if AC is on wrong side */
990 hsp->ac = swe_degnorm(hsp->ac + 180);
991 hsp->cusp[1] = hsp->ac;
992 }
993 for (i = 2; i <=12; i++) {
994 hsp->cusp[i] = swe_degnorm(hsp->cusp[1] + (i-1) * 30);
995 }
996 if (hsp->do_hspeed) {
997 for (i = 1; i <=12; i++) {
998 hsp->cusp_speed[i] = hsp->ac_speed;
999 }
1000 }
1001 break;
1002 case 'D': /* equal, begin at MC */
1003 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1004 if (acmc < 0) {
1005 /* within polar circle we swap AC/DC if AC is on wrong side */
1006 hsp->ac = swe_degnorm(hsp->ac + 180);
1007 }
1008 hsp->cusp[10] = hsp->mc;
1009 for (i = 11; i <= 12; i++)
1010 hsp->cusp[i] = swe_degnorm(hsp->cusp[10] + (i-10) * 30);
1011 for (i = 1; i <= 9; i++)
1012 hsp->cusp[i] = swe_degnorm(hsp->cusp[10] + (i + 2) * 30);
1013 if (hsp->do_hspeed) {
1014 for (i = 1; i <=12; i++) {
1015 hsp->cusp_speed[i] = hsp->mc_speed;
1016 }
1017 }
1018 break;
1019 case 'C': /* Campanus houses and Horizon or Azimut system */
1020 case 'H':
1021 if (hsy == 'H') {
1022 if (fi > 0)
1023 fi = 90 - fi;
1024 else
1025 fi = -90 - fi;
1026 /* equator */
1027 if (fabs(fabs(fi) - 90) < VERY_SMALL) {
1028 if (fi < 0)
1029 fi = -90 + VERY_SMALL;
1030 else
1031 fi = 90 - VERY_SMALL;
1032 }
1033 th = swe_degnorm(th + 180);
1034 }
1035 fh1 = asind(sind (fi) / 2);
1036 fh2 = asind(sqrt (3.0) / 2 * sind(fi));
1037 cosfi = cosd(fi);
1038 if (fabs(cosfi) == 0) { /* '==' should be save! */
1039 if (fi > 0)
1040 xh1 = xh2 = 90; /* cosfi = VERY_SMALL; */
1041 else
1042 xh1 = xh2 = 270; /* cosfi = -VERY_SMALL; */
1043 } else {
1044 xh1 = atand(sqrt (3.0) / cosfi);
1045 xh2 = atand(1 / sqrt (3.0) / cosfi);
1046 }
1047 hsp->cusp[11] = Asc1(th + 90 - xh1, fh1, sine, cose);
1048 hsp->cusp[12] = Asc1(th + 90 - xh2, fh2, sine, cose);
1049 if (hsy == 'H')
1050 hsp->cusp[1] = Asc1(th + 90, fi, sine, cose);
1051 hsp->cusp[2] = Asc1(th + 90 + xh2, fh2, sine, cose);
1052 hsp->cusp[3] = Asc1(th + 90 + xh1, fh1, sine, cose);
1053 if (hsp->do_hspeed) {
1054 hsp->cusp_speed[11] = AscDash(th + 90 - xh1, fh1, sine, cose);
1055 hsp->cusp_speed[12] = AscDash(th + 90 - xh2, fh2, sine, cose);
1056 if (hsy == 'H')
1057 hsp->cusp_speed[1] = AscDash(th + 90, fi, sine, cose);
1058 hsp->cusp_speed[2] = AscDash(th + 90 + xh2, fh2, sine, cose);
1059 hsp->cusp_speed[3] = AscDash(th + 90 + xh1, fh1, sine, cose);
1060 }
1061 /* within polar circle, when mc sinks below horizon and
1062 * ascendant changes to western hemisphere, all cusps
1063 * must be added 180 degrees.
1064 * houses will be in clockwise direction */
1065 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1066 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1067 if (acmc < 0) {
1068 hsp->ac = swe_degnorm(hsp->ac + 180);
1069 hsp->mc = swe_degnorm(hsp->mc + 180);
1070 for (i = 1; i <= 12; i++) {
1071 if (i >= 4 && i < 10) continue;
1072 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1073 }
1074 }
1075 }
1076 if (hsy == 'H') {
1077 for (i = 1; i <= 3; i++)
1078 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1079 for (i = 11; i <= 12; i++)
1080 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1081 /* restore fi and th */
1082 if (fi > 0)
1083 fi = 90 - fi;
1084 else
1085 fi = -90 - fi;
1086 th = swe_degnorm(th + 180);
1087 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1088 if (acmc < 0) {
1089 hsp->ac = swe_degnorm(hsp->ac + 180);
1090 }
1091 }
1092 break;
1093 case 'I': /* Sunshine houses, solution Treindl */
1094 case 'i': /* Sunshine houses, solution Makranski */
1095 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1096 if (acmc < 0) {
1097 /* we shift axes */
1098 hsp->ac = swe_degnorm(hsp->ac + 180);
1099 hsp->cusp[1] = hsp->ac;
1100 if (! SUNSHINE_KEEP_MC_SOUTH && hsy == 'I') {
1101 hsp->mc = swe_degnorm(hsp->mc + 180);
1102 hsp->cusp[10] = hsp->mc;
1103 }
1104 }
1105 hsp->cusp[4] = swe_degnorm(hsp->cusp[10] + 180);
1106 hsp->cusp[7] = swe_degnorm(hsp->cusp[1] + 180);
1107 if (hsy == 'I') {
1108 retc = sunshine_solution_treindl(th, fi, ekl, hsp);
1109 } else {
1110 retc = sunshine_solution_makransky(th, fi, ekl, hsp);
1111 }
1112 if (retc == ERR) { // only Makransky version does this
1113 strcpy(hsp->serr, "within polar circle, switched to Porphyry");
1114 hsy = 'O';
1115 goto porphyry;
1116 }
1117 hsp->do_interpol = hsp->do_hspeed;
1118 break;
1119 case 'K': /* Koch houses */
1120 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1121 retc = ERR;
1122 strcpy(hsp->serr, "within polar circle, switched to Porphyry");
1123 goto porphyry;
1124 }
1125 sina = sind(hsp->mc) * sine / cosd(fi);
1126 if (sina > 1) sina = 1;
1127 if (sina < -1) sina = -1;
1128 cosa = sqrt(1 - sina * sina); /* always >> 0 */
1129 c = atand(tanfi / cosa);
1130 ad3 = asind(sind(c) * sina) / 3.0;
1131 hsp->cusp[11] = Asc1(th + 30 - 2 * ad3, fi, sine, cose);
1132 hsp->cusp[12] = Asc1(th + 60 - ad3, fi, sine, cose);
1133 hsp->cusp[2] = Asc1(th + 120 + ad3, fi, sine, cose);
1134 hsp->cusp[3] = Asc1(th + 150 + 2 * ad3, fi, sine, cose);
1135 if (hsp->do_hspeed) {
1136 hsp->cusp_speed[11] = AscDash(th + 30 - 2 * ad3, fi, sine, cose);
1137 hsp->cusp_speed[12] = AscDash(th + 60 - ad3, fi, sine, cose);
1138 hsp->cusp_speed[2] = AscDash(th + 120 + ad3, fi, sine, cose);
1139 hsp->cusp_speed[3] = AscDash(th + 150 + 2 * ad3, fi, sine, cose);
1140 }
1141 break;
1142 case 'L': /* Pullen SD sinusoidal delta, ex Neo-Porphyry */
1143 {
1144 double d, q1;
1145 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1146 if (acmc < 0) {
1147 /* within polar circle we swap AC/DC if AC is on wrong side */
1148 hsp->ac = swe_degnorm(hsp->ac + 180);
1149 hsp->cusp[1] = hsp->ac;
1150 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1151 }
1152 q1 = 180 - acmc;
1153 d = (acmc - 90) / 4.0;
1154 if (acmc <= 30) { // is quadrant <= 30, house 11 = zero width.
1155 hsp->cusp[11] = hsp->cusp[12] = swe_degnorm(hsp->mc + acmc / 2);
1156 } else {
1157 hsp->cusp[11] = swe_degnorm(hsp->mc + 30 + d);
1158 hsp->cusp[12] = swe_degnorm(hsp->mc + 60 + 3 * d);
1159 }
1160 d = (q1 - 90) / 4.0;
1161 if (q1 <= 30) { // is quadrant <= 30, house 2 = zero width.
1162 hsp->cusp[2] = hsp->cusp[3] = swe_degnorm(hsp->ac + q1 / 2);
1163 } else {
1164 hsp->cusp[2] = swe_degnorm(hsp->ac + 30 + d);
1165 hsp->cusp[3] = swe_degnorm(hsp->ac + 60 + 3 * d);
1166 }
1167 }
1168 hsp->do_interpol = hsp->do_hspeed;
1169 break;
1170 case 'N': /* whole signs, begin at 0° Aries */
1171 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1172 if (acmc < 0) {
1173 /* within polar circle we swap AC/DC if AC is on wrong side */
1174 hsp->ac = swe_degnorm(hsp->ac + 180);
1175 }
1176 for (i = 1; i <= 12; i++)
1177 hsp->cusp[i] = (i - 1) * 30.0;
1178 break;
1179 case 'O': /* Porphyry houses */
1180 porphyry:
1181 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1182 if (acmc < 0) {
1183 /* within polar circle we swap AC/DC if AC is on wrong side */
1184 hsp->ac = swe_degnorm(hsp->ac + 180);
1185 hsp->cusp[1] = hsp->ac;
1186 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1187 }
1188 hsp->cusp[1] = hsp->ac; // may have been destroyed if defaulting from Gauquelin
1189 hsp->cusp[10] = hsp->mc; // dito
1190 hsp->cusp[2] = swe_degnorm(hsp->ac + (180 - acmc) / 3);
1191 hsp->cusp[3] = swe_degnorm(hsp->ac + (180 - acmc) / 3 * 2);
1192 hsp->cusp[11] = swe_degnorm(hsp->mc + acmc / 3);
1193 hsp->cusp[12] = swe_degnorm(hsp->mc + acmc / 3 * 2);
1194 if (hsp->do_hspeed) {
1195 double q1_speed = hsp->ac_speed - hsp->mc_speed; // rate of growth of quadrant 1
1196 // double q4_speed = hsp->mc_speed - hsp->ac_speed; // rate of growth of quadrant 4
1197 hsp->cusp_speed[1] = hsp->ac_speed; // may have been destroyed if defaulting from Gauquelin
1198 hsp->cusp_speed[10] = hsp->mc_speed; // dito
1199 hsp->cusp_speed[2] = hsp->ac_speed - q1_speed / 3;
1200 hsp->cusp_speed[3] = hsp->ac_speed - q1_speed / 3 * 2;
1201 hsp->cusp_speed[11] = hsp->ac_speed + q1_speed / 3;
1202 hsp->cusp_speed[12] = hsp->ac_speed + q1_speed / 3 * 2;
1203 }
1204 break;
1205 case 'Q': /* Pullen sinusoidal ratio */
1206 {
1207 double q, c, csq, ccr, cqx, two23, third, r, r1, r2, x, xr, xr3, xr4;
1208 third = 1.0 / 3.0;
1209 two23 = pow(2.0 * 2.0, third); // 2^(2/3)
1210 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1211 if (acmc < 0) {
1212 /* within polar circle we swap AC/DC if AC is on wrong side */
1213 hsp->ac = swe_degnorm(hsp->ac + 180);
1214 hsp->cusp[1] = hsp->ac;
1215 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1216 }
1217 q = acmc;
1218 if (q > 90) q = 180 - q;
1219 if (q < 1e-30) { // degenerate case of quadrant = zer0
1220 // r = INFINITY;
1221 x = xr = xr3 = 0;
1222 xr4 = 180;
1223 } else {
1224 c = (180 - q) / q;
1225 csq = c * c;
1226 ccr = pow(csq - c, third); // cuberoot(c^2 -c)
1227 cqx = sqrt(two23 * ccr + 1.0); // sqrt{2^(2/3)*cuberoot(c^2-c) + 1}
1228 r1 = 0.5 * cqx;
1229 r2 = 0.5 * sqrt(-2*(1-2*c) / cqx - two23 * ccr + 2);
1230 r = r1 + r2 - 0.5;
1231 x = q / (2 * r + 1);
1232 xr = r * x;
1233 xr3 = xr * r * r;
1234 xr4 = xr3 * r;
1235 }
1236 if (acmc > 90) {
1237 hsp->cusp[11] = swe_degnorm(hsp->mc + xr3); // house 10 and 12 size xr^3
1238 hsp->cusp[12] = swe_degnorm(hsp->cusp[11] + xr4); // house 11 size xr^4
1239 hsp->cusp[2] = swe_degnorm(hsp->ac + xr); // house 1 and 3 size xr
1240 hsp->cusp[3] = swe_degnorm(hsp->cusp[2] + x); // house 2 size x
1241 } else {
1242 hsp->cusp[11] = swe_degnorm(hsp->mc + xr); // house 10 and 12 size xr
1243 hsp->cusp[12] = swe_degnorm(hsp->cusp[11] + x); // house 11 size x
1244 hsp->cusp[2] = swe_degnorm(hsp->ac + xr3); // house 1 and 3 size xr^3
1245 hsp->cusp[3] = swe_degnorm(hsp->cusp[2] + xr4); // house 2 size xr^4
1246 }
1247 }
1248 hsp->do_interpol = hsp->do_hspeed;
1249 break;
1250 case 'R': /* Regiomontanus houses */
1251 fh1 = atand (tanfi * 0.5);
1252 fh2 = atand (tanfi * cosd(30));
1253 hsp->cusp[11] = Asc1(30 + th, fh1, sine, cose);
1254 hsp->cusp[12] = Asc1(60 + th, fh2, sine, cose);
1255 hsp->cusp[2] = Asc1(120 + th, fh2, sine, cose);
1256 hsp->cusp[3] = Asc1(150 + th, fh1, sine, cose);
1257 if (hsp->do_hspeed) {
1258 hsp->cusp_speed[11] = AscDash(30 + th, fh1, sine, cose);
1259 hsp->cusp_speed[12] = AscDash(60 + th, fh2, sine, cose);
1260 hsp->cusp_speed[2] = AscDash(120 + th, fh2, sine, cose);
1261 hsp->cusp_speed[3] = AscDash(150 + th, fh1, sine, cose);
1262 }
1263 /* within polar circle, when mc sinks below horizon and
1264 * ascendant changes to western hemisphere, all cusps
1265 * must be added 180 degrees.
1266 * houses will be in clockwise direction */
1267 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1268 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1269 if (acmc < 0) {
1270 hsp->ac = swe_degnorm(hsp->ac + 180);
1271 hsp->mc = swe_degnorm(hsp->mc + 180);
1272 for (i = 1; i <= 12; i++) {
1273 if (i >= 4 && i < 10) continue;
1274 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1275 }
1276 }
1277 }
1278 break;
1279 case 'S': /* Sripati houses */
1280 /* uses Porphyry sectors, but then takes middle of sectors as cusps */
1281 {
1282 double s1, s4, q1;
1283 acmc = swe_difdeg2n(hsp->ac, hsp->mc); // size of 4th quadrant
1284 if (acmc < 0) {
1285 /* within polar circle we swap AC/DC if AC is on wrong side */
1286 hsp->ac = swe_degnorm(hsp->ac + 180);
1287 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1288 }
1289 q1 = 180 - acmc; // size of 1st quadrant
1290 s1 = q1 / 3.0;
1291 s4 = acmc / 3.0;
1292 hsp->cusp[1] = swe_degnorm(hsp->ac - s4 * 0.5);
1293 hsp->cusp[2] = swe_degnorm(hsp->ac + s1 * 0.5);
1294 hsp->cusp[3] = swe_degnorm(hsp->ac + s1 * 1.5);
1295 hsp->cusp[10] = swe_degnorm(hsp->mc - s1 * 0.5);
1296 hsp->cusp[11] = swe_degnorm(hsp->mc + s4 * 0.5);
1297 hsp->cusp[12] = swe_degnorm(hsp->mc + s4 * 1.5);
1298 }
1299 hsp->do_interpol = hsp->do_hspeed;
1300 break;
1301 case 'T': /* 'topocentric' houses */
1302 fh1 = atand (tanfi / 3.0);
1303 fh2 = atand (tanfi * 2.0 / 3.0);
1304 hsp->cusp[11] = Asc1(30 + th, fh1, sine, cose);
1305 hsp->cusp[12] = Asc1(60 + th, fh2, sine, cose);
1306 hsp->cusp[2] = Asc1(120 + th, fh2, sine, cose);
1307 hsp->cusp[3] = Asc1(150 + th, fh1, sine, cose);
1308 if (hsp->do_hspeed) {
1309 hsp->cusp_speed[11] = AscDash(30 + th, fh1, sine, cose);
1310 hsp->cusp_speed[12] = AscDash(60 + th, fh2, sine, cose);
1311 hsp->cusp_speed[2] = AscDash(120 + th, fh2, sine, cose);
1312 hsp->cusp_speed[3] = AscDash(150 + th, fh1, sine, cose);
1313 }
1314 /* within polar circle, when mc sinks below horizon and
1315 * ascendant changes to western hemisphere, all cusps
1316 * must be added 180 degrees.
1317 * houses will be in clockwise direction */
1318 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1319 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1320 if (acmc < 0) {
1321 hsp->ac = swe_degnorm(hsp->ac + 180);
1322 hsp->mc = swe_degnorm(hsp->mc + 180);
1323 for (i = 1; i <= 12; i++)
1324 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1325 }
1326 }
1327 break;
1328 case 'V': /* equal houses after Vehlow */
1329 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1330 if (acmc < 0) {
1331 /* within polar circle we swap AC/DC if AC is on wrong side */
1332 hsp->ac = swe_degnorm(hsp->ac + 180);
1333 }
1334 hsp->cusp[1] = swe_degnorm(hsp->ac - 15);
1335 for (i = 2; i <=12; i++)
1336 hsp->cusp[i] = swe_degnorm(hsp->cusp[1] + (i-1) * 30);
1337 if (hsp->do_hspeed) {
1338 for (i = 1; i <=12; i++) {
1339 hsp->cusp_speed[i] = hsp->ac_speed;
1340 }
1341 }
1342 break;
1343 case 'W': /* equal, whole-sign houses */
1344 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1345 if (acmc < 0) {
1346 /* within polar circle we swap AC/DC if AC is on wrong side */
1347 hsp->ac = swe_degnorm(hsp->ac + 180);
1348 hsp->cusp[1] = hsp->ac;
1349 }
1350 hsp->cusp[1] = hsp->ac - fmod(hsp->ac, 30);
1351 for (i = 2; i <=12; i++)
1352 hsp->cusp[i] = swe_degnorm(hsp->cusp[1] + (i-1) * 30);
1353 break;
1354 case 'X': {
1355 /*
1356 * Meridian or axial rotation system:
1357 * ecliptic points whose rectascensions
1358 * are armc + n * 30
1359 */
1360 int j;
1361 double a = th;
1362 for (i = 1; i <= 12; i++) {
1363 j = i + 10;
1364 if (j > 12) j -= 12;
1365 a = swe_degnorm(a + 30);
1366 if (fabs(a - 90) > VERY_SMALL
1367 && fabs(a - 270) > VERY_SMALL) {
1368 tant = tand(a);
1369 hsp->cusp[j] = atand(tant / cose);
1370 if (a > 90 && a <= 270)
1371 hsp->cusp[j] = swe_degnorm(hsp->cusp[j] + 180);
1372 } else {
1373 if (fabs(a - 90) <= VERY_SMALL)
1374 hsp->cusp[j] = 90;
1375 else
1376 hsp->cusp[j] = 270;
1377 } /* if */
1378 hsp->cusp[j] = swe_degnorm(hsp->cusp[j]);
1379 }
1380 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1381 if (acmc < 0) {
1382 hsp->ac = swe_degnorm(hsp->ac + 180);
1383 }
1384 hsp->do_interpol = hsp->do_hspeed;
1385 break; }
1386 case 'M': {
1387 /*
1388 * Morinus
1389 * points of the equator (armc + n * 30) are transformed
1390 * into the ecliptic coordinate system
1391 */
1392 int j;
1393 double a = th;
1394 double x[3];
1395 for (i = 1; i <= 12; i++) {
1396 j = i + 10;
1397 if (j > 12) j -= 12;
1398 a = swe_degnorm(a + 30);
1399 x[0] = a;
1400 x[1] = 0;
1401 swe_cotrans(x, x, ekl);
1402 hsp->cusp[j] = x[0];
1403 }
1404 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1405 if (acmc < 0) {
1406 hsp->ac = swe_degnorm(hsp->ac + 180);
1407 }
1408 hsp->do_interpol = hsp->do_hspeed;
1409 break; }
1410 case 'F': {
1411 /*
1412 * Carter poli-equatorial
1413 * Rectascension a of ascendant is the starting point.
1414 * house cusps nh on the ecliptic are the points where
1415 * great circles through points of the equator (a + (nh -1) * 30)
1416 * and the poles intersect it.
1417 */
1418 double a, ra;
1419 double x[3];
1420 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1421 if (acmc < 0) {
1422 /* within polar circle we swap AC/DC if AC is on wrong side */
1423 hsp->ac = swe_degnorm(hsp->ac + 180);
1424 hsp->cusp[1] = hsp->ac;
1425 }
1426 x[0] = hsp->ac;
1427 x[1] = 0;
1428 swe_cotrans(x, x, -ekl);
1429 a = x[0]; /* rectascension of ascendant */
1430 for (i = 2; i <= 12; i++) {
1431 if (i <= 3 || i >= 10) {
1432 ra = swe_degnorm(a + (i - 1) * 30);
1433 if (fabs(ra - 90) > VERY_SMALL
1434 && fabs(ra - 270) > VERY_SMALL) {
1435 tant = tand(ra);
1436 hsp->cusp[i] = atand(tant / cose);
1437 if (ra > 90 && ra <= 270)
1438 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1439 } else {
1440 if (fabs(ra - 90) <= VERY_SMALL)
1441 hsp->cusp[i] = 90;
1442 else
1443 hsp->cusp[i] = 270;
1444 } /* if */
1445 hsp->cusp[i] = swe_degnorm(hsp->cusp[i]);
1446 }
1447 }
1448 hsp->do_interpol = hsp->do_hspeed;
1449 break; }
1450 case 'B': { /* Alcabitius */
1451 /* created by Alois 17-sep-2000, followed example in Matrix
1452 electrical library. The code reproduces the example!
1453 I think the Alcabitius code in Walter Pullen's Astrolog 5.40
1454 is wrong, because he remains in RA and forgets the transform to
1455 the ecliptic. */
1456 double dek, r, sna, sda, sn3, sd3;
1457 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1458 if (acmc < 0) {
1459 hsp->ac = swe_degnorm(hsp->ac + 180);
1460 hsp->cusp[1] = hsp->ac;
1461 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1462 }
1463 dek = asind(sind(hsp->ac) * sine); /* declination of Ascendant */
1464 /* must treat the case fi == 90 or -90 */
1465 r = -tanfi * tand(dek);
1466 /* must treat the case of abs(r) > 1; probably does not happen
1467 * because dek becomes smaller when fi is large, as ac is close to
1468 * zero Aries/Libra in that case.
1469 */
1470 sda = acos(r) * RADTODEG; /* semidiurnal arc, measured on equator */
1471 sna = 180 - sda; /* complement, seminocturnal arc */
1472 sd3 = sda / 3;
1473 sn3 = sna / 3;
1474 rectasc = swe_degnorm(th + sd3); /* cusp 11 */
1475 /* project rectasc onto eclipitic with pole height 0, i.e. along the
1476 declination circle */
1477 hsp->cusp[11] = Asc1(rectasc, 0, sine, cose);
1478 rectasc = swe_degnorm(th + 2 * sd3); /* cusp 12 */
1479 hsp->cusp[12] = Asc1(rectasc, 0, sine, cose);
1480 rectasc = swe_degnorm(th + 180 - 2 * sn3); /* cusp 2 */
1481 hsp->cusp[2] = Asc1(rectasc, 0, sine, cose);
1482 rectasc = swe_degnorm(th + 180 - sn3); /* cusp 3 */
1483 hsp->cusp[3] = Asc1(rectasc, 0, sine, cose);
1484 }
1485 hsp->do_interpol = hsp->do_hspeed;
1486 break;
1487 case 'G': /* 36 Gauquelin sectors */
1488 for (i = 1; i <= 36; i++) {
1489 hsp->cusp[i] = 0;
1490 hsp->cusp_speed[i] = 0;
1491 }
1492 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1493 retc = ERR;
1494 strcpy(hsp->serr, "within polar circle, switched to Porphyry");
1495 hsy = (int) 'O';
1496 goto porphyry;
1497 }
1498 /*************** forth/second quarter ***************/
1499 /* note: Gauquelin sectors are counted in clockwise direction */
1500 a = asind(tand(fi) * tane);
1501 for (ih = 2; ih <= 9; ih++) {
1502 ih2 = 10 - ih;
1503 fh1 = atand(sind(a * ih2 / 9) / tane);
1504 rectasc = swe_degnorm((90 / 9) * ih2 + th);
1505 tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose))));
1506 if (fabs(tant) < VERY_SMALL) {
1507 hsp->cusp[ih] = rectasc;
1508 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1509 } else {
1510 /* pole height */
1511 f = atand(sind(asind(tanfi * tant) * ih2 / 9) /tant);
1512 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1513 cuspsv = 0;
1514 for (i = 1; i <= niter_max; i++) {
1515 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1516 if (fabs(tant) < VERY_SMALL) {
1517 hsp->cusp[ih] = rectasc;
1518 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1519 break;
1520 }
1521 /* pole height */
1522 f = atand(sind(asind(tanfi * tant) * ih2 / 9) / tant);
1523 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1524 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1525 break;
1526 cuspsv = hsp->cusp[ih];
1527 }
1528 #ifdef DEBUG_PLAC_ITER
1529 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1530 #endif
1531 if (i >= niter_max) {
1532 retc = ERR;
1533 hsy = (int) 'O';
1534 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1535 goto porphyry;
1536 }
1537 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1538 }
1539 hsp->cusp[ih+18] = swe_degnorm(hsp->cusp[ih] + 180);
1540 if (hsp->do_hspeed) hsp->cusp_speed[ih + 18] = hsp->cusp_speed[ih];
1541 }
1542 /*************** first/third quarter ***************/
1543 for (ih = 29; ih <= 36; ih++) {
1544 ih2 = ih - 28;
1545 fh1 = atand(sind(a * ih2 / 9) / tane);
1546 rectasc = swe_degnorm(180 - ih2 * 90 / 9 + th);
1547 tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose))));
1548 if (fabs(tant) < VERY_SMALL) {
1549 hsp->cusp[ih] = rectasc;
1550 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1551 } else {
1552 f = atand(sind(asind(tanfi * tant) * ih2 / 9) / tant);
1553 /* pole height */
1554 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1555 cuspsv = 0;
1556 for (i = 1; i <= niter_max; i++) {
1557 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1558 if (fabs(tant) < VERY_SMALL) {
1559 hsp->cusp[ih] = rectasc;
1560 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1561 break;
1562 }
1563 f = atand(sind(asind(tanfi * tant) * ih2 / 9) / tant);
1564 /* pole height */
1565 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1566 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1567 break;
1568 cuspsv = hsp->cusp[ih];
1569 }
1570 #ifdef DEBUG_PLAC_ITER
1571 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1572 #endif
1573 if (i >= niter_max) {
1574 retc = ERR;
1575 hsy = (int) 'O';
1576 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1577 goto porphyry;
1578 }
1579 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1580 }
1581 hsp->cusp[ih-18] = swe_degnorm(hsp->cusp[ih] + 180);
1582 if (hsp->do_hspeed) hsp->cusp_speed[ih - 18] = hsp->cusp_speed[ih];
1583 }
1584 hsp->cusp[1] = hsp->ac;
1585 hsp->cusp[10] = hsp->mc;
1586 hsp->cusp[19] = swe_degnorm(hsp->ac + 180);
1587 hsp->cusp[28] = swe_degnorm(hsp->mc + 180);
1588 if (hsp->do_hspeed) {
1589 hsp->cusp_speed[1] = hsp->ac_speed;
1590 hsp->cusp_speed[10] = hsp->mc_speed;
1591 hsp->cusp_speed[19] = hsp->ac_speed;
1592 hsp->cusp_speed[28] = hsp->mc_speed;
1593 }
1594 break;
1595 case 'U': /* Krusinski-Pisa */
1596 /*
1597 * The following code was written by Bogdan Krusinski in 2006.
1598 * bogdan@astrologia.pl
1599 *
1600 * Definition:
1601 * "Krusinski - house system based on the great circle passing through
1602 * ascendant and zenith. This circle is divided into 12 equal parts
1603 * (1st cusp is ascendent, 10th cusp is zenith), then the resulting
1604 * points are projected onto the ecliptic through meridian circles.
1605 * The house cusps in space are half-circles perpendicular to the equator
1606 * and running from the north to the south celestial pole through the
1607 * resulting cusp points on the house circle. The points where they
1608 * cross the ecliptic mark the ecliptic house cusps."
1609 *
1610 * Description of the algorithm:
1611 * Transform into great circle running through Asc and zenit (where arc
1612 * between Asc and zenith is always 90 deg), and then return with
1613 * house cusps into ecliptic. Eg. solve trigonometrical triangle
1614 * with three transformations and two rotations starting from ecliptic.
1615 * House cusps in space are meridian circles.
1616 *
1617 * Notes:
1618 * 1. In this definition we assume MC on ecliptic as point where
1619 * half-meridian (from north to south pole) cuts ecliptic,
1620 * so MC may be below horizon in arctic regions.
1621 * 2. Houses could be calculated in all latitudes except the poles
1622 * themselves (-90,90) and points on arctic circle in cases where
1623 * ecliptic is equal to horizon and then ascendant is undefined.
1624 * But ascendant when 'horizon=ecliptic' could be deduced as limes
1625 * from both sides of that point and houses with that provision can
1626 * be computed also there.
1627 *
1628 * Starting values for calculations:
1629 * - Asc ecliptic longitude
1630 * - right ascension of MC (RAMC)
1631 * - geographic latitude.
1632 */
1633 /*
1634 * within polar circle we swap AC/DC if AC is on wrong side
1635 */
1636 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1637 if (acmc < 0) {
1638 hsp->ac = swe_degnorm(hsp->ac + 180);
1639 }
1640 /* A0. Start point - ecliptic coords of ascendant */
1641 x[0] = hsp->ac; /* Asc longitude */
1642 x[1] = 0.0; /* Asc declination */
1643 x[2] = 1.0; /* Radius to test validity of subsequent transformations. */
1644 swe_cotrans(x, x, -ekl); /* A1. Transform into equatorial coords */
1645 x[0] = x[0] - (th-90); /* A2. Rotate */
1646 swe_cotrans(x, x, -(90-fi)); /* A3. Transform into horizontal coords */
1647 krHorizonLon = x[0]; /* ...save asc lon on horizon to get back later with house cusp */
1648 x[0] = x[0] - x[0]; /* A4. Rotate */
1649 swe_cotrans(x, x, -90); /* A5. Transform into this house system great circle (asc-zenith) */
1650 /* As it is house circle now, simple add 30 deg increments... */
1651 for(i = 0; i < 6; i++) {
1652 /* B0. Set 'n-th' house cusp.
1653 * Note that IC/MC are also calculated here to check
1654 * if really this is the asc-zenith great circle. */
1655 x[0] = 30.0*i;
1656 x[1] = 0.0;
1657 swe_cotrans(x, x, 90); /* B1. Transform back into horizontal coords */
1658 x[0] = x[0] + krHorizonLon; /* B2. Rotate back. */
1659 swe_cotrans(x, x, 90-fi); /* B3. Transform back into equatorial coords */
1660 x[0] = swe_degnorm(x[0] + (th-90)); /* B4. Rotate back -> RA of house cusp as result. */
1661 /* B5. Where's this house cusp on ecliptic? */
1662 /* ... so last but not least - get ecliptic longitude of house cusp: */
1663 hsp->cusp[i+1] = atand(tand(x[0])/cosd(ekl));
1664 if (x[0] > 90 && x[0] <= 270)
1665 hsp->cusp[i+1] = swe_degnorm(hsp->cusp[i+1] + 180);
1666 hsp->cusp[i+1] = swe_degnorm(hsp->cusp[i+1]);
1667 hsp->cusp[i+7] = swe_degnorm(hsp->cusp[i+1]+180);
1668 }
1669 break;
1670 case 'Y': /* APC houses */
1671 for (i = 1; i <= 12; i++) {
1672 hsp->cusp[i] = apc_sector(i, fi * DEGTORAD, ekl * DEGTORAD, th * DEGTORAD);
1673 }
1674 //hsp->ac = hsp->cusp[1];
1675 //hsp->mc = hsp->cusp[10];
1676 /* note the MC provided by apc_sector() near latitude 90 is not accurate */
1677 hsp->cusp[10] = hsp->mc;
1678 hsp->cusp[4] = swe_degnorm(hsp->mc + 180);
1679 /* within polar circle, when mc sinks below horizon and
1680 * ascendant changes to western hemisphere, all cusps
1681 * must be added 180 degrees.
1682 * houses will be in clockwise direction */
1683 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1684 acmc = swe_difdeg2n(hsp->ac, hsp->mc);
1685 if (acmc < 0) {
1686 hsp->ac = swe_degnorm(hsp->ac + 180);
1687 hsp->mc = swe_degnorm(hsp->mc + 180);
1688 for (i = 1; i <= 12; i++)
1689 hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180);
1690 }
1691 }
1692 hsp->do_interpol = hsp->do_hspeed;
1693 break;
1694 default: /* Placidus houses */
1695 if (fabs(fi) >= 90 - ekl) { /* within polar circle */
1696 retc = ERR;
1697 strcpy(hsp->serr, "within polar circle, switched to Porphyry");
1698 goto porphyry;
1699 }
1700 a = asind(tand(fi) * tane);
1701 fh1 = atand(sind(a / 3) / tane);
1702 fh2 = atand(sind(a * 2 / 3) / tane);
1703 /* ************ house 11 ******************** */
1704 rectasc = swe_degnorm(30 + th);
1705 tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose))));
1706 ih = 11;
1707 if (fabs(tant) < VERY_SMALL) {
1708 hsp->cusp[ih] = rectasc;
1709 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1710 } else {
1711 /* pole height */
1712 f = atand(sind(asind(tanfi * tant) / 3) /tant);
1713 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1714 cuspsv = 0;
1715 for (i = 1; i <= niter_max; i++) {
1716 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1717 if (fabs(tant) < VERY_SMALL) {
1718 hsp->cusp[ih] = rectasc;
1719 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1720 break;
1721 }
1722 /* pole height */
1723 f = atand(sind(asind(tanfi * tant) / 3) / tant);
1724 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1725 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1726 break;
1727 cuspsv = hsp->cusp[ih];
1728 }
1729 if (i >= niter_max) {
1730 retc = ERR;
1731 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1732 goto porphyry;
1733 }
1734 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1735 #ifdef DEBUG_PLAC_ITER
1736 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1737 #endif
1738 }
1739 /* ************ house 12 ******************** */
1740 rectasc = swe_degnorm(60 + th);
1741 tant = tand(asind(sine*sind(Asc1(rectasc, fh2, sine, cose))));
1742 ih = 12;
1743 if (fabs(tant) < VERY_SMALL) {
1744 hsp->cusp[ih] = rectasc;
1745 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1746 } else {
1747 f = atand(sind(asind(tanfi * tant) / 1.5) / tant);
1748 /* pole height */
1749 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1750 cuspsv = 0;
1751 for (i = 1; i <= niter_max; i++) {
1752 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1753 if (fabs(tant) < VERY_SMALL) {
1754 hsp->cusp[ih] = rectasc;
1755 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1756 break;
1757 }
1758 f = atand(sind(asind(tanfi * tant) / 1.5) / tant);
1759 /* pole height */
1760 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1761 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1762 break;
1763 cuspsv = hsp->cusp[ih];
1764 }
1765 if (i >= niter_max) {
1766 retc = ERR;
1767 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1768 goto porphyry;
1769 }
1770 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1771 #ifdef DEBUG_PLAC_ITER
1772 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1773 #endif
1774 }
1775 /* ************ house 2 ******************** */
1776 rectasc = swe_degnorm(120 + th);
1777 tant = tand(asind(sine * sind(Asc1(rectasc, fh2, sine, cose))));
1778 ih = 2;
1779 if (fabs(tant) < VERY_SMALL) {
1780 hsp->cusp[ih] = rectasc;
1781 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1782 } else {
1783 f = atand(sind(asind(tanfi * tant) / 1.5) / tant);
1784 /* pole height */
1785 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1786 cuspsv = 0;
1787 for (i = 1; i <= niter_max; i++) {
1788 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1789 if (fabs(tant) < VERY_SMALL) {
1790 hsp->cusp[ih] = rectasc;
1791 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1792 break;
1793 }
1794 f = atand(sind(asind(tanfi * tant) / 1.5) / tant);
1795 /* pole height */
1796 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1797 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1798 break;
1799 cuspsv = hsp->cusp[ih];
1800 }
1801 if (i >= niter_max) {
1802 retc = ERR;
1803 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1804 goto porphyry;
1805 }
1806 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1807 #ifdef DEBUG_PLAC_ITER
1808 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1809 #endif
1810 }
1811 /* ************ house 3 ******************** */
1812 rectasc = swe_degnorm(150 + th);
1813 tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose))));
1814 ih = 3;
1815 if (fabs(tant) < VERY_SMALL) {
1816 hsp->cusp[ih] = rectasc;
1817 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1818 } else {
1819 f = atand(sind(asind(tanfi * tant) / 3) / tant);
1820 /* pole height */
1821 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1822 cuspsv = 0;
1823 for (i = 1; i <= niter_max; i++) {
1824 tant = tand(asind(sine * sind(hsp->cusp[ih])));
1825 if (fabs(tant) < VERY_SMALL) {
1826 hsp->cusp[ih] = rectasc;
1827 if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed;
1828 break;
1829 }
1830 f = atand(sind(asind(tanfi * tant) / 3) / tant);
1831 /* pole height */
1832 hsp->cusp[ih] = Asc1(rectasc, f, sine, cose);
1833 if (i > 1 && fabs(swe_difdeg2n(hsp->cusp[ih], cuspsv)) < VERY_SMALL_PLAC_ITER)
1834 break;
1835 cuspsv = hsp->cusp[ih];
1836 }
1837 if (i >= niter_max) {
1838 retc = ERR;
1839 strcpy(hsp->serr, "very close to polar circle, switched to Porphyry");
1840 goto porphyry;
1841 }
1842 if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose);
1843 #ifdef DEBUG_PLAC_ITER
1844 fprintf(stderr, "h=%d, niter=%d\n", ih, i);
1845 #endif
1846 }
1847 break;
1848 } /* end switch */
1849 if (hsy != 'G' && hsy != 'Y' && toupper(hsy) != 'I') {
1850 hsp->cusp[4] = swe_degnorm(hsp->cusp[10] + 180);
1851 hsp->cusp[5] = swe_degnorm(hsp->cusp[11] + 180);
1852 hsp->cusp[6] = swe_degnorm(hsp->cusp[12] + 180);
1853 hsp->cusp[7] = swe_degnorm(hsp->cusp[1] + 180);
1854 hsp->cusp[8] = swe_degnorm(hsp->cusp[2] + 180);
1855 hsp->cusp[9] = swe_degnorm(hsp->cusp[3] + 180);
1856 if (hsp->do_hspeed && ! hsp->do_interpol) {
1857 hsp->cusp_speed[4] = hsp->cusp_speed[10];
1858 hsp->cusp_speed[5] = hsp->cusp_speed[11];
1859 hsp->cusp_speed[6] = hsp->cusp_speed[12];
1860 hsp->cusp_speed[7] = hsp->cusp_speed[1];
1861 hsp->cusp_speed[8] = hsp->cusp_speed[2];
1862 hsp->cusp_speed[9] = hsp->cusp_speed[3];
1863 }
1864 }
1865 /* vertex */
1866 if (fi >= 0)
1867 f = 90 - fi;
1868 else
1869 f = -90 - fi;
1870 hsp->vertex = Asc1(th - 90, f, sine, cose);
1871 if (hsp->do_speed) hsp->vertex_speed = AscDash(th - 90, f, sine, cose);
1872 /* with tropical latitudes, the vertex behaves strange,
1873 * in a similar way as the ascendant within the polar
1874 * circle. we keep it always on the western hemisphere.*/
1875 if (fabs(fi) <= ekl) {
1876 vemc = swe_difdeg2n(hsp->vertex, hsp->mc);
1877 if (vemc > 0)
1878 hsp->vertex = swe_degnorm(hsp->vertex + 180);
1879 }
1880 /*
1881 * some strange points:
1882 */
1883 /* equasc (equatorial ascendant) */
1884 th2 = swe_degnorm(th + 90);
1885 if (fabs(th2 - 90) > VERY_SMALL
1886 && fabs(th2 - 270) > VERY_SMALL) {
1887 tant = tand(th2);
1888 hsp->equasc = atand(tant / cose);
1889 if (th2 > 90 && th2 <= 270)
1890 hsp->equasc = swe_degnorm(hsp->equasc + 180);
1891 } else {
1892 if (fabs(th2 - 90) <= VERY_SMALL)
1893 hsp->equasc = 90;
1894 else
1895 hsp->equasc = 270;
1896 } /* if */
1897 hsp->equasc = swe_degnorm(hsp->equasc);
1898 if (hsp->do_speed) hsp->equasc_speed = AscDash(th + 90, 0, sine, cose);
1899 /* "co-ascendant" W. Koch */
1900 hsp->coasc1 = swe_degnorm(Asc1(th - 90, fi, sine, cose) + 180);
1901 if (hsp->do_speed) hsp->coasc1_speed = AscDash(th - 90, fi, sine, cose);
1902 /* "co-ascendant" M. Munkasey */
1903 if (fi >= 0) {
1904 hsp->coasc2 = Asc1(th + 90, 90 - fi, sine, cose);
1905 if (hsp->do_speed) hsp->coasc2_speed = AscDash(th + 90, 90 - fi, sine, cose);
1906 } else { /* southern hemisphere */
1907 hsp->coasc2 = Asc1(th + 90, -90 - fi, sine, cose);
1908 if (hsp->do_speed) hsp->coasc2_speed = AscDash(th + 90, -90 - fi, sine, cose);
1909 }
1910 /* "polar ascendant" M. Munkasey */
1911 hsp->polasc = Asc1(th - 90, fi, sine, cose);
1912 if (hsp->do_speed) hsp->polasc_speed = AscDash(th - 90, fi, sine, cose);
1913 #if 0
1914 test_Asc1();
1915 #endif
1916 return retc;
1917 } /* procedure houses */
1918
1919 /******************************/
Asc1(double x1,double f,double sine,double cose)1920 static double Asc1(double x1, double f, double sine, double cose)
1921 {
1922 int n;
1923 double ass;
1924 x1 = swe_degnorm(x1);
1925 n = (int) ((x1 / 90) + 1); // n is quadrant 1..4
1926 if (fabs(90 - f) < VERY_SMALL) { // near north pole
1927 return 180;
1928 }
1929 if (fabs(90 + f) < VERY_SMALL) { // near south pole
1930 return 0;
1931 }
1932 if (n == 1)
1933 ass = ( Asc2(x1, f, sine, cose));
1934 else if (n == 2)
1935 ass = (180 - Asc2(180 - x1, - f, sine, cose));
1936 else if (n == 3)
1937 ass = (180 + Asc2(x1 - 180, - f, sine, cose));
1938 else
1939 ass = (360 - Asc2(360- x1, f, sine, cose));
1940 ass = swe_degnorm(ass);
1941 if (fabs(ass - 90) < VERY_SMALL) /* rounding, e.g.: if */
1942 ass = 90; /* fi = 0 & st = 0, ac = 89.999... */
1943 if (fabs(ass - 180) < VERY_SMALL)
1944 ass = 180;
1945 if (fabs(ass - 270) < VERY_SMALL) /* rounding, e.g.: if */
1946 ass = 270; /* fi = 0 & st = 0, ac = 89.999... */
1947 if (fabs(ass - 360) < VERY_SMALL)
1948 ass = 0;
1949 return ass;
1950 } /* Asc1 */
1951
1952 // derivative of Asc1, computes speed
1953 // code crontributed by Graham Dawson
AscDash(double x,double f,double sine,double cose)1954 static double AscDash(double x, double f, double sine, double cose)
1955 {
1956 double cosx = cosd(x);
1957 double sinx = sind(x);
1958 double sinx2 = sinx * sinx;
1959 double c = cose * cosx - tand(f) * sine;
1960 double d = sinx2 + c * c;
1961 double dudt;
1962 if (d > VERY_SMALL) {
1963 dudt = (cosx * c + cose * sinx2) / d;
1964 } else {
1965 dudt = 0.0; // When we are on axis of ecliptic
1966 }
1967 return dudt * ARMCS; // 360.985647366;
1968 }
1969
1970
1971 #if 0
1972 /******************************/
1973 static double Asc1_old(double x1, double f, double sine, double cose)
1974 {
1975 int n;
1976 double ass;
1977 if (f == -90) f += VERY_SMALL / 1000; // avoid exact pole 90, as tan() goes infinite
1978 if (f == 90) f -= VERY_SMALL / 1000;
1979 x1 = swe_degnorm(x1);
1980 n = (int) ((x1 / 90) + 1);
1981 if (n == 1)
1982 ass = ( Asc2(x1, f, sine, cose));
1983 else if (n == 2)
1984 ass = (180 - Asc2(180 - x1, - f, sine, cose));
1985 else if (n == 3)
1986 ass = (180 + Asc2(x1 - 180, - f, sine, cose));
1987 else
1988 ass = (360 - Asc2(360- x1, f, sine, cose));
1989 ass = swe_degnorm(ass);
1990 if (fabs(ass - 90) < VERY_SMALL) /* rounding, e.g.: if */
1991 ass = 90; /* fi = 0 & st = 0, ac = 89.999... */
1992 if (fabs(ass - 180) < VERY_SMALL)
1993 ass = 180;
1994 if (fabs(ass - 270) < VERY_SMALL) /* rounding, e.g.: if */
1995 ass = 270; /* fi = 0 & st = 0, ac = 89.999... */
1996 if (fabs(ass - 360) < VERY_SMALL)
1997 ass = 0;
1998 return ass;
1999 } /* Asc1 */
2000
2001 /******************************/
2002 static double Asc1_old_old (double x1, double f, double sine, double cose)
2003 {
2004 int n;
2005 double ass;
2006 x1 = swe_degnorm(x1);
2007 n = (int) ((x1 / 90) + 1);
2008 if (n == 1)
2009 ass = ( Asc2 (x1, f, sine, cose));
2010 else if (n == 2)
2011 ass = (180 - Asc2 (180 - x1, - f, sine, cose));
2012 else if (n == 3)
2013 ass = (180 + Asc2 (x1 - 180, - f, sine, cose));
2014 else
2015 ass = (360 - Asc2 (360- x1, f, sine, cose));
2016 ass = swe_degnorm(ass);
2017 if (fabs(ass - 90) < VERY_SMALL) /* rounding, e.g.: if */
2018 ass = 90; /* fi = 0 & st = 0, ac = 89.999... */
2019 if (fabs(ass - 180) < VERY_SMALL)
2020 ass = 180;
2021 if (fabs(ass - 270) < VERY_SMALL) /* rounding, e.g.: if */
2022 ass = 270; /* fi = 0 & st = 0, ac = 89.999... */
2023 if (fabs(ass - 360) < VERY_SMALL)
2024 ass = 0;
2025 return ass;
2026 } /* Asc1 */
2027
2028 static void test_Asc1()
2029 {
2030 double armc, dlat, eps = 23.44, asc1, asc1_old, sine, cose;
2031 sine = sind(eps);
2032 cose = cosd(eps);
2033 fprintf(stderr, "Test Asc1() <-> Asc1_old()\n");
2034 for (dlat = -90; dlat <= 90; dlat++) {
2035 for (armc = 0; armc <= 360; armc++) {
2036 asc1 = Asc1(armc, dlat, sine, cose);
2037 asc1_old = Asc1_old_old(armc, dlat, sine, cose);
2038 if (asc1 != asc1_old)
2039 fprintf(stderr, "armc=%f, lat=%f, Asc1: %.16f <-> Asc1_old: %.16f\n", armc, dlat, asc1, asc1_old);
2040 }
2041 }
2042
2043 }
2044 #endif
2045
2046 /*
2047 * x in range 0..90
2048 * f in range -90 .. +90
2049 * sine, cose around e=23°
2050 */
Asc2(double x,double f,double sine,double cose)2051 static double Asc2(double x, double f, double sine, double cose)
2052 {
2053 double ass, sinx;
2054 ass = - tand(f) * sine + cose * cosd(x);
2055 if (fabs(ass) < VERY_SMALL)
2056 ass = 0;
2057 sinx = sind(x);
2058 if (fabs(sinx) < VERY_SMALL)
2059 sinx = 0;
2060 if (sinx == 0) {
2061 if (ass < 0)
2062 ass = -VERY_SMALL;
2063 else
2064 ass = VERY_SMALL;
2065 } else if (ass == 0) {
2066 if (sinx < 0)
2067 ass = -90;
2068 else
2069 ass = 90;
2070 } else {
2071 ass = atand(sinx / ass);
2072 }
2073 if (ass < 0)
2074 ass = 180 + ass;
2075 return (ass);
2076 } /* Asc2 */
2077
armc_to_mc(double armc,double eps)2078 static double armc_to_mc(double armc, double eps)
2079 {
2080 double cose = cosd(eps);
2081 double mc, tant;
2082 if (fabs(armc - 90) > VERY_SMALL
2083 && fabs(armc - 270) > VERY_SMALL) {
2084 tant = tand(armc);
2085 mc = swe_degnorm(atand(tant / cose));
2086 if (armc > 90 && armc <= 270)
2087 mc = swe_degnorm(mc + 180);
2088 } else {
2089 if (fabs(armc - 90) <= VERY_SMALL)
2090 mc = 90;
2091 else
2092 mc = 270;
2093 }
2094 return mc;
2095 }
2096
2097 /* if ascendant is on western half of horizon, add 180° */
fix_asc_polar(double asc,double armc,double eps,double geolat)2098 static double fix_asc_polar(double asc, double armc, double eps, double geolat)
2099 {
2100 double demc = atand(sind(armc) * tand(eps));
2101 if (geolat >= 0 && 90 - geolat + demc < 0)
2102 asc = swe_degnorm(asc + 180);
2103 if (geolat < 0 && -90 - geolat + demc > 0)
2104 asc = swe_degnorm(asc + 180);
2105 return asc;
2106 }
2107
2108 /* Computes the house position of a planet or another point,
2109 * in degrees: 0 - 30 = 1st house, 30 - 60 = 2nd house, etc.
2110 * armc sidereal time in degrees
2111 * geolat geographic latitude
2112 * eps true ecliptic obliquity
2113 * hsys house system character
2114 * xpin array of 6 doubles:
2115 * only the first two of them are used: ecl. long., lat.
2116 * serr error message area
2117 *
2118 * House position is returned by function.
2119 * Currently, geometrically correct house positions are provided
2120 * for the following house methods:
2121 * A/E Equal, V Vehlow, W Whole Signs, D Equal/MC, N Equal/Zodiac,
2122 * O Porphyry, B Alcabitius, X Meridian, F Carter, M Morinus,
2123 * P Placidus, K Koch, C Campanus, R Regiomontanus, U Krusinski,
2124 * T Topocentric, H Horizon, G Gauquelin.
2125 *
2126 * A simplified house position (distance_from_cusp / house_size)
2127 * is currently provided for the following house methods:
2128 * Y APC houses, L Pullen SD, Q Pullen SR, I Sunshine, S Sripati.
2129 *
2130 * IMPORTANT: This function should NOT be used for sidereal astrology.
2131 * If you cannot avoid doing so, please note:
2132 * - The input longitudes (xpin) MUST always be tropical, even if you
2133 * are a siderealist.
2134 * - Sidereal and tropical house positions are identical for most house
2135 * systems, if a traditional definition of the sidereal zodiac is used
2136 * (sid = trop - ayanamsa).
2137 * - The function does NOT provide correct positions for Whole Sign houses.
2138 * - The function does NOT provide correct positions, if you use a
2139 * non-traditional sidereal method (where the sidereal plane is not
2140 * identical to the ecliptic of date) with a house system whose definition
2141 * is dependent on the ecliptic, such as:
2142 * equal, Porphyry, Alcabitius, Koch, Krusinski (all others should work).
2143 * The Swiss Ephemeris currently does not handle these cases.
2144 */
swe_house_pos(double armc,double geolat,double eps,int hsys,double * xpin,char * serr)2145 double CALL_CONV swe_house_pos(
2146 double armc, double geolat, double eps, int hsys, double *xpin, char *serr)
2147 {
2148 double xp[6], xeq[6], ra, de, mdd, mdn, sad, san;
2149 double hpos, sinad, ad, a, admc, adp, samc, asc, mc, acmc, tant;
2150 //double demc;
2151 double fh, ra0, tanfi, fac, dfac, tanx;
2152 double x[3], xasc[3], raep, raaz, oblaz, xtemp; /* BK 21.02.2006 */
2153 double hcusp[37], ascmc[10];
2154 double sine = sind(eps);
2155 double cose = cosd(eps);
2156 double c1, c2, d, hsize;
2157 int i, j, nloop;
2158 double dsun = 0, darmc, harmc, y, sinpsi, sa;
2159 AS_BOOL is_western_half = FALSE;
2160 hsys = toupper(hsys);
2161 if (1) {
2162 /* input is a house position: no calculation is required */
2163 ascmc[9] = 99;// dirty hack. Sunshine house system needs sun declination
2164 // which we do not know. If it sees ascmc[9] == 99, it uses
2165 // the one is saved from last call. can lead to bugs, but can
2166 // also solve many problems.
2167 if (swe_houses_armc_ex2(armc, geolat, eps, hsys, hcusp, ascmc, NULL, NULL, serr) == ERR) {
2168 if (serr != NULL)
2169 sprintf(serr, "swe_house_pos(): failed for system %c", hsys);
2170 } else {
2171 hpos = 0;
2172 for (i = 1; i <= 12; i++) {
2173 if (fabs(swe_difdeg2n(xpin[0], hcusp[i])) < MILLIARCSEC && xpin[1] == 0) {
2174 hpos = (double) i;
2175 }
2176 }
2177 for (i = 1; i <= 12; i += 3) {
2178 if (fabs(swe_difdeg2n(xpin[0], hcusp[i])) < MILLIARCSEC && xpin[1] == 0) {
2179 xp[0] = (double) i;
2180 }
2181 }
2182 if (hpos > 0)
2183 return hpos;
2184 // for Sunshine houses: declination of Sun
2185 if (hsys == 'I')
2186 dsun = ascmc[9];
2187 // for APC houses: declination of ascendant into dsun
2188 if (hsys == 'Y') {
2189 xeq[0] = ascmc[0];
2190 xeq[1] = 0;
2191 xeq[2] = 1;
2192 swe_cotrans(xeq, xeq, -eps);
2193 dsun = xeq[1];
2194 }
2195 }
2196 }
2197 AS_BOOL is_above_hor = FALSE;
2198 AS_BOOL is_invalid = FALSE;
2199 AS_BOOL is_circumpolar = FALSE;
2200 if (serr != NULL)
2201 *serr = '\0';
2202 xeq[0] = xpin[0];
2203 xeq[1] = xpin[1];
2204 xeq[2] = 1;
2205 swe_cotrans(xeq, xeq, -eps);
2206 ra = xeq[0];
2207 de = xeq[1];
2208 mdd = swe_degnorm(ra - armc);
2209 mdn = swe_degnorm(mdd + 180);
2210 if (mdd >= 180)
2211 mdd -= 360;
2212 if (mdn >= 180)
2213 mdn -= 360;
2214 /* xp[0] will contain the house position, a value between 0 and 360 */
2215 switch(hsys) {
2216 case 'N': // equal (1=Aries)
2217 xp[0] = xpin[0];
2218 hpos = xp[0] / 30.0 + 1;
2219 break;
2220 case 'A': // equal
2221 case 'E': // equal
2222 case 'D': // equal (MC)
2223 case 'V': // Vehlow
2224 case 'W': // whole signs
2225 asc = Asc1(swe_degnorm(armc + 90), geolat, sine, cose);
2226 mc = armc_to_mc(armc, eps);
2227 asc = fix_asc_polar(asc, armc, eps, geolat);
2228 xp[0] = swe_degnorm(xpin[0] - asc);
2229 if (hsys == 'V')
2230 xp[0] = swe_degnorm(xp[0] + 15);
2231 if (hsys == 'W')
2232 xp[0] = swe_degnorm(xp[0] + fmod(asc, 30));
2233 if (hsys == 'D')
2234 xp[0] = swe_degnorm(xpin[0] - mc - 90);
2235 /* to make sure that a call with a house cusp position returns
2236 * a value within the house, 0.001" is added */
2237 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2238 hpos = xp[0] / 30.0 + 1;
2239 break;
2240 case 'O': /* Porphyry */
2241 case 'B': /* Alcabitius */
2242 case 'S': /* Sripati */
2243 asc = Asc1(swe_degnorm(armc + 90), geolat, sine, cose);
2244 /* mc */
2245 mc = armc_to_mc(armc, eps);
2246 /* while MC is always south,
2247 * Asc must always be in eastern hemisphere */
2248 asc = fix_asc_polar(asc, armc, eps, geolat);
2249 if (hsys == 'O' || hsys == 'S') {
2250 xp[0] = swe_degnorm(xpin[0] - asc);
2251 /* to make sure that a call with a house cusp position returns
2252 * a value within the house, 0.001" is added */
2253 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2254 if (xp[0] < 180)
2255 hpos = 1;
2256 else {
2257 hpos = 7;
2258 xp[0] -= 180;
2259 }
2260 acmc = swe_difdeg2n(asc, mc);
2261 if (xp[0] < 180 - acmc)
2262 hpos += xp[0] * 3 / (180 - acmc);
2263 else
2264 hpos += 3 + (xp[0] - 180 + acmc) * 3 / acmc;
2265 if (hsys == 'S') {
2266 hpos += 0.5;
2267 if (hpos > 12) hpos = 1;
2268 }
2269 } else { /* Alcabitius */
2270 double dek, r, sna, sda;
2271 dek = asind(sind(asc) * sine); /* declination of Ascendant */
2272 /* must treat the case fi == 90 or -90 */
2273 tanfi = tand(geolat);
2274 r = -tanfi * tand(dek);
2275 /* must treat the case of abs(r) > 1; probably does not happen
2276 * because dek becomes smaller when fi is large, as ac is close to
2277 * zero Aries/Libra in that case.
2278 */
2279 sda = acos(r) * RADTODEG; /* semidiurnal arc, measured on equator */
2280 sna = 180 - sda; /* complement, seminocturnal arc */
2281 if (mdd > 0) {
2282 if (mdd < sda)
2283 hpos = mdd * 90 / sda;
2284 else
2285 hpos = 90 + (mdd - sda) * 90 / sna;
2286 } else {
2287 if (mdd > -sna)
2288 hpos = 360 + mdd * 90 / sna;
2289 else
2290 hpos = 270 + (mdd + sna) * 90 / sda;
2291 }
2292 hpos = swe_degnorm(hpos - 90) / 30.0 + 1.0;
2293 if (hpos >= 13.0) hpos -= 12;
2294 }
2295 break;
2296 case 'X': /* Meridian or axial rotation system */
2297 hpos = swe_degnorm(mdd - 90) / 30.0 + 1.0;
2298 break;
2299 case 'F': /* Carter poli-equatorial */
2300 x[0] = Asc1(swe_degnorm(armc + 90), geolat, sine, cose);
2301 x[0] = fix_asc_polar(x[0], armc, eps, geolat);
2302 x[1] = 0;
2303 swe_cotrans(x, x, -eps);
2304 hpos = swe_degnorm(ra - x[0]) / 30.0 + 1;
2305 break;
2306 case 'M': { /* Morinus */
2307 double a = xpin[0];
2308 if (fabs(a - 90) > VERY_SMALL
2309 && fabs(a - 270) > VERY_SMALL) {
2310 tant = tand(a);
2311 hpos = atand(tant / cose);
2312 if (a > 90 && a <= 270)
2313 hpos = swe_degnorm(hpos + 180);
2314 } else {
2315 if (fabs(a - 90) <= VERY_SMALL)
2316 hpos = 90;
2317 else
2318 hpos = 270;
2319 } /* if */
2320 hpos = swe_degnorm(hpos - armc - 90);
2321 hpos = hpos / 30.0 + 1;
2322 }
2323 break;
2324 #if 0
2325 /* old version of Koch method */
2326 case 'K':
2327 demc = atand(sind(armc) * tand(eps));
2328 /* if body is within circumpolar region, error */
2329 if (90 - fabs(geolat) <= fabs(de)) {
2330 if (serr != NULL)
2331 strcpy(serr, "no Koch house position, because planet is circumpolar.");
2332 xp[0] = 0;
2333 hpos = 0; /* Error */
2334 } else if (90 - fabs(geolat) <= fabs(demc)) {
2335 if (serr != NULL)
2336 strcpy(serr, "no Koch house position, because mc is circumpolar.");
2337 xp[0] = 0;
2338 hpos = 0; /* Error */
2339 } else {
2340 admc = asind(tand(eps) * tand(geolat) * sind(armc));
2341 adp = asind(tand(geolat) * tand(de));
2342 samc = 90 + admc;
2343 if (mdd >= 0) { /* east */
2344 xp[0] = swe_degnorm(((mdd - adp + admc) / samc - 1) * 90);
2345 } else {
2346 xp[0] = swe_degnorm(((mdd + 180 + adp + admc) / samc + 1) * 90);
2347 }
2348 /* to make sure that a call with a house cusp position returns
2349 * a value within the house, 0.001" is added */
2350 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2351 hpos = xp[0] / 30.0 + 1;
2352 }
2353 break;
2354 #endif
2355 /* version of Koch method: do calculations within circumpolar circle,
2356 * if possible; make sure house positions 4 - 9 only appear on western
2357 * hemisphere */
2358 case 'K': // Koch
2359 //demc = atand(sind(armc) * tand(eps));
2360 is_invalid = FALSE;
2361 is_circumpolar = FALSE;
2362 /* object is within a circumpolar circle */
2363 if (90 - geolat < de || -90 - geolat > de) {
2364 adp = 90;
2365 is_circumpolar = TRUE;
2366 }
2367 /* object is within a circumpolar circle, southern hemisphere */
2368 else if (geolat - 90 > de || geolat + 90 < de) {
2369 adp = -90;
2370 is_circumpolar = TRUE;
2371 }
2372 /* object does rise and set */
2373 else {
2374 adp = asind(tand(geolat) * tand(de));
2375 }
2376 #if 0
2377 if (fabs(adp) == 90)
2378 is_invalid = TRUE; /* omit this to use the above values */
2379 #endif
2380 admc = tand(eps) * tand(geolat) * sind(armc);
2381 /* midheaven is circumpolar */
2382 if (fabs(admc) > 1) {
2383 #if 0
2384 is_invalid = TRUE; /* omit this line to use the below values */
2385 #endif
2386 if (admc > 1)
2387 admc = 1;
2388 else
2389 admc = -1;
2390 is_circumpolar = TRUE;
2391 }
2392 admc = asind(admc);
2393 samc = 90 + admc;
2394 if (samc == 0)
2395 is_invalid = TRUE;
2396 if (fabs(samc) > 0) {
2397 if (mdd >= 0) { /* east */
2398 dfac = (mdd - adp + admc) / samc;
2399 xp[0] = swe_degnorm((dfac - 1) * 90);
2400 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2401 /* eastern object has longer SA than midheaven */
2402 if (dfac > 2 || dfac < 0)
2403 is_invalid = TRUE; /* if this is omitted, funny things happen */
2404 } else {
2405 dfac = (mdd + 180 + adp + admc) / samc;
2406 xp[0] = swe_degnorm((dfac + 1) * 90);
2407 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2408 /* western object has longer SA than midheaven */
2409 if (dfac > 2 || dfac < 0)
2410 is_invalid = TRUE; /* if this is omitted, funny things happen */
2411 }
2412 }
2413 if (is_invalid) {
2414 xp[0] = 0;
2415 hpos = 0;
2416 if (serr != NULL)
2417 strcpy(serr, "Koch house position failed in circumpolar area");
2418 break;
2419 }
2420 if (is_circumpolar) {
2421 if (serr != NULL)
2422 strcpy(serr, "Koch house position, doubtful result in circumpolar area");
2423 }
2424 /* to make sure that a call with a house cusp position returns
2425 * a value within the house, 0.001" is added */
2426 hpos = xp[0] / 30.0 + 1;
2427 break;
2428 case 'C': // Campanus
2429 xeq[0] = swe_degnorm(mdd - 90);
2430 swe_cotrans(xeq, xp, -geolat);
2431 /* to make sure that a call with a house cusp position returns
2432 * a value within the house, 0.001" is added */
2433 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2434 hpos = xp[0] / 30.0 + 1;
2435 break;
2436 case 'U': /* Krusinski-Pisa-Goelzer */
2437 if (fabs(geolat) < VERY_SMALL) { /* code below does not like geolat 0 */
2438 geolat = (geolat >= 0) ? VERY_SMALL : -VERY_SMALL;
2439 }
2440 /* Purpose: find point where planet's house circle (meridian)
2441 * cuts house plane, giving exact planet's house position.
2442 * Input data: ramc, geolat, asc.
2443 */
2444 asc = Asc1(swe_degnorm(armc + 90), geolat, sine, cose);
2445 /* while MC is always south,
2446 * Asc must always be in eastern hemisphere */
2447 asc = fix_asc_polar(asc, armc, eps, geolat);
2448 /*
2449 * Descr: find the house plane 'asc-zenith' - where it intersects
2450 * with equator and at what angle, and then simple find arc
2451 * from asc on that plane to planet's meridian intersection
2452 * with this plane.
2453 */
2454 /* I. find plane of 'asc-zenith' great circle relative to equator:
2455 * solve spherical triangle 'EP-asc-intersection of house circle with equator' */
2456 /* Ia. Find intersection of house plane with equator: */
2457 x[0] = asc; x[1] = 0.0; x[2] = 1.0; /* 1. Start with ascendent on ecliptic */
2458 swe_cotrans(x, x, -eps); /* 2. Transform asc into equatorial coords */
2459 raep = swe_degnorm(armc + 90); /* 3. RA of east point */
2460 x[0] = swe_degnorm(raep - x[0]); /* 4. Rotation - found arc raas-raep */
2461 swe_cotrans(x, x, -(90-geolat)); /* 5. Transform into horizontal coords - arc EP-asc on horizon */
2462 tanx = tand(x[0]);
2463 if (geolat == 0) {
2464 xtemp = (tanx >= 0) ? 90 : -90;
2465 } else {
2466 xtemp = atand(tanx/cosd((90-geolat))); /* 6. Rotation from horizon on circle perpendicular to equator */
2467 }
2468 if (x[0] > 90 && x[0] <= 270)
2469 xtemp = swe_degnorm(xtemp + 180);
2470 x[0] = swe_degnorm(xtemp);
2471 raaz = swe_degnorm(raep - x[0]); /* result: RA of intersection 'asc-zenith' great circle with equator */
2472 /* Ib. Find obliquity to equator of 'asc-zenith' house plane: */
2473 x[0] = raaz; x[1] = 0.0;
2474 x[0] = swe_degnorm(raep - x[0]); /* 1. Rotate start point relative to EP */
2475 swe_cotrans(x, x, -(90-geolat)); /* 2. Transform into horizontal coords */
2476 x[1] = x[1] + 90; /* 3. Add 90 deg do decl - so get the point on house plane most distant from equ. */
2477 swe_cotrans(x, x, 90-geolat); /* 4. Rotate back to equator */
2478 oblaz = x[1]; /* 5. Obliquity of house plane to equator */
2479 /* II. Next find asc and planet position on house plane,
2480 * so to find relative distance of planet from
2481 * coords beginning. */
2482 /* IIa. Asc on house plane relative to intersection
2483 * of equator with 'asc-zenith' plane. */
2484 xasc[0] = asc; xasc[1] = 0.0; xasc[2] = 1.0;
2485 swe_cotrans(xasc, xasc, -eps);
2486 xasc[0] = swe_degnorm(xasc[0] - raaz);
2487 xtemp = atand(tand(xasc[0])/cosd(oblaz));
2488 if (xasc[0] > 90 && xasc[0] <= 270)
2489 xtemp = swe_degnorm(xtemp + 180);
2490 xasc[0] = swe_degnorm(xtemp);
2491 /* IIb. Planet on house plane relative to intersection
2492 * of equator with 'asc-zenith' plane */
2493 xp[0] = swe_degnorm(xeq[0] - raaz); /* Rotate on equator */
2494 xtemp = atand(tand(xp[0])/cosd(oblaz)); /* Find arc on house plane from equator */
2495 if (xp[0] > 90 && xp[0] <= 270)
2496 xtemp = swe_degnorm(xtemp + 180);
2497 xp[0] = swe_degnorm(xtemp);
2498 xp[0] = swe_degnorm(xp[0]-xasc[0]); /* find arc between asc and planet, and get planet house position */
2499 /* IIc. Distance from planet to house plane on declination circle: */
2500 x[0] = xeq[0];
2501 x[1] = xeq[1];
2502 swe_cotrans(x, x, oblaz);
2503 xp[1] = xeq[1] - x[1]; /* How many degrees is the point on declination circle from house circle */
2504 /* to make sure that a call with a house cusp position returns
2505 * a value within the house, 0.001" is added */
2506 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2507 hpos = xp[0] / 30.0 + 1;
2508 break;
2509 case 'H': // horizon / azimuth
2510 xeq[0] = swe_degnorm(mdd - 90);
2511 swe_cotrans(xeq, xp, 90 - geolat);
2512 /* to make sure that a call with a house cusp position returns
2513 * a value within the house, 0.001" is added */
2514 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2515 hpos = xp[0] / 30.0 + 1;
2516 break;
2517 case 'R': // Regiomontanus
2518 if (fabs(mdd) < VERY_SMALL)
2519 xp[0] = 270;
2520 else if (180 - fabs(mdd) < VERY_SMALL)
2521 xp[0] = 90;
2522 else {
2523 if (90 - fabs(geolat) < VERY_SMALL) {
2524 if (geolat > 0)
2525 geolat = 90 - VERY_SMALL;
2526 else
2527 geolat = -90 + VERY_SMALL;
2528 }
2529 if (90 - fabs(de) < VERY_SMALL) {
2530 if (de > 0)
2531 de = 90 - VERY_SMALL;
2532 else
2533 de = -90 + VERY_SMALL;
2534 }
2535 a = tand(geolat) * tand(de) + cosd(mdd);
2536 xp[0] = swe_degnorm(atand(-a / sind(mdd)));
2537 if (mdd < 0)
2538 xp[0] += 180;
2539 xp[0] = swe_degnorm(xp[0]);
2540 /* to make sure that a call with a house cusp position returns
2541 * a value within the house, 0.001" is added */
2542 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2543 }
2544 hpos = xp[0] / 30.0 + 1;
2545 break;
2546 /* with Sunshine and APC houses, the method is the same,
2547 * except that Sunshine users dsun = (declination of Sun),
2548 * whereas APC uses dsun = (declination of ascendant).
2549 * The Sunshine method has more problems within the polar
2550 * circles, if the Sun is circumpolar. The ascendant is never
2551 * circumpolar except if it coincides with the MC or the IC.
2552 */
2553 case 'I': case 'i': // sunshine houses (Makransky)
2554 case 'Y': // APC houses (Knegt)
2555 if (geolat > 90 - MILLIARCSEC)
2556 geolat = 90 - MILLIARCSEC;
2557 if (geolat < -90 + MILLIARCSEC)
2558 geolat = -90 + MILLIARCSEC;
2559 //fprintf(stdout, "in=%f, mdd=%f\n", xpin[0], mdd);
2560 if (90 - fabs(de) < VERY_SMALL) {
2561 if (de > 0)
2562 de = 90 - VERY_SMALL;
2563 else
2564 de = -90 + VERY_SMALL;
2565 }
2566 a = tand(geolat) * tand(de) + cosd(mdd);
2567 xp[0] = swe_degnorm(atand(-a / sind(mdd)));
2568 if (mdd < 0)
2569 xp[0] += 180;
2570 xp[0] = swe_degnorm(xp[0]); // house position with hsys = 'R'
2571 /* is object above horizon? */
2572 sinad = tand(de) * tand(geolat);
2573 a = sinad + cosd(mdd);
2574 if (a >= 0)
2575 is_above_hor = TRUE;
2576 /* height of armc above horizon */
2577 harmc = 90 - geolat;
2578 if (geolat < 0)
2579 harmc = 90 + geolat;
2580 /* meridian distance of crossing of house position line with equator */
2581 darmc = swe_degnorm(xp[0] - 270);
2582 if (darmc > 180) {
2583 is_western_half = TRUE;
2584 darmc = (360 - darmc);
2585 }
2586 /* semi-diurnal arc of sun */
2587 sinad = tand(dsun) * tand(geolat);
2588 if (sinad >= 1)
2589 ad = 90;
2590 //ad = 90 - VERY_SMALL;
2591 else if (sinad <= -1)
2592 ad = -90;
2593 //ad = -(90 - VERY_SMALL);
2594 else
2595 ad = asind(sinad);
2596 sad = 90 + ad;
2597 san = 90 - ad;
2598 //fprintf(stdout, "in=%f, above=%d, sad=%f, san=%f, sinad=%f\n", xpin[0], (int) is_above_hor, sad, san, sinad);
2599 /* circumpolar sun has diurnal arc = 0 and object is above the horizon:
2600 * house position = 10 (270°) */
2601 if (sad == 0 && is_above_hor) {
2602 xp[0] = 270;
2603 /* circumpolar sun has nocturnal arc = 0 and object is below the horizon:
2604 * house position = 4 (90°) */
2605 } else if (san == 0 && !is_above_hor) {
2606 xp[0] = 90;
2607 /* otherwise we can calculate the house position */
2608 } else {
2609 sa = sad;
2610 if (!is_above_hor) {
2611 dsun = -dsun;
2612 sa = san;
2613 darmc = 180 - darmc;
2614 is_western_half = !is_western_half;
2615 }
2616 /* length of position line between south point and equator */
2617 a = acosd(cosd(harmc) * cosd(darmc));
2618 if (a < VERY_SMALL)
2619 a = VERY_SMALL;
2620 /* sine of angle between position line and equator */
2621 sinpsi = sind(harmc) / sind(a);
2622 if (sinpsi > 1) sinpsi = 1;
2623 if (sinpsi < -1) sinpsi = -1;
2624 /* meridian distance of crossing of house position line with solar diurnal arc */
2625 y = sind(dsun) / sinpsi;
2626 if (y > 1)
2627 y = 90 - VERY_SMALL;
2628 else if (y < -1)
2629 y = - (90 - VERY_SMALL);
2630 else
2631 y = asind(y);
2632 d = acosd(cosd(y) / cosd(dsun));
2633 if (dsun < 0) d = -d;
2634 if (geolat < 0) d = -d;
2635 darmc += d;
2636 if (is_western_half)
2637 xp[0] = 270 - (darmc / sa) * 90;
2638 else
2639 xp[0] = 270 + (darmc / sa) * 90;
2640 if (!is_above_hor)
2641 xp[0] = swe_degnorm(xp[0] + 180);
2642 }
2643 /* to make sure that a call with a house cusp position returns
2644 * a value within the house, 0.001" is added */
2645 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2646 hpos = xp[0] / 30.0 + 1;
2647 break;
2648 case 'T': // Polich-Page ("topocentric")
2649 fh = geolat;
2650 if (fh > 89.999)
2651 fh = 89.999;
2652 if (fh < -89.999)
2653 fh = -89.999;
2654 mdd = swe_degnorm(mdd);
2655 if (de > 90 - VERY_SMALL)
2656 de = 90 - VERY_SMALL;
2657 if (de < -90 + VERY_SMALL)
2658 de = -90 + VERY_SMALL;
2659 sinad = tand(de) * tand(fh);
2660 if (sinad > 1.0) sinad = 1.0;
2661 if (sinad < -1.0) sinad = -1.0;
2662 a = sinad + cosd(mdd);
2663 if (a >= 0)
2664 is_above_hor = TRUE;
2665 /* mirror everything below the horizon to the opposite point
2666 * above the horizon */
2667 if (!is_above_hor) {
2668 ra = swe_degnorm(ra + 180);
2669 de = -de;
2670 mdd = swe_degnorm(mdd + 180);
2671 }
2672 /* mirror everything on western hemisphere to eastern hemisphere */
2673 if (mdd > 180) {
2674 ra = swe_degnorm(armc - mdd);
2675 }
2676 /* binary search for "topocentric" position line of body */
2677 tanfi = tand(fh);
2678 ra0 = swe_degnorm(armc + 90);
2679 xp[1] = 1;
2680 xeq[1] = de;
2681 fac = 2;
2682 nloop = 0;
2683 while (fabs(xp[1]) > 0.000001 && nloop < 1000) {
2684 if (xp[1] > 0) {
2685 fh = atand(tand(fh) - tanfi / fac);
2686 ra0 -= 90 / fac;
2687 } else {
2688 fh = atand(tand(fh) + tanfi / fac);
2689 ra0 += 90 / fac;
2690 }
2691 xeq[0] = swe_degnorm(ra - ra0);
2692 swe_cotrans(xeq, xp, 90 - fh);
2693 fac *= 2;
2694 nloop++;
2695 }
2696 hpos = swe_degnorm(ra0 - armc);
2697 /* mirror back to west */
2698 if (mdd > 180)
2699 hpos = swe_degnorm(-hpos);
2700 /* mirror back to below horizon */
2701 if (!is_above_hor)
2702 hpos = swe_degnorm(hpos + 180);
2703 hpos = swe_degnorm(hpos - 90) / 30 + 1;
2704 break;
2705 case 'P': // Placidus
2706 case 'G': // Gauquelin
2707 /* circumpolar region */
2708 if (90 - fabs(de) <= fabs(geolat)) {
2709 if (de * geolat < 0)
2710 xp[0] = swe_degnorm(90 + mdn / 2);
2711 else
2712 xp[0] = swe_degnorm(270 + mdd / 2);
2713 if (serr != NULL)
2714 strcpy(serr, "Otto Ludwig procedure within circumpolar regions.");
2715 } else {
2716 sinad = tand(de) * tand(geolat);
2717 ad = asind(sinad);
2718 a = sinad + cosd(mdd);
2719 if (a >= 0)
2720 is_above_hor = TRUE;
2721 sad = 90 + ad;
2722 san = 90 - ad;
2723 if (is_above_hor)
2724 xp[0] = (mdd / sad + 3) * 90;
2725 else
2726 xp[0] = (mdn / san + 1) * 90;
2727 /* to make sure that a call with a house cusp position returns
2728 * a value within the house, 0.001" is added */
2729 xp[0] = swe_degnorm(xp[0] + MILLIARCSEC);
2730 }
2731 if (hsys == 'G') {
2732 xp[0] = 360 - xp[0]; /* Gauquelin sectors are in clockwise direction */
2733 hpos = xp[0] / 10.0 + 1;
2734 } else {
2735 hpos = xp[0] / 30.0 + 1;
2736 }
2737 break;
2738 default:
2739 hpos = 0;
2740 if (swe_houses_armc_ex2(armc, geolat, eps, hsys, hcusp, ascmc, NULL, NULL, serr) == ERR) {
2741 if (serr != NULL)
2742 sprintf(serr, "swe_house_pos(): failed for system %c", hsys);
2743 break;
2744 }
2745 if (swe_difdeg2n(hcusp[6], hcusp[1]) > 0) {
2746 d = swe_degnorm(xpin[0] - hcusp[1]);
2747 for (i = 1; i <= 12; i++) {
2748 j = i + 1;
2749 if (j > 12)
2750 c2 = 360;
2751 else
2752 c2 = swe_degnorm(hcusp[j] - hcusp[1]);
2753 if (d < c2) break;
2754 }
2755 c1 = swe_degnorm(hcusp[i] - hcusp[1]);
2756 } else { // houses retrograde
2757 d = swe_degnorm(hcusp[1] - xpin[0]);
2758 for (i = 1; i <= 12; i++) {
2759 j = i + 1;
2760 if (j > 12)
2761 c2 = 360;
2762 else
2763 c2 = swe_degnorm(hcusp[1] - hcusp[j]);
2764 if (d < c2) break;
2765 }
2766 c1 = swe_degnorm(hcusp[1] - hcusp[i]);
2767 }
2768 hsize = c2 - c1;
2769 if (hsize == 0) {
2770 hpos = i;
2771 } else {
2772 hpos = i + (d - c1) / hsize;
2773 }
2774 if (serr != NULL)
2775 sprintf(serr, "swe_house_pos(): using simplified algorithm for system %c\n", hsys);
2776 break;
2777 }
2778 return hpos;
2779 }
2780
sunshine_init(double lat,double dec,double xh[])2781 static int sunshine_init(double lat, double dec, double xh[])
2782 {
2783 double ad, nsa, dsa, arg;
2784 // ascensional difference: sin ad = tan dec tan lat
2785 // or near +- 90 if Sun circumpolar
2786 arg = tand(dec) * tand(lat);
2787 if (arg >= 1) {
2788 ad = 90 - VERY_SMALL;
2789 } else if (arg <= -1) {
2790 ad = -90 + VERY_SMALL;
2791 } else {
2792 ad = asind(arg);
2793 }
2794 nsa = 90 - ad;
2795 dsa = 90 + ad;
2796 xh[2] = -2 * nsa / 3;
2797 xh[3] = -1 * nsa / 3;
2798 xh[5] = 1 * nsa / 3;
2799 xh[6] = 2 * nsa / 3;
2800 xh[8] = -2 * dsa / 3;
2801 xh[9] = -1 * dsa / 3;
2802 xh[11] = 1 * dsa / 3;
2803 xh[12] = 2 * dsa / 3;
2804 if (fabs(arg) >= 1)
2805 return ERR;
2806 return OK;
2807 }
2808
sunshine_solution_makransky(double ramc,double lat,double ecl,struct houses * hsp)2809 static int sunshine_solution_makransky(double ramc, double lat, double ecl, struct houses *hsp)
2810 {
2811 double xh[13];
2812 double md;
2813 double zd; // zenith distance of house circle, along prime vertical
2814 double pole, q, w, a, b, c, f, cu, r = 0, rah;
2815 double sinlat, coslat, tanlat, tandec, sinecl;
2816 double dec = hsp->sundec;
2817 sinlat = sind(lat);
2818 coslat = cosd(lat);
2819 tanlat = tand(lat);
2820 tandec = tand(dec);
2821 sinecl = sind(ecl);
2822 int ih;
2823 // if (90 - fabs(lat) <= ecl) {
2824 // strcpy(hsp->serr, "Sunshine in polar circle not allowed");
2825 // return ERR;
2826 // }
2827 if (sunshine_init(lat, dec, xh) == ERR)
2828 return ERR;
2829 for (ih = 1; ih <= 12; ih++) {
2830 double z = 0;
2831 if ((ih - 1) % 3 == 0) continue; // skip 1,4,7,10
2832 md = fabs(xh[ih]);
2833 if (ih <= 6)
2834 rah = swe_degnorm(ramc + 180 + xh[ih]);
2835 else
2836 rah = swe_degnorm(ramc + xh[ih]);
2837 if (lat < 0) { // Makransky deals with southern latitude this way
2838 rah = swe_degnorm(180 + rah);
2839 }
2840 // HP is the house point on the semidiurnal arc
2841 // CP = intersection house meridian with prime vertical
2842 // MP = intersection house meridian with equator
2843 // XP = intersection house circle with prime meridian
2844 if (md == 90) {
2845 // CP = east point (or west point),
2846 // HP is on meridian east point - north pole
2847 // use triangle CP - HP - XP with long side dec
2848 // and angle 90 - lat.
2849 // use tan b = cos alph tan c = sin lat tan dec
2850 zd = 90.0 - atand(sinlat * tandec);
2851 } else {
2852 if (md < 90) {
2853 // triangle 1) CP, Zenith, north pole: side 90-lat, angle md at
2854 // north pole.
2855 // tan a = cos lat * tan md
2856 // a is distance of CP from zenith on prime vertical
2857 a = atand(coslat * tand(md));
2858 } else {
2859 // triangle 1) MP - east point - CP : side b = 90-md, angle lat
2860 // at east point
2861 // tan c = tan md / cos lat
2862 // a is distance of CP from zenith on prime vertical
2863 a = atand(tand(md - 90) / coslat); // lat = 90 not allowed
2864 }
2865 // triangle 2) CP, MP, east point: side 90 - md, angle lat
2866 // tan b = tan lat * cos md
2867 // b is distance of CP from equator
2868 b = atand(tanlat * cosd(md));
2869 // c is distance of HP house point from CP, along its meridian.
2870 if (ih <= 6)
2871 c = b + dec;
2872 else
2873 c = b - dec;
2874 // triangle 3) HP - CP - XP, side c, angle from triangle 1)
2875 // tan f = sin lat * sin md * tan c;
2876 // f is the distance from CP to XP
2877 f = atand(sinlat * sind(md) * tand(c));
2878 // a + f give zd, the zenith distance
2879 // of house circle measured on prime vertical.
2880 zd = a + f;
2881 }
2882 pole = asind(sind(zd) * sinlat);
2883 q = asind(tandec * tand(pole));
2884 if (ih <= 3 || ih >= 11)
2885 w = swe_degnorm(rah - q);
2886 else
2887 w = swe_degnorm(rah + q);
2888 if (w == 90) {
2889 r = atand(sind(ecl) * tand(pole));
2890 if (ih <= 3 || ih >= 11)
2891 cu = 90 + r;
2892 else
2893 cu = 90 - r;
2894 } else if (w == 270) {
2895 r = atand(sinecl * tand(pole));
2896 if (ih <= 3 || ih >= 11)
2897 cu = 270 - r;
2898 else
2899 cu = 270 + r;
2900 } else {
2901 double m;
2902 m = atand(fabs(tand(pole) / cosd(w)));
2903 if (ih <= 3 || ih >= 11) {
2904 if (w > 90 && w < 270)
2905 z = m - ecl;
2906 else
2907 z = m + ecl;
2908 } else {
2909 if (w > 90 && w < 270)
2910 z = m + ecl;
2911 else
2912 z = m - ecl;
2913 }
2914 if (z == 90) {
2915 if (w < 180)
2916 cu = 90;
2917 else
2918 cu = 270;
2919 } else {
2920 // r is between 0 and 90
2921 r = atand(fabs(cosd(m) * tand(w) / cosd(z)));
2922 if (w < 90)
2923 cu = r;
2924 else if (w > 90 && w < 180)
2925 cu = 180 - r;
2926 else if (w > 180 && w < 270)
2927 cu = 180 + r;
2928 else
2929 cu = 360 - r;
2930 }
2931 if (z > 90) {
2932 // i am not sure if I understood the remark 'value will fall away from cancer..
2933 // on page 146 correctly.
2934 if (w < 90)
2935 cu = 180 - r;
2936 else if (w > 90 && w < 180)
2937 cu = + r;
2938 else if (w > 180 && w < 270)
2939 cu = 360 - r;
2940 else
2941 cu = 180 + r;
2942 }
2943 if (lat < 0) // Makransky deals with southern latitude this way
2944 cu = swe_degnorm(cu + 180);
2945 }
2946 hsp->cusp[ih] = cu;
2947 }
2948 return OK;
2949 }
2950
sunshine_solution_treindl(double ramc,double lat,double ecl,struct houses * hsp)2951 static int sunshine_solution_treindl(double ramc, double lat, double ecl, struct houses *hsp)
2952 {
2953 double xh[13];
2954 double mcdec, sinlat, coslat, cosdec, tandec, sinecl, cosecl;
2955 double xhs, pole, a, cosa, alph, alpha2, c, cosc, b, sinzd, zd, rax, hc;
2956 int ih, retval = OK;
2957 AS_BOOL mc_under_horizon;
2958 double dec = hsp->sundec;
2959 // if (90 - fabs(lat) <= ecl) {
2960 // strcpy(hsp->serr, "Sunshine in polar circle not allowed");
2961 // return ERR;
2962 // }
2963 sinlat = sind(lat);
2964 coslat = cosd(lat);
2965 cosdec = cosd(dec);
2966 tandec = tand(dec);
2967 sinecl = sind(ecl);
2968 cosecl = cosd(ecl);
2969 sunshine_init(lat, dec, xh);
2970 // find out if MC under horizon
2971 mcdec = atand(sind(ramc) * tand(ecl));
2972 mc_under_horizon = fabs(lat - mcdec) > 90;
2973 if (mc_under_horizon && SUNSHINE_KEEP_MC_SOUTH) {
2974 // we have switched ac/mc, invert offsets on diurnal arcs
2975 for (ih = 2; ih <= 12; ih++) {
2976 xh[ih] = -xh[ih];
2977 }
2978 }
2979 //if (sunshine_init(lat, dec, xh) == ERR)
2980 // return ERR;
2981 // HP is the house point on the semidiurnal arc
2982 // CP = intersection house meridian with prime vertical
2983 // MP = intersection house meridian with equator
2984 // XP = intersection house circle with prime vertical
2985 // EP = intersection house circle with equator
2986 // MP0 = intersection semiarc with meridian
2987 for (ih = 1; ih <= 12; ih++) {
2988 if ((ih - 1) % 3 == 0) continue; // skip 1,4,7,10
2989 xhs = 2 * asind(cosdec * sind(xh[ih] / 2)); // x' great-circle length of x
2990 // compute triangle north pole - mp0 -hp
2991 // we have two sides 90 - dec, base xhs, ange at pole x
2992 // derive from cosine rule
2993 cosa = tandec * tand(xhs / 2);
2994 alph = acosd(cosa);
2995 // compute triangle south point - mp0 - hp
2996 // we have: side x', side b = 90 - lat + dec, angle alpha2 between the sides.
2997 // we want: angle zd at south point.
2998 // we compute first the other side with Seitencosinus-Satz
2999 // cos c = cos x' cos b + sin x' sin b cos alpha2
3000 // for nocturnal side: use alpha, side b = 90 - lat - dec;
3001 // zd will be angle at north point.
3002 if (ih > 7) {
3003 // complementary angle
3004 alpha2 = 180 - alph;
3005 b = 90 - lat + dec;
3006 } else { // nocturnal side
3007 alpha2 = alph;
3008 b = 90 - lat - dec;
3009 }
3010 // b can be zero, xhs can 90, c can get small.
3011 cosc = cosd(xhs) * cosd(b) + sind(xhs) * sind(b) * cosd(alpha2);
3012 c = acosd(cosc);
3013 // now Sinussatz
3014 if (c < 1e-6) {
3015 sprintf(hsp->serr, "Sunshine house %d c=%le very small", ih, c);
3016 retval = ERR;
3017 }
3018 sinzd = sind(xhs) * sind(alpha2) / sind(c);
3019 zd = asind(sinzd);
3020 // compute intersection house circle with equator, point rax
3021 // day side triangle south point - meridian point - EP
3022 // night side: triangle north point
3023 // sides: 90 - lat, angle zd
3024 rax = atand(coslat * tand(zd));
3025 // compute pole height (distance of house circle pole from equator
3026 // with triangle at west point
3027 pole = asind(sinzd * sinlat);
3028 if (ih <= 6) {
3029 pole = -pole;
3030 a = swe_degnorm(rax + ramc + 180);
3031 } else {
3032 a = swe_degnorm(ramc + rax);
3033 }
3034 // with pole and a = rectascension of equator intersection, we use Asc1()
3035 // like with many other house systems, to intersect house circle with eclitpic
3036 hc = Asc1(a, pole, sinecl, cosecl);
3037 hsp->cusp[ih] = hc;
3038 }
3039 if (mc_under_horizon && ! SUNSHINE_KEEP_MC_SOUTH) {
3040 for (ih = 2; ih <= 12; ih++) {
3041 if ((ih - 1) % 3 == 0) continue; // skip 1,4,7,10
3042 hsp->cusp[ih] = swe_degnorm(hsp->cusp[ih] + 180);
3043 }
3044 }
3045 return retval;
3046 }
3047