diff --git a/deps/swisseph/swedll.h b/deps/swisseph/swedll.h index 296a777..0d74c4e 100644 --- a/deps/swisseph/swedll.h +++ b/deps/swisseph/swedll.h @@ -153,10 +153,18 @@ DllImport int CALL_CONV_IMP swe_houses_ex( double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *hcusps, double *ascmc); +DllImport int CALL_CONV_IMP swe_houses_ex2( + double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, + double *hcusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr); + DllImport int CALL_CONV_IMP swe_houses_armc( double armc, double geolat, double eps, int hsys, double *hcusps, double *ascmc); +DllImport int CALL_CONV_IMP swe_houses_armc_ex2( + double armc, double geolat, double eps, int hsys, + double *hcusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr); + DllImport double CALL_CONV_IMP swe_house_pos( double armc, double geolon, double eps, int hsys, double *xpin, char *serr); diff --git a/deps/swisseph/swehouse.c b/deps/swisseph/swehouse.c index cd2c62c..8237e52 100644 --- a/deps/swisseph/swehouse.c +++ b/deps/swisseph/swehouse.c @@ -66,8 +66,11 @@ house and (simple) aspect calculation #include #define MILLIARCSEC (1.0 / 3600000.0) +#define SOLAR_YEAR 365.24219893 +#define ARMCS ((SOLAR_YEAR+1) / SOLAR_YEAR * 360) static double Asc1(double, double, double, double); +static double AscDash(double, double, double, double); static double Asc2(double, double, double, double); static int CalcH(double th, double fi, double ekl, char hsy, struct houses *hsp); static int sidereal_houses_ecl_t0(double tjde, @@ -77,7 +80,10 @@ static int sidereal_houses_ecl_t0(double tjde, double lat, int hsys, double *cusp, - double *ascmc); + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr); static int sidereal_houses_trad(double tjde, int32 iflag, double armc, @@ -86,7 +92,10 @@ static int sidereal_houses_trad(double tjde, double lat, int hsys, double *cusp, - double *ascmc); + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr); static int sidereal_houses_ssypl(double tjde, double armc, double eps, @@ -94,7 +103,10 @@ static int sidereal_houses_ssypl(double tjde, double lat, int hsys, double *cusp, - double *ascmc); + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr); static int sunshine_solution_makransky(double ramc, double lat, double ecl, struct houses *hsp); static int sunshine_solution_treindl(double ramc, double lat, double ecl, struct houses *hsp); #if 0 @@ -136,7 +148,7 @@ int CALL_CONV swe_houses(double tjd_ut, int result = swe_calc_ut(tjd_ut, SE_SUN, flags, xp, NULL); if (result < 0) { // in case of failure, Porphyry houses - result = swe_houses_armc(armc, geolat, eps + nutlo[1], 'O', cusp, ascmc); + result = swe_houses_armc_ex2(armc, geolat, eps + nutlo[1], 'O', cusp, ascmc, NULL, NULL, NULL); return ERR; } ascmc[9] = xp[1]; // declination in ascmc[9]; @@ -158,31 +170,50 @@ int CALL_CONV swe_houses(double tjd_ut, } } #endif - retc = swe_houses_armc(armc, geolat, eps + nutlo[1], hsys, cusp, ascmc); + retc = swe_houses_armc_ex2(armc, geolat, eps + nutlo[1], hsys, cusp, ascmc, NULL, NULL, NULL); return retc; } -/* housasp.c +// For explanation see function swe_houses_ex2() below. +int CALL_CONV swe_houses_ex(double tjd_ut, + int32 iflag, + double geolat, + double geolon, + int hsys, + double *cusp, + double *ascmc) +{ + return swe_houses_ex2(tjd_ut, iflag, geolat, geolon, hsys, cusp, ascmc, NULL, NULL, NULL); +} + +/* + * Function returns OK or ERR. * cusps are returned in double cusp[13], * or cusp[37] with house system 'G'. - * cusp[1...12] houses 1 - 12 - * additional points are returned in ascmc[10]. - * ascmc[0] = ascendant - * ascmc[1] = mc - * ascmc[2] = armc - * ascmc[3] = vertex - * ascmc[4] = equasc * "equatorial ascendant" * - * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) * - * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) * - * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) * + * cusp[1...12] houses 1 - 12 + * ascmc[0...10] additional points: + * ascmc[0] = ascendant + * ascmc[1] = mc + * ascmc[2] = armc + * ascmc[3] = vertex + * ascmc[4] = equasc * "equatorial ascendant" * + * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) * + * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) * + * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) * + * cusp_speed[1...12] speeds (daily motions) of the cusps. + * ascmc_speed[0...10] speeds (daily motions) of the additional points. + * serr error message or warning */ -int CALL_CONV swe_houses_ex(double tjd_ut, +int CALL_CONV swe_houses_ex2(double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusp, - double *ascmc) + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr) { int i, retc = 0; double armc, eps_mean, nutlo[2]; @@ -237,13 +268,13 @@ int CALL_CONV swe_houses_ex(double tjd_ut, } if (iflag & SEFLG_SIDEREAL) { if (sip->sid_mode & SE_SIDBIT_ECL_T0) - retc = sidereal_houses_ecl_t0(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc); + retc = sidereal_houses_ecl_t0(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); else if (sip->sid_mode & SE_SIDBIT_SSY_PLANE) - retc = sidereal_houses_ssypl(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc); + retc = sidereal_houses_ssypl(tjde, armc, eps_mean + nutlo[1], nutlo, geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); else - retc = sidereal_houses_trad(tjde, iflag, armc, eps_mean + nutlo[1], nutlo[0], geolat, hsys, cusp, ascmc); + retc = sidereal_houses_trad(tjde, iflag, armc, eps_mean + nutlo[1], nutlo[0], geolat, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); } else { - retc = swe_houses_armc(armc, geolat, eps_mean + nutlo[1], hsys, cusp, ascmc); + retc = swe_houses_armc_ex2(armc, geolat, eps_mean + nutlo[1], hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); if (toupper(hsys) == 'I') ascmc[9] = xp[1]; // declination in ascmc[9]; } @@ -291,7 +322,10 @@ static int sidereal_houses_ecl_t0(double tjde, double lat, int hsys, double *cusp, - double *ascmc) + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr) { int i, j, retc = OK; double x[6], xvpx[6], x2[6], epst0, xnorm[6]; @@ -347,7 +381,7 @@ static int sidereal_houses_ecl_t0(double tjde, /* auxiliary armc */ armcx = swe_degnorm(armc - dvpx); /* 3 */ /* compute axes and houses: */ - retc = swe_houses_armc(armcx, lat, epsx, hsys, cusp, ascmc); /* 4 */ + retc = swe_houses_armc_ex2(armcx, lat, epsx, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); /* 4 */ /* distance between auxiliary vernal point and * vernal point of t0 (a section on the sidereal plane) */ dvpxe = acos(swi_dot_prod_unit(x, xvpx)) * RADTODEG; /* 5 */ @@ -395,7 +429,10 @@ static int sidereal_houses_ssypl(double tjde, double lat, int hsys, double *cusp, - double *ascmc) + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr) { int i, j, retc = OK; double x[6], x0[6], xvpx[6], x2[6], xnorm[6]; @@ -454,7 +491,7 @@ static int sidereal_houses_ssypl(double tjde, /* auxiliary armc */ armcx = swe_degnorm(armc - dvpx); /* 3 */ /* compute axes and houses: */ - retc = swe_houses_armc(armcx, lat, epsx, hsys, cusp, ascmc); /* 4 */ + retc = swe_houses_armc_ex2(armcx, lat, epsx, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr); /* 4 */ /* distance between the auxiliary vernal point at t and * the sidereal zero point of 2000 at t * (a section on the sidereal plane). @@ -503,7 +540,10 @@ static int sidereal_houses_trad(double tjde, double lat, int hsys, double *cusp, - double *ascmc) + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr) { int i, retc = OK; double ay; @@ -523,7 +563,7 @@ static int sidereal_houses_trad(double tjde, ihs2 = 'E'; //fprintf(stderr, "armc=%f\n", armc); //if (hsys == 'P') fprintf(stderr, "ay=%f, t=%f %c", ay, tjde, (char) hsys); - retc = swe_houses_armc(armc, lat, eps, ihs2, cusp, ascmc); + retc = swe_houses_armc_ex2(armc, lat, eps, ihs2, cusp, ascmc, cusp_speed, ascmc_speed, serr); //if (hsys == 'P') fprintf(stderr, " h1=%f", cusp[1]); for (i = 1; i <= ito; i++) { //cusp[i] = swe_degnorm(cusp[i] - ay - nutl); @@ -546,33 +586,52 @@ static int sidereal_houses_trad(double tjde, return retc; } +// For explanation see function swe_houses_armc_ex2() below. +int CALL_CONV swe_houses_armc( + double armc, + double geolat, + double eps, + int hsys, + double *cusp, + double *ascmc) +{ + return swe_houses_armc_ex2(armc, geolat, eps, hsys, cusp, ascmc, NULL, NULL, NULL); +} + /* + * Function returns OK or ERR. * this function is required for very special computations * where no date is given for house calculation, * e.g. for composite charts or progressive charts. * cusps are returned in double cusp[13], * or cusp[37] with house system 'G'. - * cusp[1...12] houses 1 - 12 - * additional points are returned in ascmc[10]. - * ascmc[0] = ascendant - * ascmc[1] = mc - * ascmc[2] = armc - * ascmc[3] = vertex - * ascmc[4] = equasc * "equatorial ascendant" * - * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) * - * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) * - * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) * + * cusp[1...12] houses 1 - 12 + * ascmc[0...10] additional points: + * ascmc[0] = ascendant + * ascmc[1] = mc + * ascmc[2] = armc + * ascmc[3] = vertex + * ascmc[4] = equasc * "equatorial ascendant" * + * ascmc[5] = coasc1 * "co-ascendant" (W. Koch) * + * ascmc[6] = coasc2 * "co-ascendant" (M. Munkasey) * + * ascmc[7] = polasc * "polar ascendant" (M. Munkasey) * + * cusp_speed[1...12] speeds (daily motions) of the cusps. + * ascmc_speed[0...10] speeds (daily motions) of the additional points. + * serr error message or warning */ -int CALL_CONV swe_houses_armc( +int CALL_CONV swe_houses_armc_ex2( double armc, double geolat, double eps, int hsys, double *cusp, - double *ascmc) + double *ascmc, + double *cusp_speed, + double *ascmc_speed, + char *serr) { - struct houses h; - int i, retc = 0; + struct houses h, hm1, hp1; + int i, retc = 0, rm1, rp1; int ito; static double saved_sundec = 99; if (toupper(hsys) == 'G') @@ -580,6 +639,12 @@ int CALL_CONV swe_houses_armc( else ito = 12; armc = swe_degnorm(armc); + h.do_speed = FALSE; + h.do_hspeed = FALSE; + if (ascmc_speed != NULL || cusp_speed != NULL) + h.do_speed = TRUE; // is needed if cusp_speed wanted + if (cusp_speed != NULL) + h.do_hspeed = TRUE; if (toupper(hsys) == 'I') { // declination for sunshine houses if (ascmc[9] == 99) { h.sundec = 0; @@ -592,10 +657,13 @@ int CALL_CONV swe_houses_armc( retc = CalcH(armc, geolat, eps, (char)hsys, &h); cusp[0] = 0; // on failure, we only have 12 Porphyry cusps - if (retc < 0) + if (retc < 0) { ito = 12; + if (serr != NULL) strcpy(serr, h.serr); + } for (i = 1; i <= ito; i++) { cusp[i] = h.cusp[i]; + if (h.do_hspeed) cusp_speed[i] = h.cusp_speed[i]; } ascmc[0] = h.ac; /* Asc */ ascmc[1] = h.mc; /* Mid */ @@ -609,29 +677,74 @@ int CALL_CONV swe_houses_armc( ascmc[i] = 0; if (toupper(hsys) == 'I') // declination for sunshine houses ascmc[9] = h.sundec ; + if (h.do_speed && ascmc_speed != NULL) { + ascmc_speed[0] = h.ac_speed; /* Asc */ + ascmc_speed[1] = h.mc_speed; /* Mid */ + ascmc_speed[2] = h.armc_speed; + ascmc_speed[3] = h.vertex_speed; + ascmc_speed[4] = h.equasc_speed; + ascmc_speed[5] = h.coasc1_speed; /* "co-ascendant" (W. Koch) */ + ascmc_speed[6] = h.coasc2_speed; /* "co-ascendant" (M. Munkasey) */ + ascmc_speed[7] = h.polasc_speed; /* "polar ascendant" (M. Munkasey) */ + for (i = SE_NASCMC; i < 10; i++) + ascmc_speed[i] = 0; + } + if (h.do_interpol) { // must compute cusp_speed via interpolation + double dt = 1.0 / 86400; + double darmc = dt * ARMCS; + hm1.do_speed = FALSE; + hm1.do_hspeed = FALSE; + hp1.do_speed = FALSE; + hp1.do_hspeed = FALSE; + if (toupper(hsys) == 'I') { + hm1.sundec = h.sundec; + hp1.sundec = h.sundec; + } + rm1 = CalcH(armc - darmc, geolat, eps, (char)hsys, &hm1); + rp1 = CalcH(armc + darmc, geolat, eps, (char)hsys, &hp1); + if (rp1 >= 0 && rm1 >=0) { + if (fabs(swe_difdeg2n(hp1.ac, h.ac)) > 90) { + hp1 = h; // use only upper interval + dt = dt / 2; + } else if (fabs(swe_difdeg2n(hm1.ac, h.ac)) > 90) { + hm1 = h; // use only lower interval + dt = dt / 2; + } + for (i = 1; i <= 12; i++) { + double dx = swe_difdeg2n(hp1.cusp[i], hm1.cusp[i]); + cusp_speed[i] = dx / 2 / dt ; + } + } + } #ifdef TRACE swi_open_trace(NULL); if (swi_trace_count <= TRACE_COUNT_MAX) { if (swi_fp_trace_c != NULL) { - fputs("\n/*SWE_HOUSES_ARMC*/\n", swi_fp_trace_c); + fputs("\n/*SWE_HOUSES_ARMC_EX2*/\n", swi_fp_trace_c); fprintf(swi_fp_trace_c, " armc = %.9f;", armc); fprintf(swi_fp_trace_c, " geolat = %.9f;", geolat); fprintf(swi_fp_trace_c, " eps = %.9f;", eps); fprintf(swi_fp_trace_c, " hsys = %d;\n", hsys); - fprintf(swi_fp_trace_c, " retc = swe_houses_armc(armc, geolat, eps, hsys, cusp, ascmc);\n"); - fputs(" printf(\"swe_houses_armc: %f\\t%f\\t%f\\t%c\\t\\n\", ", swi_fp_trace_c); + fprintf(swi_fp_trace_c, " retc = swe_houses_armc_ex2(armc, geolat, eps, hsys, cusp, ascmc, cusp_speed, ascmc_speed, serr);\n"); + fputs(" printf(\"swe_houses_armc_ex2: %f\\t%f\\t%f\\t%c\\t\\n\", ", swi_fp_trace_c); fputs(" armc, geolat, eps, hsys);\n", swi_fp_trace_c); fputs(" printf(\"retc = %d\\n\", retc);\n", swi_fp_trace_c); fputs(" printf(\"cusp:\\n\");\n", swi_fp_trace_c); - fputs(" for (i = 0; i < 12; i++)\n", swi_fp_trace_c); + fputs(" for (i = 1; i <= 12; i++)\n", swi_fp_trace_c); fputs(" printf(\" %d\\t%f\\n\", i, cusp[i]);\n", swi_fp_trace_c); fputs(" printf(\"ascmc:\\n\");\n", swi_fp_trace_c); fputs(" for (i = 0; i < 10; i++)\n", swi_fp_trace_c); fputs(" printf(\" %d\\t%f\\n\", i, ascmc[i]);\n", swi_fp_trace_c); + fputs(" printf(\"cusp_speed:\\n\");\n", swi_fp_trace_c); + fputs(" for (i = 1; i <= 12; i++)\n", swi_fp_trace_c); + fputs(" printf(\" %d\\t%f\\n\", i, cusp_speed[i]);\n", swi_fp_trace_c); + fputs(" printf(\"ascmc_speed:\\n\");\n", swi_fp_trace_c); + fputs(" for (i = 0; i < 10; i++)\n", swi_fp_trace_c); + fputs(" printf(\" %d\\t%f\\n\", i, ascmc_speed[i]);\n", swi_fp_trace_c); fflush(swi_fp_trace_c); } if (swi_fp_trace_out != NULL) { - fprintf(swi_fp_trace_out, "swe_houses_armc: %f\t%f\t%f\t%c\t\n", armc, geolat, eps, hsys); + fprintf(swi_fp_trace_out, "swe_houses_armc_ex2: %f\t%f\t%f\t%c\t\n", armc, geolat, eps, hsys); fprintf(swi_fp_trace_out, "retc = %d\n", retc); fputs("cusp:\n", swi_fp_trace_out); for (i = 1; i <= 12; i++) @@ -820,6 +933,7 @@ static int CalcH( int niter_max = 100; // maximum iterations allowed with Placidus double cuspsv; *hsp->serr = '\0'; + hsp->do_interpol = 0; cose = cosd(ekl); sine = sind(ekl); tane = tand(ekl); @@ -845,10 +959,23 @@ static int CalcH( hsp->mc = 270; } /* if */ hsp->mc = swe_degnorm(hsp->mc); + if (hsp->do_speed) hsp->mc_speed = AscDash(th, 0, sine, cose); /* ascendant */ hsp->ac = Asc1(th + 90, fi, sine, cose); + if (hsp->do_speed) + hsp->ac_speed = AscDash(th + 90, fi, sine, cose); + if (hsp->do_hspeed) { + for (i = 0; i <= 12; i++) + hsp->cusp_speed[i] = 0; + } + hsp->armc_speed = ARMCS; + // these cusp[1] and cusp[10] values may be changed further down for some house systems hsp->cusp[1] = hsp->ac; hsp->cusp[10] = hsp->mc; + if (hsp->do_hspeed) { + hsp->cusp_speed[1] = hsp->ac_speed; + hsp->cusp_speed[10] = hsp->mc_speed; + } /* we respect smaller case letter for i, otherwise they are deprecated */ if (hsy > 95 && hsy != 'i') { sprintf(hsp->serr, "use of lower case letters like %c for house systems is deprecated", hsy); @@ -863,8 +990,14 @@ static int CalcH( hsp->ac = swe_degnorm(hsp->ac + 180); hsp->cusp[1] = hsp->ac; } - for (i = 2; i <=12; i++) + for (i = 2; i <=12; i++) { hsp->cusp[i] = swe_degnorm(hsp->cusp[1] + (i-1) * 30); + } + if (hsp->do_hspeed) { + for (i = 1; i <=12; i++) { + hsp->cusp_speed[i] = hsp->ac_speed; + } + } break; case 'D': /* equal, begin at MC */ acmc = swe_difdeg2n(hsp->ac, hsp->mc); @@ -877,6 +1010,11 @@ static int CalcH( hsp->cusp[i] = swe_degnorm(hsp->cusp[10] + (i-10) * 30); for (i = 1; i <= 9; i++) hsp->cusp[i] = swe_degnorm(hsp->cusp[10] + (i + 2) * 30); + if (hsp->do_hspeed) { + for (i = 1; i <=12; i++) { + hsp->cusp_speed[i] = hsp->mc_speed; + } + } break; case 'C': /* Campanus houses and Horizon or Azimut system */ case 'H': @@ -912,6 +1050,14 @@ static int CalcH( hsp->cusp[1] = Asc1(th + 90, fi, sine, cose); hsp->cusp[2] = Asc1(th + 90 + xh2, fh2, sine, cose); hsp->cusp[3] = Asc1(th + 90 + xh1, fh1, sine, cose); + if (hsp->do_hspeed) { + hsp->cusp_speed[11] = AscDash(th + 90 - xh1, fh1, sine, cose); + hsp->cusp_speed[12] = AscDash(th + 90 - xh2, fh2, sine, cose); + if (hsy == 'H') + hsp->cusp_speed[1] = AscDash(th + 90, fi, sine, cose); + hsp->cusp_speed[2] = AscDash(th + 90 + xh2, fh2, sine, cose); + hsp->cusp_speed[3] = AscDash(th + 90 + xh1, fh1, sine, cose); + } /* within polar circle, when mc sinks below horizon and * ascendant changes to western hemisphere, all cusps * must be added 180 degrees. @@ -968,6 +1114,7 @@ static int CalcH( hsy = 'O'; goto porphyry; } + hsp->do_interpol = hsp->do_hspeed; break; case 'K': /* Koch houses */ if (fabs(fi) >= 90 - ekl) { /* within polar circle */ @@ -985,6 +1132,12 @@ static int CalcH( hsp->cusp[12] = Asc1(th + 60 - ad3, fi, sine, cose); hsp->cusp[2] = Asc1(th + 120 + ad3, fi, sine, cose); hsp->cusp[3] = Asc1(th + 150 + 2 * ad3, fi, sine, cose); + if (hsp->do_hspeed) { + hsp->cusp_speed[11] = AscDash(th + 30 - 2 * ad3, fi, sine, cose); + hsp->cusp_speed[12] = AscDash(th + 60 - ad3, fi, sine, cose); + hsp->cusp_speed[2] = AscDash(th + 120 + ad3, fi, sine, cose); + hsp->cusp_speed[3] = AscDash(th + 150 + 2 * ad3, fi, sine, cose); + } break; case 'L': /* Pullen SD sinusoidal delta, ex Neo-Porphyry */ { @@ -1012,6 +1165,7 @@ static int CalcH( hsp->cusp[3] = swe_degnorm(hsp->ac + 60 + 3 * d); } } + hsp->do_interpol = hsp->do_hspeed; break; case 'N': /* whole signs, begin at 0° Aries */ acmc = swe_difdeg2n(hsp->ac, hsp->mc); @@ -1037,6 +1191,16 @@ static int CalcH( hsp->cusp[3] = swe_degnorm(hsp->ac + (180 - acmc) / 3 * 2); hsp->cusp[11] = swe_degnorm(hsp->mc + acmc / 3); hsp->cusp[12] = swe_degnorm(hsp->mc + acmc / 3 * 2); + if (hsp->do_hspeed) { + double q1_speed = hsp->ac_speed - hsp->mc_speed; // rate of growth of quadrant 1 + // double q4_speed = hsp->mc_speed - hsp->ac_speed; // rate of growth of quadrant 4 + hsp->cusp_speed[1] = hsp->ac_speed; // may have been destroyed if defaulting from Gauquelin + hsp->cusp_speed[10] = hsp->mc_speed; // dito + hsp->cusp_speed[2] = hsp->ac_speed - q1_speed / 3; + hsp->cusp_speed[3] = hsp->ac_speed - q1_speed / 3 * 2; + hsp->cusp_speed[11] = hsp->ac_speed + q1_speed / 3; + hsp->cusp_speed[12] = hsp->ac_speed + q1_speed / 3 * 2; + } break; case 'Q': /* Pullen sinusoidal ratio */ { @@ -1081,6 +1245,7 @@ static int CalcH( hsp->cusp[3] = swe_degnorm(hsp->cusp[2] + xr4); // house 2 size xr^4 } } + hsp->do_interpol = hsp->do_hspeed; break; case 'R': /* Regiomontanus houses */ fh1 = atand (tanfi * 0.5); @@ -1089,6 +1254,12 @@ static int CalcH( hsp->cusp[12] = Asc1(60 + th, fh2, sine, cose); hsp->cusp[2] = Asc1(120 + th, fh2, sine, cose); hsp->cusp[3] = Asc1(150 + th, fh1, sine, cose); + if (hsp->do_hspeed) { + hsp->cusp_speed[11] = AscDash(30 + th, fh1, sine, cose); + hsp->cusp_speed[12] = AscDash(60 + th, fh2, sine, cose); + hsp->cusp_speed[2] = AscDash(120 + th, fh2, sine, cose); + hsp->cusp_speed[3] = AscDash(150 + th, fh1, sine, cose); + } /* within polar circle, when mc sinks below horizon and * ascendant changes to western hemisphere, all cusps * must be added 180 degrees. @@ -1125,6 +1296,7 @@ static int CalcH( hsp->cusp[11] = swe_degnorm(hsp->mc + s4 * 0.5); hsp->cusp[12] = swe_degnorm(hsp->mc + s4 * 1.5); } + hsp->do_interpol = hsp->do_hspeed; break; case 'T': /* 'topocentric' houses */ fh1 = atand (tanfi / 3.0); @@ -1133,6 +1305,12 @@ static int CalcH( hsp->cusp[12] = Asc1(60 + th, fh2, sine, cose); hsp->cusp[2] = Asc1(120 + th, fh2, sine, cose); hsp->cusp[3] = Asc1(150 + th, fh1, sine, cose); + if (hsp->do_hspeed) { + hsp->cusp_speed[11] = AscDash(30 + th, fh1, sine, cose); + hsp->cusp_speed[12] = AscDash(60 + th, fh2, sine, cose); + hsp->cusp_speed[2] = AscDash(120 + th, fh2, sine, cose); + hsp->cusp_speed[3] = AscDash(150 + th, fh1, sine, cose); + } /* within polar circle, when mc sinks below horizon and * ascendant changes to western hemisphere, all cusps * must be added 180 degrees. @@ -1158,6 +1336,11 @@ static int CalcH( for (i = 2; i <=12; i++) hsp->cusp[i] = swe_degnorm(hsp->cusp[1] + (i-1) * 30); break; + if (hsp->do_hspeed) { + for (i = 1; i <=12; i++) { + hsp->cusp_speed[i] = hsp->ac_speed; + } + } case 'W': /* equal, whole-sign houses */ acmc = swe_difdeg2n(hsp->ac, hsp->mc); if (acmc < 0) { @@ -1199,6 +1382,7 @@ static int CalcH( if (acmc < 0) { hsp->ac = swe_degnorm(hsp->ac + 180); } + hsp->do_interpol = hsp->do_hspeed; break; } case 'M': { /* @@ -1222,6 +1406,7 @@ static int CalcH( if (acmc < 0) { hsp->ac = swe_degnorm(hsp->ac + 180); } + hsp->do_interpol = hsp->do_hspeed; break; } case 'F': { /* @@ -1261,6 +1446,7 @@ static int CalcH( hsp->cusp[i] = swe_degnorm(hsp->cusp[i]); } } + hsp->do_interpol = hsp->do_hspeed; break; } case 'B': { /* Alcabitius */ /* created by Alois 17-sep-2000, followed example in Matrix @@ -1297,10 +1483,12 @@ static int CalcH( rectasc = swe_degnorm(th + 180 - sn3); /* cusp 3 */ hsp->cusp[3] = Asc1(rectasc, 0, sine, cose); } + hsp->do_interpol = hsp->do_hspeed; break; case 'G': /* 36 Gauquelin sectors */ for (i = 1; i <= 36; i++) { hsp->cusp[i] = 0; + hsp->cusp_speed[i] = 0; } if (fabs(fi) >= 90 - ekl) { /* within polar circle */ retc = ERR; @@ -1318,6 +1506,7 @@ static int CalcH( tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose)))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { /* pole height */ f = atand(sind(asind(tanfi * tant) * ih2 / 9) /tant); @@ -1327,6 +1516,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } /* pole height */ @@ -1345,8 +1535,10 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); } hsp->cusp[ih+18] = swe_degnorm(hsp->cusp[ih] + 180); + if (hsp->do_hspeed) hsp->cusp_speed[ih + 18] = hsp->cusp_speed[ih]; } /*************** first/third quarter ***************/ for (ih = 29; ih <= 36; ih++) { @@ -1356,6 +1548,7 @@ static int CalcH( tant = tand(asind(sine * sind(Asc1(rectasc, fh1, sine, cose)))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { f = atand(sind(asind(tanfi * tant) * ih2 / 9) / tant); /* pole height */ @@ -1365,6 +1558,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } f = atand(sind(asind(tanfi * tant) * ih2 / 9) / tant); @@ -1383,13 +1577,21 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); } hsp->cusp[ih-18] = swe_degnorm(hsp->cusp[ih] + 180); + if (hsp->do_hspeed) hsp->cusp_speed[ih - 18] = hsp->cusp_speed[ih]; } hsp->cusp[1] = hsp->ac; hsp->cusp[10] = hsp->mc; hsp->cusp[19] = swe_degnorm(hsp->ac + 180); hsp->cusp[28] = swe_degnorm(hsp->mc + 180); + if (hsp->do_hspeed) { + hsp->cusp_speed[1] = hsp->ac_speed; + hsp->cusp_speed[10] = hsp->mc_speed; + hsp->cusp_speed[19] = hsp->ac_speed; + hsp->cusp_speed[28] = hsp->mc_speed; + } break; case 'U': /* Krusinski-Pisa */ /* @@ -1488,6 +1690,7 @@ static int CalcH( hsp->cusp[i] = swe_degnorm(hsp->cusp[i] + 180); } } + hsp->do_interpol = hsp->do_hspeed; break; default: /* Placidus houses */ if (fabs(fi) >= 90 - ekl) { /* within polar circle */ @@ -1504,6 +1707,7 @@ static int CalcH( ih = 11; if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { /* pole height */ f = atand(sind(asind(tanfi * tant) / 3) /tant); @@ -1513,6 +1717,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } /* pole height */ @@ -1527,6 +1732,7 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); #ifdef DEBUG_PLAC_ITER fprintf(stderr, "h=%d, niter=%d\n", ih, i); #endif @@ -1537,6 +1743,7 @@ static int CalcH( ih = 12; if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { f = atand(sind(asind(tanfi * tant) / 1.5) / tant); /* pole height */ @@ -1546,6 +1753,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } f = atand(sind(asind(tanfi * tant) / 1.5) / tant); @@ -1560,6 +1768,7 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); #ifdef DEBUG_PLAC_ITER fprintf(stderr, "h=%d, niter=%d\n", ih, i); #endif @@ -1570,6 +1779,7 @@ static int CalcH( ih = 2; if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { f = atand(sind(asind(tanfi * tant) / 1.5) / tant); /* pole height */ @@ -1579,6 +1789,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } f = atand(sind(asind(tanfi * tant) / 1.5) / tant); @@ -1593,6 +1804,7 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); #ifdef DEBUG_PLAC_ITER fprintf(stderr, "h=%d, niter=%d\n", ih, i); #endif @@ -1603,6 +1815,7 @@ static int CalcH( ih = 3; if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; } else { f = atand(sind(asind(tanfi * tant) / 3) / tant); /* pole height */ @@ -1612,6 +1825,7 @@ static int CalcH( tant = tand(asind(sine * sind(hsp->cusp[ih]))); if (fabs(tant) < VERY_SMALL) { hsp->cusp[ih] = rectasc; + if (hsp->do_hspeed) hsp->cusp_speed[ih] = hsp->armc_speed; break; } f = atand(sind(asind(tanfi * tant) / 3) / tant); @@ -1626,6 +1840,7 @@ static int CalcH( strcpy(hsp->serr, "very close to polar circle, switched to Porphyry"); goto porphyry; } + if (hsp->do_hspeed) hsp->cusp_speed[ih] = AscDash(rectasc, f, sine, cose); #ifdef DEBUG_PLAC_ITER fprintf(stderr, "h=%d, niter=%d\n", ih, i); #endif @@ -1639,6 +1854,14 @@ static int CalcH( hsp->cusp[7] = swe_degnorm(hsp->cusp[1] + 180); hsp->cusp[8] = swe_degnorm(hsp->cusp[2] + 180); hsp->cusp[9] = swe_degnorm(hsp->cusp[3] + 180); + if (hsp->do_hspeed && ! hsp->do_interpol) { + hsp->cusp_speed[4] = hsp->cusp_speed[10]; + hsp->cusp_speed[5] = hsp->cusp_speed[11]; + hsp->cusp_speed[6] = hsp->cusp_speed[12]; + hsp->cusp_speed[7] = hsp->cusp_speed[1]; + hsp->cusp_speed[8] = hsp->cusp_speed[2]; + hsp->cusp_speed[9] = hsp->cusp_speed[3]; + } } /* vertex */ if (fi >= 0) @@ -1646,6 +1869,7 @@ static int CalcH( else f = -90 - fi; hsp->vertex = Asc1(th - 90, f, sine, cose); + if (hsp->do_speed) hsp->vertex_speed = AscDash(th - 90, f, sine, cose); /* with tropical latitudes, the vertex behaves strange, * in a similar way as the ascendant within the polar * circle. we keep it always on the western hemisphere.*/ @@ -1672,15 +1896,21 @@ static int CalcH( hsp->equasc = 270; } /* if */ hsp->equasc = swe_degnorm(hsp->equasc); + if (hsp->do_speed) hsp->equasc_speed = AscDash(th + 90, 0, sine, cose); /* "co-ascendant" W. Koch */ hsp->coasc1 = swe_degnorm(Asc1(th - 90, fi, sine, cose) + 180); + if (hsp->do_speed) hsp->coasc1_speed = AscDash(th - 90, fi, sine, cose); /* "co-ascendant" M. Munkasey */ - if (fi >= 0) + if (fi >= 0) { hsp->coasc2 = Asc1(th + 90, 90 - fi, sine, cose); - else /* southern hemisphere */ + if (hsp->do_speed) hsp->coasc2_speed = AscDash(th + 90, 90 - fi, sine, cose); + } else { /* southern hemisphere */ hsp->coasc2 = Asc1(th + 90, -90 - fi, sine, cose); + if (hsp->do_speed) hsp->coasc2_speed = AscDash(th + 90, -90 - fi, sine, cose); + } /* "polar ascendant" M. Munkasey */ hsp->polasc = Asc1(th - 90, fi, sine, cose); + if (hsp->do_speed) hsp->polasc_speed = AscDash(th - 90, fi, sine, cose); #if 0 test_Asc1(); #endif @@ -1720,6 +1950,25 @@ static double Asc1(double x1, double f, double sine, double cose) return ass; } /* Asc1 */ +// derivative of Asc1, computes speed +// code crontributed by Graham Dawson +static double AscDash(double x, double f, double sine, double cose) +{ + double cosx = cosd(x); + double sinx = sind(x); + double sinx2 = sinx * sinx; + double c = cose * cosx - tand(f) * sine; + double d = sinx2 + c * c; + double dudt; + if (d > VERY_SMALL) { + dudt = (cosx * c + cose * sinx2) / d; + } else { + dudt = 0.0; // When we are on axis of ecliptic + } + return dudt * ARMCS; // 360.985647366; +} + + #if 0 /******************************/ static double Asc1_old(double x1, double f, double sine, double cose) @@ -1916,7 +2165,7 @@ if (1) { // which we do not know. If it sees ascmc[9] == 99, it uses // the one is saved from last call. can lead to bugs, but can // also solve many problems. - if (swe_houses_armc(armc, geolat, eps, hsys, hcusp, ascmc) == ERR) { + if (swe_houses_armc_ex2(armc, geolat, eps, hsys, hcusp, ascmc, NULL, NULL, serr) == ERR) { if (serr != NULL) sprintf(serr, "swe_house_pos(): failed for system %c", hsys); } else { @@ -2489,7 +2738,7 @@ if (1) { break; default: hpos = 0; - if (swe_houses_armc(armc, geolat, eps, hsys, hcusp, ascmc) == ERR) { + if (swe_houses_armc_ex2(armc, geolat, eps, hsys, hcusp, ascmc, NULL, NULL, serr) == ERR) { if (serr != NULL) sprintf(serr, "swe_house_pos(): failed for system %c", hsys); break; diff --git a/deps/swisseph/swehouse.h b/deps/swisseph/swehouse.h index a75ccdd..0acf5c2 100644 --- a/deps/swisseph/swehouse.h +++ b/deps/swisseph/swehouse.h @@ -60,14 +60,26 @@ house and (simple) aspect calculation struct houses { double cusp[37]; + double cusp_speed[37]; double ac; + double ac_speed; // speed of ac double mc; + double mc_speed; // speed of mc + double armc_speed; // speed of armc double vertex; + double vertex_speed; // speed of vertex double equasc; + double equasc_speed; // speed double coasc1; + double coasc1_speed; // speed double coasc2; + double coasc2_speed; // speed double polasc; + double polasc_speed; // speed double sundec; // declination of Sun for Sunshine houses + AS_BOOL do_speed; + AS_BOOL do_hspeed; + AS_BOOL do_interpol; char serr[AS_MAXCH]; }; diff --git a/deps/swisseph/sweph.c b/deps/swisseph/sweph.c index 97b9869..e76ae31 100644 --- a/deps/swisseph/sweph.c +++ b/deps/swisseph/sweph.c @@ -321,6 +321,8 @@ int32 CALL_CONV swe_calc(double tjd, int ipl, int32 iflag, struct save_positions *sd; double x[6], *xs, x0[24], x2[24]; double dt; + if (serr != NULL) + *serr = '\0'; #ifdef TRACE #ifdef FORCE_IFLAG /* @@ -343,8 +345,6 @@ int32 CALL_CONV swe_calc(double tjd, int ipl, int32 iflag, FILE *fp; char s[AS_MAXCH], *sp; memset(x, 0, sizeof(double) * 6); - if (serr != NULL) - *serr = '\0'; /* if the following file exists, flag is read from it and or'ed into iflag */ if (!force_flag_checked) { if ((fp = fopen(fname_force_flg, BFILE_R_ACCESS)) != NULL) { @@ -1021,7 +1021,7 @@ static int32 swecalc(double tjd, int ipl, int32 iflag, double *x, char *serr) /* internal planet number */ if (ipl < SE_NPLANETS) ipli = pnoext2int[ipl]; - else if (ipl <= SE_AST_OFFSET + MPC_VESTA) { + else if (ipl <= SE_AST_OFFSET + MPC_VESTA && ipl > SE_AST_OFFSET) { ipli = SEI_CERES + ipl - SE_AST_OFFSET - 1; ipl = SE_CERES + ipl - SE_AST_OFFSET - 1; #if 0 @@ -3134,7 +3134,7 @@ int32 swi_get_ayanamsa_ex(double tjd_et, int32 iflag, double *daya, char *serr) t0 += swe_deltat_ex(t0, iflag, serr); swi_precess(x, t0, 0, J2000_TO_J); /* to ecliptic t0 */ - eps = swi_epsiln(t0, 0); + eps = swi_epsiln(t0, iflag); swi_coortrf(x, x, eps); /* to polar */ swi_cartpol(x, x); @@ -3153,7 +3153,7 @@ int32 swi_get_ayanamsa_ex(double tjd_et, int32 iflag, double *daya, char *serr) t0 = sip->t0; if (sip->t0_is_UT) t0 += swe_deltat_ex(t0, iflag, serr); - eps = swi_epsiln(t0, 0); + eps = swi_epsiln(t0, iflag); // to polar equatorial relative to equinox t0 swi_polcart(x, x); swi_coortrf(x, x, -eps); @@ -3163,7 +3163,7 @@ int32 swi_get_ayanamsa_ex(double tjd_et, int32 iflag, double *daya, char *serr) // precess to date swi_precess(x, tjd_et, 0, J2000_TO_J); // epsilon of date - eps = swi_epsiln(tjd_et, 0); + eps = swi_epsiln(tjd_et, iflag); // to polar swi_coortrf(x, x, eps); swi_cartpol(x, x); @@ -6735,27 +6735,27 @@ static AS_BOOL get_builtin_star(char *star, char *sstar, char *srecord) /* some stars are built-in, because they are required for Hindu * sidereal ephemerides */ /* Ayanamsha SE_SIDM_TRUE_CITRA */ - if (strncmp(star, "spica", 5) == 0) { + if (strncmp(star, "spica", 5) == 0 || strncmp(star, "Spica", 5) == 0) { strcpy(srecord, "Spica,alVir,ICRS,13,25,11.57937,-11,09,40.7501,-42.35,-30.67,1,13.06,0.97,-10,3672"); strcpy(sstar, "spica"); return TRUE; /* Ayanamsha SE_SIDM_TRUE_REVATI */ - } else if (strstr(star, ",zePsc") != NULL || strncmp(star, "revati", 6) == 0) { + } else if (strstr(star, ",zePsc") != NULL || strncmp(star, "revati", 6) == 0 || strncmp(star, "Revati", 6) == 0) { strcpy(srecord, "Revati,zePsc,ICRS,01,13,43.88735,+07,34,31.2745,145,-55.69,15,18.76,5.187,06,174"); strcpy(sstar, "revati"); return TRUE; /* Ayanamsha SE_SIDM_TRUE_PUSHYA */ - } else if (strstr(star, ",deCnc") != NULL || strncmp(star, "pushya", 6) == 0) { + } else if (strstr(star, ",deCnc") != NULL || strncmp(star, "pushya", 6) == 0 || strncmp(star, "Pushya", 6) == 0 ) { strcpy(srecord, "Pushya,deCnc,ICRS,08,44,41.09921,+18,09,15.5034,-17.67,-229.26,17.14,24.98,3.94,18,2027"); strcpy(sstar, "pushya"); return TRUE; /* Ayanamsha SE_SIDM_TRUE_SHEORAN */ - } else if (strstr(star, ",deCnc") != NULL || strncmp(star, "pushya", 6) == 0) { + } else if (strstr(star, ",deCnc") != NULL) { strcpy(srecord, "Pushya,deCnc,ICRS,08,44,41.09921,+18,09,15.5034,-17.67,-229.26,17.14,24.98,3.94,18,2027"); strcpy(sstar, "pushya"); return TRUE; /* Ayanamsha SE_SIDM_TRUE_MULA */ - } else if (strstr(star, ",laSco") != NULL || strncmp(star, "mula", 6) == 0) { + } else if (strstr(star, ",laSco") != NULL || strncmp(star, "mula", 6) == 0 || strncmp(star, "Mula", 6) == 0) { strcpy(srecord, "Mula,laSco,ICRS,17,33,36.52012,-37,06,13.7648,-8.53,-30.8,-3,5.71,1.62,-37,11673"); strcpy(sstar, "mula"); return TRUE; @@ -7059,7 +7059,7 @@ char *CALL_CONV swe_get_planet_name(int ipl, char *s) * 2. asteroid name * The asteroid number may or may not be in brackets */ - if (s[0] == '?' || isdigit((int) s[1])) { + if (ipl > SE_AST_OFFSET && (s[0] == '?' || isdigit((int) s[1]))) { int ipli = (int) (ipl - SE_AST_OFFSET), iplf = 0; FILE *fp; char si[AS_MAXCH], *sp, *sp2; @@ -7257,7 +7257,7 @@ void swi_force_app_pos_etc() swed.pldat[i].xflgs = -1; for (i = 0; i < SEI_NNODE_ETC; i++) swed.nddat[i].xflgs = -1; - for (i = 0; i < SE_NPLANETS; i++) { + for (i = 0; i <= SE_NPLANETS; i++) { // "=" because save area for asteroids > SE_AST_OFFSET is at i == SE_NPLANETS swed.savedat[i].tsave = 0; swed.savedat[i].iflgsave = -1; } diff --git a/deps/swisseph/sweph.h b/deps/swisseph/sweph.h index 830d069..b39ef34 100644 --- a/deps/swisseph/sweph.h +++ b/deps/swisseph/sweph.h @@ -62,7 +62,7 @@ * move over from swephexp.h */ -#define SE_VERSION "2.09.01" // "2.08.00b" +#define SE_VERSION "2.09.03" #define J2000 2451545.0 /* 2000 January 1.5 */ #define B1950 2433282.42345905 /* 1950 January 0.923 */ @@ -175,6 +175,7 @@ #define SEI_FILE_MAIN_AST 2 #define SEI_FILE_ANY_AST 3 #define SEI_FILE_FIXSTAR 4 +#define SEI_FILE_PLMOON 5 #if 0 #define SEI_FILE_TEST_ENDIAN (97L * 65536L + 98L * 256L + 99L) /*abc*/ @@ -466,7 +467,7 @@ static const struct aya_init ayanamsa[] = { {1684532.5, -4.44138598, TRUE, 0}, /*************************/ /* 15: Hipparchos */ -{1674484, -9.33333, TRUE, -1}, // 15: Hipparchos +{1674484.0, -9.33333, TRUE, -1}, // 15: Hipparchos /*************************/ /* 16: Sassanian */ {1927135.8747793, 0, TRUE, -1}, // 16: Sassanian diff --git a/deps/swisseph/swephexp.h b/deps/swisseph/swephexp.h index 75caf86..7352c01 100644 --- a/deps/swisseph/swephexp.h +++ b/deps/swisseph/swephexp.h @@ -124,6 +124,7 @@ extern "C" { #define SE_NPLANETS 23 +#define SE_PLMOON_OFFSET 9000 #define SE_AST_OFFSET 10000 #define SE_VARUNA (SE_AST_OFFSET + 20000) @@ -804,10 +805,18 @@ ext_def( int ) swe_houses_ex( double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusps, double *ascmc); +ext_def( int ) swe_houses_ex2( + double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, + double *cusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr); + ext_def( int ) swe_houses_armc( double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc); +ext_def( int ) swe_houses_armc_ex2( + double armc, double geolat, double eps, int hsys, + double *cusps, double *ascmc, double *cusp_speed, double *ascmc_speed, char *serr); + ext_def(double) swe_house_pos( double armc, double geolat, double eps, int hsys, double *xpin, char *serr); diff --git a/deps/swisseph/swephlib.c b/deps/swisseph/swephlib.c index 7013c5d..1c81225 100644 --- a/deps/swisseph/swephlib.c +++ b/deps/swisseph/swephlib.c @@ -2550,7 +2550,7 @@ static int32 calc_deltat(double tjd, int32 iflag, double *deltat, char *serr) int32 retc; int deltat_model = swed.astro_models[SE_MODEL_DELTAT]; double tid_acc; - int32 denumret; + int32 denum, denumret; int32 epheflag, otherflag; //fprintf(stderr, "dmod=%f, %.f\n", (double) deltat_model, (double) SEMOD_DELTAT_DEFAULT); if (deltat_model == 0) deltat_model = SEMOD_DELTAT_DEFAULT; @@ -2561,12 +2561,14 @@ static int32 calc_deltat(double tjd, int32 iflag, double *deltat, char *serr) retc = swi_get_tid_acc(tjd, 0, 9999, &denumret, &tid_acc, serr); /* for default tid_acc */ /* otherwise we use tid_acc consistent with epheflag */ } else { + denum = swed.jpldenum; + if (epheflag & SEFLG_SWIEPH) denum = swed.fidat[SEI_FILE_MOON].sweph_denum; if (swi_init_swed_if_start() == 1 && !(epheflag & SEFLG_MOSEPH)) { if (serr != NULL) strcpy(serr, "Please call swe_set_ephe_path() or swe_set_jplfile() before calling swe_deltat_ex()"); - retc = swi_set_tid_acc(tjd, epheflag, 0, NULL); /* _set_ saves tid_acc in swed */ + retc = swi_set_tid_acc(tjd, epheflag, denum, NULL); /* _set_ saves tid_acc in swed */ } else { - retc = swi_set_tid_acc(tjd, epheflag, 0, serr); /* _set_ saves tid_acc in swed */ + retc = swi_set_tid_acc(tjd, epheflag, denum, serr); /* _set_ saves tid_acc in swed */ } tid_acc = swed.tid_acc; } @@ -3195,12 +3197,13 @@ int32 swi_guess_ephe_flag() int32 swi_get_tid_acc(double tjd_ut, int32 iflag, int32 denum, int32 *denumret, double *tid_acc, char *serr) { - double xx[6], tjd_et; + //double xx[6], tjd_et; iflag &= SEFLG_EPHMASK; if (swed.is_tid_acc_manual) { *tid_acc = swed.tid_acc; return iflag; } +//denum = 0; if (denum == 0) { if (iflag & SEFLG_MOSEPH) { *tid_acc = SE_TIDAL_DE404; @@ -3210,33 +3213,12 @@ int32 swi_get_tid_acc(double tjd_ut, int32 iflag, int32 denum, int32 *denumret, if (iflag & SEFLG_JPLEPH) { if (swed.jpl_file_is_open) { denum = swed.jpldenum; - } else { - tjd_et = tjd_ut; /* + swe_deltat_ex(tjd_ut, 0, NULL); we do not add - delta t, because it would result in a recursive - call of swi_set_tid_acc() */ - iflag = SEFLG_JPLEPH|SEFLG_J2000|SEFLG_TRUEPOS|SEFLG_ICRS|SEFLG_BARYCTR; - iflag = swe_calc(tjd_et, SE_JUPITER, iflag, xx, serr); - if (swed.jpl_file_is_open && (iflag & SEFLG_JPLEPH)) { - denum = swed.jpldenum; - } } } /* SEFLG_SWIEPH wanted or SEFLG_JPLEPH failed: */ - if (denum == 0) { - tjd_et = tjd_ut; /* + swe_deltat_ex(tjd_ut, 0, NULL); we do not add - delta t, because it would result in a recursive - call of swi_set_tid_acc() */ - if (swed.fidat[SEI_FILE_MOON].fptr == NULL || - tjd_et < swed.fidat[SEI_FILE_MOON].tfstart + 1 || - tjd_et > swed.fidat[SEI_FILE_MOON].tfend - 1) { - iflag = SEFLG_SWIEPH|SEFLG_J2000|SEFLG_TRUEPOS|SEFLG_ICRS; - iflag = swe_calc(tjd_et, SE_MOON, iflag, xx, serr); - } + if (iflag & SEFLG_SWIEPH) { if (swed.fidat[SEI_FILE_MOON].fptr != NULL) { denum = swed.fidat[SEI_FILE_MOON].sweph_denum; - /* Moon ephemeris file is not available, default to Moshier ephemeris */ - } else { - denum = 404; /* DE number of Moshier ephemeris */ } } } @@ -3261,6 +3243,7 @@ int32 swi_set_tid_acc(double tjd_ut, int32 iflag, int32 denum, char *serr) { int32 retc = iflag; int32 denumret; + //fprintf(stderr, "ifl=%d\n", iflag); /* manual tid_acc overrides automatic tid_acc */ if (swed.is_tid_acc_manual) return retc; diff --git a/src/house.cc b/src/house.cc index 5fe2fe8..fc87bb7 100644 --- a/src/house.cc +++ b/src/house.cc @@ -162,6 +162,113 @@ NAN_METHOD(node_swe_houses_ex) { info.GetReturnValue().Set (result); }; +/** + * int swe_houses_ex2(double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) + * => + * int swe_houses_ex2(double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) { + * house: [double], + * ascendant: double, + * mc: double, + * armc: double, + * vertex: double, + * equatorialAscendant: double, + * kochCoAscendant: double, + * munkaseyCoAscendant: double, + * munkaseyPolarAscendant: double, + * houseSpeed: [double], + * ascendantSpeed: double, + * mcSpeed: double, + * armcSpeed: double, + * vertexSpeed: double, + * equatorialAscendantSpeed: double, + * kochCoAscendantSpeed: double, + * munkaseyCoAscendantSpeed: double, + * munkaseyPolarAscendantSpeed: double, + * error: string + * } + */ +NAN_METHOD(node_swe_houses_ex2) { + Nan::HandleScope scope; + + if (info.Length () < 5) { + Nan::ThrowTypeError ("Wrong number of arguments"); + }; + + if ( + !info [0]->IsNumber () || + !info [1]->IsNumber () || + !info [2]->IsNumber () || + !info [3]->IsNumber () || + !info [4]->IsString () + ) { + Nan::ThrowTypeError ("Wrong type of arguments"); + }; + + double cusps [40] = {0}; + double ascmc [40] = {0}; + double cusps_speed [40] = {0}; + double ascmc_speed [40] = {0}; + int rflag; + size_t cuspsCount; + int hsys; + char serr [AS_MAXCH]; + + Local result = Nan::New (); + Local house = Nan::New (); + Local houseSpeed = Nan::New (); + + hsys = (* String::Utf8Value (Isolate::GetCurrent(), info [4]->ToString (Nan::GetCurrentContext()).FromMaybe(v8::Local()))) [0]; + + if (hsys == 'G') { + cuspsCount = 36; + } else { + cuspsCount = 12; + } + + rflag = ::swe_houses_ex2 ( + info [0]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + (int)info [1]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + info [2]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + info [3]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + hsys, + cusps, ascmc, cusps_speed, ascmc_speed, + serr + ); + + if (rflag < 0) { + Nan::Set(result,Nan::New ("error").ToLocalChecked(), Nan::New ("Can't calculate houses.").ToLocalChecked()); + Nan::Set(result,Nan::New ("message").ToLocalChecked(), Nan::New (serr).ToLocalChecked()); + } else { + for (size_t i = 0; i < cuspsCount; i ++) { + Nan::Set(house,Nan::New (i), Nan::New (cusps [i + 1])); + Nan::Set(houseSpeed,Nan::New (i), Nan::New (cusps_speed [i + 1])); + }; + + Nan::Set(result,Nan::New ("house").ToLocalChecked(), house); + Nan::Set(result,Nan::New ("ascendant").ToLocalChecked(), Nan::New (ascmc [SE_ASC])); + Nan::Set(result,Nan::New ("mc").ToLocalChecked(), Nan::New (ascmc [SE_MC])); + Nan::Set(result,Nan::New ("armc").ToLocalChecked(), Nan::New (ascmc [SE_ARMC])); + Nan::Set(result,Nan::New ("vertex").ToLocalChecked(), Nan::New (ascmc [SE_VERTEX])); + Nan::Set(result,Nan::New ("equatorialAscendant").ToLocalChecked(), Nan::New (ascmc [SE_EQUASC])); + Nan::Set(result,Nan::New ("kochCoAscendant").ToLocalChecked(), Nan::New (ascmc [SE_COASC1])); + Nan::Set(result,Nan::New ("munkaseyCoAscendant").ToLocalChecked(), Nan::New (ascmc [SE_COASC2])); + Nan::Set(result,Nan::New ("munkaseyPolarAscendant").ToLocalChecked(), Nan::New (ascmc [SE_POLASC])); + + Nan::Set(result,Nan::New ("houseSpeed").ToLocalChecked(), houseSpeed); + Nan::Set(result,Nan::New ("ascendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_ASC])); + Nan::Set(result,Nan::New ("mcSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_MC])); + Nan::Set(result,Nan::New ("armcSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_ARMC])); + Nan::Set(result,Nan::New ("vertexSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_VERTEX])); + Nan::Set(result,Nan::New ("equatorialAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_EQUASC])); + Nan::Set(result,Nan::New ("kochCoAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_COASC1])); + Nan::Set(result,Nan::New ("munkaseyCoAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_COASC2])); + Nan::Set(result,Nan::New ("munkaseyPolarAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_POLASC])); + }; + + HandleCallback (info, result); + info.GetReturnValue().Set (result); +}; + /** * int swe_houses_armc(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc) * => @@ -241,6 +348,111 @@ NAN_METHOD(node_swe_houses_armc) { info.GetReturnValue().Set (result); }; +/** + * int swe_houses_armc_ex2(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) + * => + * int swe_houses_armc_ex2(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) { + * house: [double], + * ascendant: double, + * mc: double, + * armc: double, + * vertex: double, + * equatorialAscendant: double, + * kochCoAscendant: double, + * munkaseyCoAscendant: double, + * munkaseyPolarAscendant: double, + * houseSpeed: [double], + * ascendantSpeed: double, + * mcSpeed: double, + * armcSpeed: double, + * vertexSpeed: double, + * equatorialAscendantSpeed: double, + * kochCoAscendantSpeed: double, + * munkaseyCoAscendantSpeed: double, + * munkaseyPolarAscendantSpeed: double, + * error: string + * } + */ +NAN_METHOD(node_swe_houses_armc_ex2) { + Nan::HandleScope scope; + + if (info.Length () < 4) { + Nan::ThrowTypeError ("Wrong number of arguments"); + }; + + if ( + !info [0]->IsNumber () || + !info [1]->IsNumber () || + !info [2]->IsNumber () || + !info [3]->IsString () + ) { + Nan::ThrowTypeError ("Wrong type of arguments"); + }; + + double cusps [40] = {0}; + double ascmc [40] = {0}; + double cusps_speed [40] = {0}; + double ascmc_speed [40] = {0}; + int rflag; + size_t cuspsCount; + int hsys; + char serr [AS_MAXCH]; + + Local result = Nan::New (); + Local house = Nan::New (); + Local houseSpeed = Nan::New (); + + hsys = (* String::Utf8Value (Isolate::GetCurrent(), info [3]->ToString (Nan::GetCurrentContext()).FromMaybe(v8::Local()))) [0]; + + if (hsys == 'G') { + cuspsCount = 36; + } else { + cuspsCount = 12; + } + + rflag = ::swe_houses_armc_ex2 ( + info [0]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + info [1]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + info [2]->NumberValue (Nan::GetCurrentContext()).ToChecked(), + hsys, + cusps, ascmc, cusps_speed, ascmc_speed, + serr + ); + + if (rflag < 0) { + Nan::Set(result,Nan::New ("error").ToLocalChecked(), Nan::New ("Can't calculate houses.").ToLocalChecked()); + Nan::Set(result,Nan::New ("message").ToLocalChecked(), Nan::New (serr).ToLocalChecked()); + } else { + for (size_t i = 0; i < cuspsCount; i ++) { + Nan::Set(house,Nan::New (i), Nan::New (cusps [i + 1])); + Nan::Set(houseSpeed,Nan::New (i), Nan::New (cusps_speed [i + 1])); + }; + + Nan::Set(result,Nan::New ("house").ToLocalChecked(), house); + Nan::Set(result,Nan::New ("ascendant").ToLocalChecked(), Nan::New (ascmc [SE_ASC])); + Nan::Set(result,Nan::New ("mc").ToLocalChecked(), Nan::New (ascmc [SE_MC])); + Nan::Set(result,Nan::New ("armc").ToLocalChecked(), Nan::New (ascmc [SE_ARMC])); + Nan::Set(result,Nan::New ("vertex").ToLocalChecked(), Nan::New (ascmc [SE_VERTEX])); + Nan::Set(result,Nan::New ("equatorialAscendant").ToLocalChecked(), Nan::New (ascmc [SE_EQUASC])); + Nan::Set(result,Nan::New ("kochCoAscendant").ToLocalChecked(), Nan::New (ascmc [SE_COASC1])); + Nan::Set(result,Nan::New ("munkaseyCoAscendant").ToLocalChecked(), Nan::New (ascmc [SE_COASC2])); + Nan::Set(result,Nan::New ("munkaseyPolarAscendant").ToLocalChecked(), Nan::New (ascmc [SE_POLASC])); + + Nan::Set(result,Nan::New ("houseSpeed").ToLocalChecked(), houseSpeed); + Nan::Set(result,Nan::New ("ascendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_ASC])); + Nan::Set(result,Nan::New ("mcSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_MC])); + Nan::Set(result,Nan::New ("armcSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_ARMC])); + Nan::Set(result,Nan::New ("vertexSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_VERTEX])); + Nan::Set(result,Nan::New ("equatorialAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_EQUASC])); + Nan::Set(result,Nan::New ("kochCoAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_COASC1])); + Nan::Set(result,Nan::New ("munkaseyCoAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_COASC2])); + Nan::Set(result,Nan::New ("munkaseyPolarAscendantSpeed").ToLocalChecked(), Nan::New (ascmc_speed [SE_POLASC])); + }; + + HandleCallback (info, result); + info.GetReturnValue().Set (result); +}; + /** * double swe_house_pos(double armc, double geolat, double eps, int hsys, double *xpin, char *serr) * => diff --git a/src/house.h b/src/house.h index 638259b..787163d 100644 --- a/src/house.h +++ b/src/house.h @@ -42,6 +42,33 @@ NAN_METHOD(node_swe_houses); */ NAN_METHOD(node_swe_houses_ex); +/** + * int swe_houses_ex2(double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) + * => + * int swe_houses_ex2(double tjd_ut, int32 iflag, double geolat, double geolon, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) { + * house: [double], + * ascendant: double, + * mc: double, + * armc: double, + * vertex: double, + * equatorialAscendant: double, + * kochCoAscendant: double, + * munkaseyCoAscendant: double, + * munkaseyPolarAscendant: double, + * houseSpeed: [double], + * ascendantSpeed: double, + * mcSpeed: double, + * armcSpeed: double, + * vertexSpeed: double, + * equatorialAscendantSpeed: double, + * kochCoAscendantSpeed: double, + * munkaseyCoAscendantSpeed: double, + * munkaseyPolarAscendantSpeed: double, + * error: string + * } + */ +NAN_METHOD(node_swe_houses_ex2); + /** * int swe_houses_armc(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc) * => @@ -60,6 +87,33 @@ NAN_METHOD(node_swe_houses_ex); */ NAN_METHOD(node_swe_houses_armc); +/** + * int swe_houses_armc2(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) + * => + * int swe_houses_armc2(double armc, double geolat, double eps, int hsys, double *cusps, double *ascmc, double *cusps_speed, double *ascmc_speed, char *serr) { + * house: [double], + * ascendant: double, + * mc: double, + * armc: double, + * vertex: double, + * equatorialAscendant: double, + * kochCoAscendant: double, + * munkaseyCoAscendant: double, + * munkaseyPolarAscendant: double, + * houseSpeed: [double], + * ascendantSpeed: double, + * mcSpeed: double, + * armcSpeed: double, + * vertexSpeed: double, + * equatorialAscendantSpeed: double, + * kochCoAscendantSpeed: double, + * munkaseyCoAscendantSpeed: double, + * munkaseyPolarAscendantSpeed: double, + * error: string + * } + */ +NAN_METHOD(node_swe_houses_armc_ex2); + /** * double swe_house_pos(double armc, double geolat, double eps, int hsys, double *xpin, char *serr) * => diff --git a/src/swisseph.cc b/src/swisseph.cc index 8426eee..96be6ed 100644 --- a/src/swisseph.cc +++ b/src/swisseph.cc @@ -73,7 +73,9 @@ void Initialize (Local exports) { // house Nan::SetMethod (exports, "swe_houses", node_swe_houses); Nan::SetMethod (exports, "swe_houses_ex", node_swe_houses_ex); + Nan::SetMethod (exports, "swe_houses_ex2", node_swe_houses_ex2); Nan::SetMethod (exports, "swe_houses_armc", node_swe_houses_armc); + Nan::SetMethod (exports, "swe_houses_armc_ex2", node_swe_houses_armc_ex2); Nan::SetMethod (exports, "swe_houses_pos", node_swe_houses_pos); // eclipse