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

#include "placalc.h"

#ifdef PLACALC
#ifdef ASTROLOG
/* Begin contents of placalc.c */
#endif
/*****************************************************
$Header: placalc.c,v 1.14 93/07/19 22:13:07 alois Exp $


  ---------------------------------------------------------------
  | Copyright Astrodienst AG and Alois Treindl, 1989,1991,1993  |
  | The use of this source code is subject to regulations made	|
  | by Astrodienst Zurich. The code is NOT in the public domain.|
  |								|
  | This copyright notice must not be changed or removed	|
  | by any user of this program.				|
  ---------------------------------------------------------------

Important changes:
11-jun-93 revision 1.12: fixed error which affected Mercury between -2100 and
			-3100 (it jumped wildly).

*******************************************************/

#ifndef ASTROLOG
#include "placalc.h"	/* includes ourdef.h */
#include "astrolib.h"	/* includes ourdef.h */
#include "helconst.c"	/* orbital elements and disturbations
			for inner planets and moon */
#include "deltat.c"
#else /* ASTROLOG */
#ifdef ASTROLOG
/* Begin contents of helconst.c */
#endif
/***********************************************************
 * $Header$
 
 definition module for planetary elements
 and disturbation coefficients
 version HP-UX C  for new version with stored outer planets
 31-jul-88
 by Alois Treindl 

  ---------------------------------------------------------------
  | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.	|
  | The use of this source code is subject to regulations made	|
  | by Astrodienst Zurich. The code is NOT in the public domain.|
  |								|
  | This copyright notice must not be changed or removed	|
  | by any user of this program.				|
  ---------------------------------------------------------------

***********************************************************/
/* In the elements degrees were kept as the units for the constants. This
requires conversion to radians, when the actual calculations are performed.
This approach is not the most efficient, but safer for development.
Constant conversion could be done by writing all degree constants with
 value * DEGTORAD */

# define TIDAL_26	TRUE	/* decide wheter to use new or old lunar tidal
				term; a consistent system of delta t must be 
				used */
# define MOON_TEST_CORR FALSE	/* to include more lunar terms in longitude */

REAL8  ekld [4] = { 23.452294, -46.845, -.0059, 0.00181 };
	/* ecliptic with epoch1900, Ekd(0..3) in basic */
