584 lines
18 KiB
C
584 lines
18 KiB
C
/*****************************************************
|
|
$Header: /home/dieter/sweph/RCS/swepcalc.c,v 1.74 2008/06/16 10:07:20 dieter Exp $
|
|
Placalc compatibility interface for Swiss Ephemeris.
|
|
|
|
*******************************************************/
|
|
/* Copyright (C) 1997 - 2008 Astrodienst AG, Switzerland. All rights reserved.
|
|
|
|
License conditions
|
|
------------------
|
|
|
|
This file is part of Swiss Ephemeris.
|
|
|
|
Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
|
|
or distributor accepts any responsibility for the consequences of using it,
|
|
or for whether it serves any particular purpose or works at all, unless he
|
|
or she says so in writing.
|
|
|
|
Swiss Ephemeris is made available by its authors under a dual licensing
|
|
system. The software developer, who uses any part of Swiss Ephemeris
|
|
in his or her software, must choose between one of the two license models,
|
|
which are
|
|
a) GNU public license version 2 or later
|
|
b) Swiss Ephemeris Professional License
|
|
|
|
The choice must be made before the software developer distributes software
|
|
containing parts of Swiss Ephemeris to others, and before any public
|
|
service using the developed software is activated.
|
|
|
|
If the developer choses the GNU GPL software license, he or she must fulfill
|
|
the conditions of that license, which includes the obligation to place his
|
|
or her whole software project under the GNU GPL or a compatible license.
|
|
See http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
|
|
|
|
If the developer choses the Swiss Ephemeris Professional license,
|
|
he must follow the instructions as found in http://www.astro.com/swisseph/
|
|
and purchase the Swiss Ephemeris Professional Edition from Astrodienst
|
|
and sign the corresponding license contract.
|
|
|
|
The License grants you the right to use, copy, modify and redistribute
|
|
Swiss Ephemeris, but only under certain conditions described in the License.
|
|
Among other things, the License requires that the copyright notices and
|
|
this notice be preserved on all copies.
|
|
|
|
Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
|
|
|
|
The authors of Swiss Ephemeris have no control or influence over any of
|
|
the derived works, i.e. over software or services created by other
|
|
programmers which use Swiss Ephemeris functions.
|
|
|
|
The names of the authors or of the copyright holder (Astrodienst) must not
|
|
be used for promoting any software, product or service which uses or contains
|
|
the Swiss Ephemeris. This copyright notice is the ONLY place where the
|
|
names of the authors can legally appear, except in cases where they have
|
|
given special permission in writing.
|
|
|
|
The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
|
|
for promoting such software, products or services.
|
|
*/
|
|
|
|
/*
|
|
* This file is the PLACALC compatibility interface for Swiss Ephemeris.
|
|
* It allows very easy porting of older Placalc application to the SwissEph.
|
|
* A user has to replace #include "placalc.h" and "housasp.h" with
|
|
* #include "swepcalc.h"
|
|
* If he has used "ourdef.h" he replaces it with "sweodef.h".
|
|
* Then he links his application with swepcalc.o and runs it against the
|
|
* Swiss Ephemeris DLL or linkable library.
|
|
*
|
|
* All calls which were present in the placalc sources are contained
|
|
* here, and either implemented directly or translated into Swiss Ephemeris
|
|
* calls.
|
|
*
|
|
*
|
|
*/
|
|
|
|
#include "swepcalc.h"
|
|
#include "swephexp.h"
|
|
|
|
/************************************************************
|
|
local globals, not exported
|
|
************************************************************/
|
|
|
|
static const int plac2swe[] = {SE_SUN, SE_MOON, SE_MERCURY, SE_VENUS, SE_MARS, SE_JUPITER, SE_SATURN, SE_URANUS, SE_NEPTUNE, SE_PLUTO, SE_MEAN_NODE, SE_TRUE_NODE, SE_CHIRON, SE_MEAN_APOG, SE_CERES, SE_PALLAS, SE_JUNO, SE_VESTA,};
|
|
|
|
/* If there occurs an internal error in placalc, a message is
|
|
* written into the string variable perrtx.
|
|
* The message can be read with placalc_get_errtext();
|
|
*/
|
|
static TLS char perrtx[AS_MAXCH];
|
|
static double ekl, nut;
|
|
|
|
/**********************************************************
|
|
function nacalc ()
|
|
calculates an array of planet longitudes and speeds,
|
|
as needed for complete nathan data records.
|
|
The function knows itself how many planets and in which mode
|
|
they have to be calculated for Nathan.
|
|
|
|
return OK or ERR
|
|
|
|
The returned positions are in centiseconds, our standard
|
|
coordinate format for fast mathematics with planetary positions.
|
|
|
|
This function is just a template of how the calc() package
|
|
can be used.
|
|
**********************************************************/
|
|
int nacalc (double jd_ad, /* universal time relative julian date */
|
|
centisec *plon, /* returned longitudes */
|
|
centisec *pspe /* returned speeds, if not NULL pointer */
|
|
)
|
|
{
|
|
char err[AS_MAXCH];
|
|
int planet, flag;
|
|
double rlng, rrad, rlat, rspeed;
|
|
int result = OK;
|
|
flag = CALC_BIT_SPEED; /* same, with speed */
|
|
jd_ad += deltat( jd_ad ); /* ET = UT + Delta_T */
|
|
for (planet = SUN; planet <= MAXPL_NACALC; planet++) {
|
|
if (calc (planet, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
|
|
plon [planet] = swe_csnorm(swe_d2l (rlng * DEG));
|
|
if (pspe != NULL) pspe [planet] = swe_d2l (rspeed * DEG);
|
|
} else {
|
|
plon [planet] = -1;
|
|
if (pspe != NULL) pspe [planet] = 0;
|
|
if (result != ERR) { /* save first error message */
|
|
strcpy(err, placalc_get_errtext());
|
|
}
|
|
result = ERR;
|
|
}
|
|
}
|
|
if (result == ERR)
|
|
strcpy(perrtx, err);
|
|
return result;
|
|
} /* end nacalc */
|
|
|
|
/******************************************************************
|
|
function calc():
|
|
This is the main routine for computing a planets position.
|
|
The function has several modes, which are controlled by bits in
|
|
the parameter 'flag'. The normal mode (flag == 0) computes
|
|
a planets apparent geocentric position in ecliptic coordinates relative to
|
|
the true equinox of date, without speed
|
|
|
|
Explanation of the arguments: see the functions header.
|
|
|
|
Returns OK or ERR (if some planet out of time range). OK and ERR are
|
|
defined in ourdef.h and must not be confused with TRUE and FALSE.
|
|
OK and ERR are of type int, not of type AS_BOOL.
|
|
|
|
Bits used in flag:
|
|
CALC_BIT_HELIO 0 = geocentric, 1 = heliocentric
|
|
CALC_BIT_NOAPP 0 = apparent positions, 1 = true positions
|
|
CALC_BIT_NONUT 0 = do nutation (true equinox of date)
|
|
1 = don't do nutation (mean equinox of date).
|
|
|
|
CALC_BIT_SPEED 0 = don't calc speed,
|
|
1 = calc speed
|
|
|
|
Time range:
|
|
The function can be used savely in the time range 3000 BC to
|
|
3000 AD.
|
|
|
|
Getting ecliptic and nutation:
|
|
calc(CALC_ONLY_ECL_NUT,teph,0,&nutv,&meaneklv,&eklv,NULL);
|
|
will return the values for time teph.
|
|
|
|
******************************************************************/
|
|
int calc(int planet, /* planet index as defined in placalc.h,
|
|
SUN = 0, MOON = 1 etc.
|
|
planet == -1 calc calculates only nut, ekl, meanekl */
|
|
double jd_ad, /* relative Astrodienst Juldate, ephemeris time.
|
|
Astrodienst Juldate is relative 31 Dec 1949, noon. */
|
|
int flag, /* See definition of flag bits above */
|
|
double *alng,
|
|
double *arad,
|
|
double *alat,
|
|
double *alngspeed)
|
|
/* pointers to the return variables:
|
|
alng = ecliptic longitude in degrees
|
|
arad = radius vector in AU (astronomic units)
|
|
alat = ecliptic latitude in degrees
|
|
alngspeed = speed of planet in degrees per day
|
|
*/
|
|
{
|
|
double tjd = jd_ad + JUL_OFFSET;
|
|
double x[6];
|
|
int32 iflagret = 0, iflag = 0;
|
|
int ipl;
|
|
/* planet number
|
|
*/
|
|
/* ecliptic and nutation */
|
|
if (planet == CALC_ONLY_ECL_NUT)
|
|
ipl = SE_ECL_NUT;
|
|
/* earth: placalc makes no difference between sun and earth,
|
|
* swisseph does */
|
|
else if (planet == SUN && (flag & CALC_BIT_HELIO))
|
|
ipl = SE_EARTH;
|
|
else if (planet >= SUN && planet <= VESTA)
|
|
ipl = plac2swe[planet];
|
|
else {
|
|
sprintf(perrtx, "invalid planet number %d. ", planet);
|
|
return ERR;
|
|
}
|
|
/* flag */
|
|
if (flag & CALC_BIT_HELIO)
|
|
if (ipl != SE_MEAN_NODE && ipl != SE_TRUE_NODE && ipl != SE_MEAN_APOG)
|
|
iflag |= SEFLG_HELCTR; /* lunar node and apogee is always geocentric */
|
|
if (flag & CALC_BIT_NOAPP)
|
|
iflag |= SEFLG_TRUEPOS;
|
|
if (flag & CALC_BIT_NONUT)
|
|
iflag |= SEFLG_NONUT;
|
|
if (flag & CALC_BIT_SPEED)
|
|
iflag |= SEFLG_SPEED;
|
|
/* ecliptic and nutation */
|
|
if ((iflagret = swe_calc(tjd, ipl, iflag, x, perrtx)) == ERR)
|
|
return iflagret;
|
|
if (ipl == SE_ECL_NUT) {
|
|
*alng = nut = x[2];
|
|
*arad = x[1];
|
|
*alat = ekl = x[0];
|
|
} else {
|
|
*alng = x[0];
|
|
*arad = x[2];
|
|
*alat = x[1];
|
|
*alngspeed = x[3];
|
|
}
|
|
return (OK);
|
|
} /* end calc */
|
|
|
|
/*
|
|
* get the planet index for an AFL letter
|
|
* returns -1 if the letter does not correspond to a planet.
|
|
*/
|
|
int afl2planet(int afl)
|
|
{
|
|
int p;
|
|
switch (afl) {
|
|
case AFL_SUN : p = SUN; break;
|
|
case AFL_MON : p = MOON; break;
|
|
case AFL_MER : p = MERCURY; break;
|
|
case AFL_VEN : p = VENUS; break;
|
|
case AFL_MAR : p = MARS; break;
|
|
case AFL_JUP : p = JUPITER; break;
|
|
case AFL_SAT : p = SATURN; break;
|
|
case AFL_URA : p = URANUS; break;
|
|
case AFL_NEP : p = NEPTUNE; break;
|
|
case AFL_PLU : p = PLUTO; break;
|
|
case AFL_MNODE : p = MEAN_NODE; break;
|
|
case AFL_TNODE : p = TRUE_NODE; break;
|
|
case AFL_CHI : p = CHIRON; break;
|
|
case AFL_LIL : p = LILITH; break;
|
|
case AFL_CER : p = CERES; break;
|
|
case AFL_PAL : p = PALLAS; break;
|
|
case AFL_JUN : p = JUNO; break;
|
|
case AFL_VES : p = VESTA; break;
|
|
case AFL_AC : p = AC; break;
|
|
case AFL_MC : p = MC; break;
|
|
default : p = -1; break;
|
|
}
|
|
return p;
|
|
}
|
|
|
|
/*
|
|
* get the AFL letter for a planet
|
|
* returns -1 if no letter corresponds to a planet.
|
|
*/
|
|
int planet2afl(int planet)
|
|
{
|
|
switch (planet) {
|
|
case SUN: return AFL_SUN;
|
|
case MOON: return AFL_MON;
|
|
case MERCURY: return AFL_MER;
|
|
case VENUS: return AFL_VEN;
|
|
case MARS: return AFL_MAR;
|
|
case JUPITER: return AFL_JUP;
|
|
case SATURN: return AFL_SAT;
|
|
case URANUS: return AFL_URA;
|
|
case NEPTUNE: return AFL_NEP;
|
|
case PLUTO: return AFL_PLU;
|
|
case MEAN_NODE: return AFL_MNODE;
|
|
case TRUE_NODE: return AFL_TNODE;
|
|
case CHIRON: return AFL_CHI;
|
|
case LILITH: return AFL_LIL;
|
|
case CERES: return AFL_CER;
|
|
case PALLAS: return AFL_PAL;
|
|
case JUNO: return AFL_JUN;
|
|
case VESTA: return AFL_VES;
|
|
case AC: return AFL_AC;
|
|
case MC: return AFL_MC;
|
|
}
|
|
if (planet >= FIRST_HSNR && planet <= LAST_HSNR)
|
|
return AFL_HOUSE;
|
|
else
|
|
return -1;
|
|
}
|
|
|
|
/*
|
|
* get the 2-letter abbreviation for a planet
|
|
* returns ?? if not defined
|
|
*/
|
|
char *planet2abbr2(int planet)
|
|
{
|
|
switch (planet) {
|
|
case SUN: return "su";
|
|
case MOON: return "mo";
|
|
case MERCURY: return "me";
|
|
case VENUS: return "ve";
|
|
case MARS: return "ma";
|
|
case JUPITER: return "ju";
|
|
case SATURN: return "sa";
|
|
case URANUS: return "ur";
|
|
case NEPTUNE: return "ne";
|
|
case PLUTO: return "pl";
|
|
case MEAN_NODE: return "mn";
|
|
case TRUE_NODE: return "tn";
|
|
case CHIRON: return "ch";
|
|
case LILITH: return "li";
|
|
case CERES: return "ce";
|
|
case PALLAS: return "pa";
|
|
case JUNO: return "jn";
|
|
case VESTA: return "vs";
|
|
case AC: return "ac";
|
|
case MC: return "mc";
|
|
}
|
|
return "??";
|
|
}
|
|
|
|
/*
|
|
* get the 3-letter abbreviation for a planet
|
|
* returns ??? if not defined
|
|
*/
|
|
char *planet2abbr3(int planet)
|
|
{
|
|
switch (planet) {
|
|
case SUN: return "sun";
|
|
case MOON: return "mon";
|
|
case MERCURY: return "mer";
|
|
case VENUS: return "ven";
|
|
case MARS: return "mar";
|
|
case JUPITER: return "jup";
|
|
case SATURN: return "sat";
|
|
case URANUS: return "ura";
|
|
case NEPTUNE: return "nep";
|
|
case PLUTO: return "plu";
|
|
case MEAN_NODE: return "mno";
|
|
case TRUE_NODE: return "tno";
|
|
case CHIRON: return "chi";
|
|
case LILITH: return "lil";
|
|
case CERES: return "cer";
|
|
case PALLAS: return "pal";
|
|
case JUNO: return "jun";
|
|
case VESTA: return "ves";
|
|
case AC: return "asc";
|
|
case MC: return "mc ";
|
|
}
|
|
return "???";
|
|
}
|
|
|
|
char *placalc_set_ephepath(char *path)
|
|
{
|
|
static TLS char *epath;
|
|
if (path == NULL) return epath;
|
|
if (epath != NULL)
|
|
free((void *) epath);
|
|
epath = malloc(strlen(path) + 1);
|
|
if (epath != NULL) {
|
|
strcpy(epath, path);
|
|
swe_set_ephe_path(epath);
|
|
}
|
|
return epath;
|
|
}
|
|
|
|
char *placalc_get_errtext()
|
|
{
|
|
return perrtx;
|
|
}
|
|
|
|
/*****************************************************
|
|
deltat(t): returns delta t (in julian days) from universal time t
|
|
is included by users
|
|
ET = UT + deltat
|
|
******************************************************/
|
|
double deltat (double jd_ad) /* Astrodienst relative julian date */
|
|
{
|
|
return swe_deltat(jd_ad + JUL_OFFSET);
|
|
}
|
|
|
|
/******************************************************************/
|
|
double diff8360 (double a, double b)
|
|
/* a - b on a 360 degree circle, result -180..180*/
|
|
{
|
|
double d;
|
|
d = a - b;
|
|
if ( d >= 180.0 ) return( d - 360.0 );
|
|
if ( d < -180.0 ) return( d + 360.0 );
|
|
return( d );
|
|
} /* diff8360 */
|
|
|
|
/*
|
|
* originally in swephous.c
|
|
*/
|
|
|
|
/*************************************
|
|
return in which house pp is;
|
|
houses are numbered from 1 .. 12
|
|
*************************************/
|
|
int HouseNr(struct houses *hsp, CSEC pp)
|
|
{
|
|
CSEC cx;
|
|
int i = 2;
|
|
cx = swe_difcsn(pp, hsp->cusp [1]); /* distance from cusp 1 */
|
|
while (i < 13 && cx >= difcsn(hsp->cusp[i], hsp->cusp[1])) i++;
|
|
return (i - 1);
|
|
}
|
|
|
|
/************************************
|
|
returns the inp-house number, where pp is in
|
|
houses are numbered from 1 .. 12
|
|
************************************/
|
|
int InpHouseNr (struct houses *hsp, CSEC pp, CSEC *coff)
|
|
{
|
|
CSEC cx;
|
|
int i = 2;
|
|
cx = swe_difcsn(pp, hsp->cusp [1] + coff [1]);
|
|
while(i<13 && cx >= swe_difcsn(hsp->cusp[i] + coff[i], hsp->cusp[1] + coff[1]))
|
|
i++;
|
|
return (i - 1);
|
|
}
|
|
|
|
/* variation of InpHouseNr(). Able to handle house pre-orbs that are
|
|
* proportional to house size.
|
|
* value 1 in doff[0] means that the offset is proportional to house size,
|
|
* e.g. doff[ihs] = -5 means here:
|
|
* doff[ihs] = -5 / 30 * preceding_house_size;
|
|
* We first calculate the absolute offsets for each house of our birth chart,
|
|
* then call the function InpHouseNr() with those values.
|
|
*/
|
|
int InpHouseNr2 (struct houses *hsp, CSEC pp, CSEC *coff)
|
|
{
|
|
int i, j;
|
|
CSEC myoff[13];
|
|
for (i = 0; i < 13; i++)
|
|
myoff[i] = coff[i];
|
|
if (myoff[0] == 1) {
|
|
for (i = 1; i < 13; i++) {
|
|
j = i + 1;
|
|
if (j > 12) j = 1;
|
|
myoff[j] = swe_degnorm((hsp->cusp[j] - hsp->cusp[i]) / 360000.0) / 30.0 * myoff[j];
|
|
}
|
|
}
|
|
return InpHouseNr(hsp, pp, myoff);
|
|
}
|
|
|
|
// *********************************************************
|
|
// Houses:
|
|
// *********************************************************
|
|
// Koch and Placidus don't work in the polar circle.
|
|
// We swap MC/IC so that MC is always before AC in the zodiac
|
|
// We then divide the quadrants into 3 equal parts.
|
|
// *********************************************************
|
|
// All angles are expressed in centiseconds (1/100th of a
|
|
// second of arc) and integer arithmetic is used for these.
|
|
// *********************************************************
|
|
// Arguments: th = sidereal time (angle 0..360 degrees
|
|
// hsy = letter code for house system; implemen-
|
|
// ted are P,K,C,R,E,V.
|
|
// fi = geographic latitude
|
|
// ekl = obliquity of the ecliptic
|
|
// iteration_count = number of iterations in
|
|
// Placidus calculation; is obsolete but kept for API continuity
|
|
// *********************************************************
|
|
void CalcHouses(CSEC th, CSEC fi, CSEC mekl, char hsy, int iteration_count,
|
|
struct houses *hsp )
|
|
{
|
|
int retc = 0, i;
|
|
double cs2deg = 360000;
|
|
double cusps[13];
|
|
double ascmc[10];
|
|
/* iteration_count is always 2 */
|
|
retc = swe_houses_armc(th / cs2deg, fi / cs2deg, mekl / cs2deg, (int) hsy,
|
|
cusps, ascmc);
|
|
for (i = 0; i < 13; i++)
|
|
hsp->cusp[i] = swe_d2l(cusps[i] * cs2deg);
|
|
hsp->ac = swe_d2l(ascmc[0] * cs2deg);
|
|
hsp->mc = swe_d2l(ascmc[1] * cs2deg);
|
|
/*
|
|
* this is just to shut up lint
|
|
*/
|
|
retc += iteration_count;
|
|
iteration_count = retc;
|
|
} /* procedure houses */
|
|
|
|
/******************************/
|
|
void RecalcAspects(struct AspectType *a)
|
|
{
|
|
centisec diff,orbis;
|
|
int p1, p2, i;
|
|
struct aspRec *arp;
|
|
if (a->ppos2 == NULL) { /* no set ppos2, no interaspects */
|
|
for (p1 = 0; p1 < a->NrOfPlanets; p1++) {
|
|
a->Asp[p1][p1].index = 0; /* ignore p1 conjunct p1 */
|
|
for (p2 = p1 + 1; p2 < a->NrOfPlanets; p2++) {
|
|
arp = &(a->Asp[p1][p2]);
|
|
diff = a->PlanetPos [p2] - a->PlanetPos [p1];
|
|
if (diff >= DEG180)
|
|
diff -= DEG360;
|
|
else if (diff < - DEG180)
|
|
diff += DEG360;
|
|
i = 1;
|
|
arp->index = 0;
|
|
while (i <= a->NrOfAspects) {
|
|
orbis = ABS4 (diff) - a->Angle [i];
|
|
if (ABS4 (orbis) <= a->Maxorb [i]) {
|
|
arp->index = i;
|
|
arp->orb = orbis;
|
|
break; /* out of while */
|
|
}
|
|
i++;
|
|
}
|
|
a->Asp [p2][p1].index = arp->index;
|
|
a->Asp [p2][p1].orb = arp->orb;
|
|
} /* for p2 */
|
|
} /* for p1 */
|
|
} else { /* interaspects between set 1 and set 2 */
|
|
for (p1 = 0; p1 < a->NrOfPlanets; p1++) {
|
|
for (p2 = 0; p2 < a->NrOfPlanets; p2++) {
|
|
arp = &(a->Asp[p1][p2]);
|
|
diff = a->ppos2 [p2] - a->PlanetPos [p1];
|
|
if (diff >= DEG180)
|
|
diff -= DEG360;
|
|
else if (diff < - DEG180)
|
|
diff += DEG360;
|
|
i = 1;
|
|
arp->index = 0;
|
|
while (i <= a->NrOfAspects) {
|
|
orbis = ABS4 (diff) - a->Angle [i];
|
|
if (ABS4 (orbis) <= a->Maxorb [i]) {
|
|
arp->index = i;
|
|
arp->orb = orbis;
|
|
break; /* out of while */
|
|
}
|
|
i++;
|
|
}
|
|
} /* for p2 */
|
|
} /* for p1 */
|
|
} /* else */
|
|
a->dataValid = TRUE;
|
|
}
|
|
|
|
/***********************************************************
|
|
function sidtime (t): returns sidereal time at greenwich;
|
|
Parameters differ from ASYS version! after AESuppl. 1961, page 75
|
|
version 24-oct-87
|
|
***********************************************************/
|
|
double sidtime (double jd_ad, double ecl, double nuta)
|
|
/* jd_ad relative julian date */
|
|
/* ecl, nuta ecliptic and nutation of date, in degrees */
|
|
{
|
|
return swe_sidtime0(jd_ad + JUL_OFFSET, ecl, nuta);
|
|
}
|
|
|
|
# ifdef INTEL_BYTE_ORDER
|
|
/********************************************************************/
|
|
void longreorder (UCHAR *p, int n)
|
|
/* p points to memory filled with int32 values; for
|
|
each of the values the seqeuence of the four bytes
|
|
has to be reversed, to translate HP-UX and VAX
|
|
ordering to MSDOS/Turboc ordering */
|
|
{
|
|
int i;
|
|
unsigned char c0, c1, c2, c3;
|
|
for (i = 0; i < n; i += 4, p += 4) {
|
|
c0 = *p;
|
|
c1 = *(p + 1);
|
|
c2 = *(p + 2);
|
|
c3 = *(p + 3);
|
|
*p = c3;
|
|
*(p + 1) = c2;
|
|
*(p + 2) = c1;
|
|
*(p + 3) = c0;
|
|
}
|
|
}
|
|
# endif
|