struct eledata {
  REAL8 axis,		/* mean distance, in a.u., A(N) in basic */
	period,		/* days for one revolution, P(N) in basic */
	epoch,		/* relative juldate of epoch, Ep(N) in basic */
			/* T = distance to epoch in jul.centuries 36525 day*/
	lg0,lg1,lg2,lg3,/* deg(epoch), degree/day, seconds/T^2,seconds/T^3 */
			/* Pd(N,0..2) in basic, lg3 was not present */
	pe0,pe1,pe2,pe3,/* deg(epoch), seconds/T,  seconds/T^2,seconds/T^3 */
			/* Pd(N,3..5) in basic, pe3 was not present */
	ex0,ex1,ex2,	/* ecl(epoch), 1/T, 1/T^2	*/
			/* Pd(N,6..8) in basic */
	kn0,kn1,kn2,kn3,/* node(epoch),seconds/T,  seconds/T^2,seconds/T^3 */
			/* Pd(N,9..11) in basic, kn3 was not present */
	in0,in1,in2;    /* incl(epoch),1/T, 1/T^2	*/
			/* Pd(N,12..14) in basic */
	} pd [MARS + 1] = {
	{/*earth*/	1.00000023,	365.25636042,	EPOCH1900,
	  99.696678,	.9856473354,	1.089,		0,	
	  101.220833,	6189.03,	1.63,		0.012,
	  0.01675104,	-0.00004180,	-0.000000126,
	  0,		0,		0,		0,
	  0,		0,		0},
	  /* note 29 June 88 by Alois: G.M.Clemence, Astronomical Journal
	  vol.53,p. 178 (1948) gives a correction to the perihel motion
	  of -4.78" T, giving 6184.25 for the linear Term above. We have
	  not yet applied this correction. It has been used in APAE 22,4
	  on the motion of mars and does make an official impression. */
	{/*moon*/	0.0025955307,	27.321661,	EPOCH1900,
# if ! TIDAL_26 
	/* values from Improved Lunar Ephemeris, corresponding to tidal
	   term -22.44"/cy and  consistent with delta t ~ 29.949 T*T
	*/
	  270.4341638,	13.176396526808121, -4.08,	0.0068,
# endif
# if TIDAL_26 
	/* new values from Morrison 1979, with tidal term -26"/cy as
	   stated in A.E. 1986 onwards, consistent with delta t ~ 44.3 T*T 
	   correction: -1.54" + 2.33" T - 1.78" T*T
	*/
	  270.4337361,	13.176396544528099, -5.86,	0.0068,
# endif
	  334.329556,	14648522.52,	-37.17,		-0.045,
	  0.054900489,	0,		0,
	  259.183275,	-6962911.23,	7.48 ,		0.008,
	  5.145388889,	0,		0},
	{/*mercury*/	.3870986,	87.969252,	EPOCH1900,
	  178.179078,	4.0923770233,	1.084,		0,
	  75.89969722,	5599.76,	1.061,		0,
	  0.20561421,	.00002046,	 -.000000030,
	  47.145944444,	4266.75,	.626,		0,
	  7.0028805555,	6.699,		-.066},
	{/*venus*/	.72333162,	224.700726,	EPOCH1900,
	  342.767053,	1.6021687039,	1.1148,		0,
	  130.16383333,	5068.93,	-3.515,		0,
	  0.00682069,	-.00004774,	.000000091,
	  75.7796472223,3239.46,	1.476,		0,
	  3.3936305555,	3.621,		.0035},
	{/*mars*/	1.5236914620,	686.9296097,	EPOCH1900,
	  /* These are the corrected elements by Ross */
	  293.74762778,	.524071163814,	1.1184,		0,
	  334.21820278,	6626.73,	.4675,		-0.0043,
	  0.09331290,	.000092064,	-.000000077,
	  48.786441667,	2775.57,	-.005,		-0.0192,
	  1.85033333,	-2.430,		.0454}
};

/*
 * mimimum and maximum distances computed over 1000 years with plamimax,
 * required for relative distances rgeo, where the distance is given
 * as 100 when a planet is closest and as 0 when farthest from earth.
 */
REAL8 rmima[CALC_N][2] = {	
	{ 0.98296342,  1.01704665},
	{ 0.00238267,  0.00271861},
	{ 0.54900496,  1.45169607},
	{ 0.26411287,  1.73597885},
	{ 0.37289847,  2.67626927},
	{ 3.94877993,  6.45627627},
	{ 7.99362824, 11.09276636},
	{17.28622633, 21.10714104},
	{28.81374786, 31.33507284},
	{28.67716748, 50.29208774},
	{ 0.00,	 0.00259553},	/* nodes don't get a real value*/
	{ 0.00,  0.00259553},
	{ 7.36277475, 19.86585062}};
      
	

#define SDNUM 20
struct sdat {		/* 0..19 mean anomalies of disturbing planets 
			Sd(0..19,0..1) in basic */
	REAL8 sd0,	/* mean anomaly at epoch 1850 */
	      sd1;	/* degrees/year */
	} sd [SDNUM] = {
	 114.50,	585.17493,
	 109.856,	191.39977,
	 148.031,	30.34583,
	 284.716,	12.21794,
	 114.508,	585.17656,
	 -0.56,		359.99213,
	 148.03,	30.34743,
	 284.72,	12.2196,
	 248.07,	1494.726615,
	 359.44,	359.993595,
	 109.86,	191.402867,
	 148.02,	30.348930,
	 114.503,	585.173715,
	 359.444,	359.989285,
	 148.021,	30.344620,
	 284.716,	12.21669,
	 148.0315,	30.34906264,
	 284.7158,	12.22117085,
	 220.1695, 	4.284931111,
	 291.8024, 	2.184704167
};

REAL8  sa [SDNUM];

struct kor {
	int     j, i;
	REAL8   lampl;	/* amplitude of disturbation in long, seconds of arc */
	REAL8	lphase; /* phase of disturbation in long, degrees */
	INT4    rampl;  /* ampl. of disturbation in radius, 9th place of log */
	REAL8   rphase; /* phase of disturbation in radius, degrees */
	int     k;	/* index into disturbing planet anomaly table sa[] */
};
       /* delta long = lampl * COS (lphase - arg) in seconds of arc
	  delta rad  = rampl * COS (rphase - arg) in ninth place of log
	  arg = j * sa (k) + i * ma (this planet)
	  ma = mean anomaly
	  sa = mean anomaly of disturbing planet, where this
	       is taken from the aproximate value in sa[]
	  For the COS (phase - arg) it is good enough to compute
	  with 32 bit reals, because ampl and phase have only 
	  four to five significant digits. 
	  While saving constant space, it is costing execution time due
	  to float/double conversions.
	*/
        /* In basic, all correction terms for sun, mercury, venus and mars
	   were contained in one array K(0..142,0..6); Nk(N,0) contained
	   the index of the first term of planet N and Nk(N,1) the number
	   of terms for this planet. Here, we use a  0 in the first column
	   kor.j to indicate the end of the table for a planet.
	   K(*) was a basic INTEGER array, therefore the amplitudes and phases
	   had to be expressed as
	   K(i,2) = ampl. of longitude in 0.001 seconds of arc
	   K(i,3) = phase of longitude in 0.01 degrees
	   K(i,4) = ampl. of radius in 9th place of log
	   K(i,5) = phase of radius in 0.01 degrees.
	   Here we have converted the amplitude of long. to seconds of arc
	   and the phases to degrees.
	 */

struct kor earthkor[] = {	/* 11-jul-88 all terms to 0.020" longitude */
      /*  j     i       lampl   lphase  rampl   rphase  k */
	 -1,	1,	0.013,	243,	28,	335,	8,	/* mercury */
	 -1,	3,	0.015,	357,	18,	267,	8,
	 -1,	4,	0.023,	326,	5,	239,	8,
	 -1,	0,	0.075,	296.6,	94,	205.0,	0,	/* Venus */
	 -1,	1,	4.838,	299.10,	2359,	209.08,	0,	
	 -1,	2,	0.074,	207.9,	69,	348.5,	0,
	 -1,	3,	0.009,	249,	16,	330,	0,
   	 -2,	1,	.116,	148.90,	160,	58.40,	0,	
	 -2,	2,	5.526,	148.31,	6842,	58.32,	0,	
	 -2,	3,	2.497,	315.94,	869,	226.70,	0,	
	 -2,	4,	0.044,	311.4,	52,	38.8,	0,
	 -3,	2,	0.013,	176,	21,	90,	0,
	 -3,	3,	.666,	177.71,	1045,	87.57,	0,	
	 -3,	4,	1.559,	-14.75,	1497,	255.25,	0,	
	 -3,	5,	1.024,	318.15,	194,	49.50,	0,	
	 -3,	6,	0.017,	315,	19,	43,	0,
	 -4,	4,	.210,	206.20,	376,	116.28,	0,	
	 -4,	5,	.144,	195.40,	196,	105.20,	0,	
	 -4,	6,	.152,	-16.20,	94,	254.80,	0,	
	 -5,	5,	0.084,	235.6,	163,	145.4,	0,
	 -5,	6,	0.037,	221.8,	59,	132.2,	0,
	 -5,	7,	.123,	195.30,	141,	105.40,	0,	
	 -5,	8,	.154,	-.40,	26,	270.00,	0,	
	 -6,	6,	0.038,	264.1,	80,	174.3,	0,
	 -6,	7,	0.014,	253,	25,	164,	0,
	 -6,	8,	0.01,	230,	14,	135,	0,
	 -6,	9,	0.014,	12,	12,	284,	0,
	 -7,	7,	0.020,	294,	42,	203.5,	0,
	 -7,	8,	0.006,	279,	12,	194,	0,
	 -8,	8,	0.011,	322,	24,	234,	0,
	 -8,	12,	0.042,	259.2,	44,	169.7,	0,
	 -8,	14,	0.032,	48.8,	33,	138.7,	0,
	 -9,	9,	0.006,	351,	13,	261,	0,
	 1,	-1,	.273,	217.70,	150,	127.70,	1,	/* mars */
	 1,	0,	0.048,	260.3,	28,	347,	1,
	 2,	-3,	0.041,	346,	52,	255.4,	1,
	 2,	-2,	2.043,	343.89,	2057,	253.83,	1,
	 2,	-1,	1.770,	200.40,	151,	295.00,	1,
	 2,	0,	0.028,	148,	31,	234.3,	1,  
	 3,	-3,	.129,	294.20,	168,	203.50,	1,
	 3,	-2,	.425,	-21.12,	215,	249.00,	1,
	 4,	-4,	0.034,	71,	49,	339.7,	1,
	 4,	-3,	.500,	105.18,	478,	15.17,	1,
	 4,	-2,	.585,	-25.94,	105,	65.90,	1,
	 5,	-4,	0.085,	54.6,	107,	324.6,	1,
	 5,	-3,	.204,	100.80,	89,	11.00,	1,
	 6,	-5,	0.020,	186,	30,	95.7,	1,
	 6,	-4,	.154,	227.40,	139,	137.30,	1,
	 6,	-3,	.101,	96.30,	27,	188.00,	1,
	 7,	-5,	0.049,	176.5,	60,	86.2,	1,
	 7,	-4,	.106,	222.70,	38,	132.90,	1,
	 8,	-5,	0.052,	348.9,	45,	259.7,	1,
	 8,	-4,	0.021,	215.2,	8,	310,	1,
	 8,	-6,	0.010,	307,	15,	217,	1,
	 9,	-6,	0.028,	298,	34,	208.1,	1,
	 9,	-5,	0.062,	346,	17,	257,	1,
	 10,	-6,	0.019,	111,	15,	23,	1,
	 11,	-7,	0.017,	59,	20,	330,	1,
	 11,	-6,	0.044,	105.9,	9,	21,	1,
	 13,	-8,	0.013,	184,	15,	94,	1,
	 13,	-7,	0.045,	227.8,	5,	143,	1,
	 15,	-9,	0.021,	309,	22,	220,	1,
	 17,	-9,	0.026,	113,	0,	0,	1,
	 1,	-2,	.163,	198.60,	208,	112.00,	2,	/* jupiter */
	 1,	-1,	7.208,	179.53,	7067,	89.55,	2,
	 1,	0,	2.600,	263.22,	244,	-21.40,	2,
	 1,	1,	0.073,	276.3,	80,	6.5,	2,
	 2,	-3,	0.069,	80.8,	103,	350.5,	2,
	 2,	-2,	2.731,	87.15,	4026,	-2.89,	2,
	 2,	-1,	1.610,	109.49,	1459,	19.47,	2,
	 2,	0,	0.073,	252.6,	8,	263,	2,
	 3,	-3,	.164,	170.50,	281,	81.20,	2,
	 3,	-2,	.556,	82.65,	803,	-7.44,	2,
	 3,	-1,	.210,	98.50,	174,	8.60,	2,
	 4,	-4,	0.016,	259,	29,	170,	2,
	 4,	-3,	0.044,	168.2,	74,	79.9,	2,
	 4,	-2,	0.080,	77.7,	113,	347.7,	2,
	 4,	-1,	0.023,	93,	17,	3,	2,
	 5,	-2,	0.009,	71,	14,	343,	2,
	 1,	-2,	0.011,	105,	15,	11,	3,	/* saturn */
	 1,	-1,	.419,	100.58,	429,	10.60,	3,
	 1,	0,	.320,	269.46,	8,	-7.00,	3,
	 2,	-2,	.108,	290.60,	162,	200.60,	3,
	 2,	-1,	.112,	293.60,	112,	203.10,	3,
	 3,	-2,	0.021,	289,	32,	200.1,	3,
	 3,	-1,	0.017,	291,	17,	201,	3,
	 ENDMARK
};

struct kor mercurykor[] = {