/*
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 *
 *  Copyright (C) 2000 - 2005 Liam Girdwood
 */

module nova.pluto;

import std.math;

import nova.vsop87;
import nova.solar;
import nova.earth;
import nova.transform;
import nova.rise_set;
import nova.utility;

enum PLUTO_COEFFS = 43;

struct pluto_argument {
	double J, S, P;
}

struct pluto_longitude {
	double A, B;
}

struct pluto_latitude {
	double A, B;
}

struct pluto_radius {
	double A, B;
}

/* cache variables */
static double cJD = 0.0, cL = 0.0, cB = 0.0, cR = 0.0;

static const pluto_argument[PLUTO_COEFFS] argument = [
	{0, 0, 1},
	{0, 0, 2},
	{0, 0, 3},
	{0, 0, 4},
	{0, 0, 5},
	{0, 0, 6},
	{0, 1, -1},
	{0, 1, 0},
	{0, 1, 1},
	{0, 1, 2},
	{0, 1, 3},
	{0, 2, -2},
	{0, 2, -1},
	{0, 2, 0},
	{1, -1, 0},
	{1, -1, 1},
	{1, 0, -3},
	{1, 0, -2},
	{1, 0, -1},
	{1, 0, 0},
	{1, 0, 1},
	{1, 0, 2},
	{1, 0, 3},
	{1, 0, 4},
	{1, 1, -3},
	{1, 1, -2},
	{1, 1, -1},
	{1, 1, 0},
	{1, 1, 1},
	{1, 1, 3},
	{2, 0, -6},
	{2, 0, -5},
	{2, 0, -4},
	{2, 0, -3},
	{2, 0, -2},
	{2, 0, -1},
	{2, 0, 0},
	{2, 0, 1},
	{2, 0, 2},
	{2, 0, 3},
	{3, 0, -2},
	{3, 0, -1},
	{3, 0, 0}
];


static const pluto_longitude[PLUTO_COEFFS] longitude = [
	{-19799805, 19850055},
	{897144, -4954829},
	{611149, 1211027},
	{-341243, -189585},
	{129287, -34992},
	{-38164, 30893},
	{20442, -9987},
	{-4063, -5071},
	{-6016, -3336},
	{-3956, 3039},
	{-667, 3572},
	{1276, 501},
	{1152, -917},
	{630, -1277},
	{2571, -459},
	{899, -1449},
	{-1016, 1043},
	{-2343, -1012},
	{7042, 788},
	{1199, -338},
	{418, -67},
	{120, -274},
	{-60, -159},
	{-82, -29},
	{-36, -20},
	{-40, 7},
	{-14, 22},
	{4, 13},
	{5,2},
	{-1,0},
	{2,0},
	{-4, 5},
	{4, -7},
	{14, 24},
	{-49, -34},
	{163, -48},
	{9, 24},
	{-4, 1},
	{-3,1},
	{1,3},
	{-3, -1},
	{5, -3},
	{0,0}
];

static const pluto_latitude[PLUTO_COEFFS] latitude = [
	{-5452852, -14974862},
	{3527812, 1672790},
	{-1050748, 327647},
	{178690, -292153},
	{18650, 100340},
	{-30697, -25823},
	{4878, 11248},
	{226, -64},
	{2030, -836},
	{69, -604},
	{-247, -567},
	{-57, 1},
	{-122, 175},
	{-49, -164},
	{-197, 199},
	{-25, 217},
	{589, -248},
	{-269, 711},
	{185, 193},
	{315, 807},
	{-130, -43},
	{5, 3},
	{2, 17},
	{2, 5},
	{2, 3},
	{3, 1},
	{2, -1},
	{1, -1},
	{0, -1},
	{0, 0},
	{0, -2},
	{2, 2},
	{-7, 0},
	{10, -8},
	{-3, 20},
	{6, 5},
	{14, 17},
	{-2, 0},
	{0, 0},
	{0, 0},
	{0, 1},
	{0, 0},
	{1, 0}
];

static const pluto_radius[PLUTO_COEFFS] radius = [
	{66865439, 68951812},
	{-11827535, -332538},
	{1593179, -1438890},
	{-18444, 483220},
	{-65977, -85431},
	{31174, -6032},
	{-5794, 22161},
	{4601, 4032},
	{-1729, 234},
	{-415, 702},
	{239, 723},
	{67, -67},
	{1034, -451},
	{-129, 504},
	{480, -231},
	{2, -441},
	{-3359, 265},
	{7856, -7832},
	{36, 45763},
	{8663, 8547},
	{-809, -769},
	{263, -144},
	{-126, 32},
	{-35, -16},
	{-19, -4},
	{-15, 8},
	{-4, 12},
	{5, 6},
	{3, 1},
	{6, -2},
	{2, 2},
	{-2, -2},
	{14, 13},
	{-63, 13},
	{136, -236},
	{273, 1065},
	{251, 149},
	{-25, -9},
	{9, -2},
	{-8, 7},
	{2, -10},
	{19, 35},
	{10, 2}
];

extern (C) {

/*! \fn void ln_get_pluto_equ_coords(double JD, struct ln_equ_posn *position);
* \param JD julian Day
* \param position Pointer to store position
*
* Calculates Pluto's equatorial position for the given julian day.
*/
@nogc void ln_get_pluto_equ_coords(double JD, ref ln_equ_posn position) nothrow
{
	ln_helio_posn h_sol, h_pluto;
	ln_rect_posn g_sol, g_pluto;
	double a,b,c;
	double ra, dec, delta, diff, last, t = 0;

	/* need typdef for solar heliocentric coords */
	ln_get_solar_geom_coords(JD, h_sol);
	ln_get_rect_from_helio(h_sol, g_sol);

	do {
		last = t;
		ln_get_pluto_helio_coords(JD - t, h_pluto);
		ln_get_rect_from_helio(h_pluto, g_pluto);

		/* equ 33.10 pg 229 */
		a = g_sol.X + g_pluto.X;
		b = g_sol.Y + g_pluto.Y;
		c = g_sol.Z + g_pluto.Z;

		delta = a * a + b * b + c * c;
		delta = sqrt(delta);
		t = delta * 0.0057755183;
		diff = t - last;
	} while (diff > 0.0001 || diff < -0.0001);

	ra = atan2(b, a);
	dec = c / delta;
	dec = asin(dec);

	/* back to hours, degrees */
	position.ra = ln_range_degrees(ln_rad_to_deg(ra));
	position.dec = ln_rad_to_deg(dec);
}


/*! \fn void ln_get_pluto_helio_coords(double JD, struct ln_helio_posn *position)
* \param JD Julian Day
* \param position Pointer to store new heliocentric position
*
* Calculate Pluto's heliocentric coordinates for the given julian day.
* This function is accurate to within 0.07" in longitude, 0.02" in latitude
* and 0.000006 AU in radius vector.
*
* Note: This function is not valid outside the period of 1885-2099.
*/
/* Chap 37. Equ 37.1
*/

@nogc void ln_get_pluto_helio_coords(double JD, ref ln_helio_posn position) nothrow
{
	double sum_longitude = 0, sum_latitude = 0, sum_radius = 0;
	double J, S, P;
	double t, a, sin_a, cos_a;
	int i;

	/* check cache first */
	if(JD == cJD) {
		/* cache hit */
		position.L = cL;
		position.B = cB;
		position.R = cR;
		return;
	}

	/* get julian centuries since J2000 */
	t =(JD - 2451545.0) / 36525.0;

	/* calculate mean longitudes for jupiter, saturn and pluto */
	J =  34.35 + 3034.9057 * t;
   	S =  50.08 + 1222.1138 * t;
   	P = 238.96 +  144.9600 * t;

	/* calc periodic terms in table 37.A */
	for (i = 0; i < PLUTO_COEFFS; i++) {
		a = argument[i].J * J + argument[i].S * S + argument[i].P * P;
		sin_a = sin(ln_deg_to_rad(a));
		cos_a = cos(ln_deg_to_rad(a));

		/* longitude */
		sum_longitude += longitude[i].A * sin_a + longitude[i].B * cos_a;

		/* latitude */
		sum_latitude += latitude[i].A * sin_a + latitude[i].B * cos_a;

		/* radius */
		sum_radius += radius[i].A * sin_a + radius[i].B * cos_a;
	}

	/* calc L, B, R */
	position.L = 238.958116 + 144.96 * t + sum_longitude * 0.000001;
	position.B = -3.908239 + sum_latitude * 0.000001;
	position.R = 40.7241346 + sum_radius * 0.0000001;

	/* save cache */
	cJD = JD;
	cL = position.L;
	cB = position.B;
	cR = position.R;
}

/*! \fn double ln_get_pluto_earth_dist(double JD);
* \param JD Julian day
* \return Distance in AU
*
* Calculates the distance in AU between the Earth and Pluto for the
* given julian day.
*/
@nogc double ln_get_pluto_earth_dist(double JD) nothrow
{
	ln_helio_posn h_pluto, h_earth;
	ln_rect_posn g_pluto, g_earth;
	double x, y, z;

	/* get heliocentric positions */
	ln_get_pluto_helio_coords(JD, h_pluto);
	ln_get_earth_helio_coords(JD, h_earth);

	/* get geocentric coords */
	ln_get_rect_from_helio(h_pluto, g_pluto);
	ln_get_rect_from_helio(h_earth, g_earth);

	/* use pythag */
	x = g_pluto.X - g_earth.X;
	y = g_pluto.Y - g_earth.Y;
	z = g_pluto.Z - g_earth.Z;
	x = x * x;
	y = y * y;
	z = z * z;

	return sqrt(x + y + z);
}

/*! \fn double ln_get_pluto_solar_dist(double JD);
* \param JD Julian day
* \return Distance in AU
*
* Calculates the distance in AU between the Sun and Pluto for the
* given julian day.
*/
@nogc double ln_get_pluto_solar_dist(double JD) nothrow
{
	ln_helio_posn h_pluto;

	/* get heliocentric position */
	ln_get_pluto_helio_coords(JD, h_pluto);
	return h_pluto.R;
}

/*! \fn double ln_get_pluto_magnitude(double JD);
* \param JD Julian day
* \return Visible magnitude of Pluto
*
* Calculate the visible magnitude of Pluto for the given
* julian day.
*/
@nogc double ln_get_pluto_magnitude(double JD) nothrow
{
	double delta, r;

	/* get distances */
	r = ln_get_pluto_solar_dist(JD);
	delta = ln_get_pluto_earth_dist(JD);

	return -1.0 + 5.0 * log10(r * delta);
}

/*! \fn double ln_get_pluto_disk(double JD);
* \param JD Julian day
* \return Illuminated fraction of Plutos disk
*
* Calculate the illuminated fraction of Pluto's disk for
* the given julian day.
*/
/* Chapter 41 */
@nogc double ln_get_pluto_disk(double JD) nothrow
{
	double r,delta,R;

	/* get distances */
	R = ln_get_earth_solar_dist(JD);
	r = ln_get_pluto_solar_dist(JD);
	delta = ln_get_pluto_earth_dist(JD);

	/* calc fraction angle */
	return (((r + delta) * (r + delta)) - R * R) / (4.0 * r * delta);
}

/*! \fn double ln_get_pluto_phase(double JD);
* \param JD Julian day
* \return Phase angle of Pluto (degrees)
*
* Calculate the phase angle of Pluto (Sun - Pluto - Earth)
* for the given julian day.
*/
/* Chapter 41 */
@nogc double ln_get_pluto_phase(double JD) nothrow
{
	double i,r,delta,R;

	/* get distances */
	R = ln_get_earth_solar_dist(JD);
	r = ln_get_pluto_solar_dist(JD);
	delta = ln_get_pluto_earth_dist(JD);

	/* calc phase */
	i = (r * r + delta * delta - R * R) / (2.0 * r * delta);
	i = acos(i);
	return ln_rad_to_deg(i);
}


/*! \fn double ln_get_pluto_rst(double JD, struct ln_lnlat_posn *observer, struct ln_rst_time *rst);
* \param JD Julian day
* \param observer Observers position
* \param rst Pointer to store Rise, Set and Transit time in JD
* \return 0 for success, else 1 for circumpolar.
*
* Calculate the time the rise, set and transit (crosses the local meridian at upper culmination)
* time of Pluto for the given Julian day.
*
* Note: this functions returns 1 if Pluto is circumpolar, that is it remains the whole
* day either above or below the horizon.
*/
@nogc int ln_get_pluto_rst(double JD, const ref ln_lnlat_posn observer,
		ref ln_rst_time rst) nothrow
{
	return ln_get_body_rst_horizon(JD, observer, &ln_get_pluto_equ_coords,
		LN_STAR_STANDART_HORIZON, rst);
}


/*! \fn double ln_get_pluto_sdiam(double JD)
* \param JD Julian day
* \return Semidiameter in arc seconds
*
* Calculate the semidiameter of Pluto in arc seconds for the
* given julian day.
*/
@nogc double ln_get_pluto_sdiam(double JD) nothrow
{
	double So = 2.07; /* at 1 AU */
	double dist;

	dist = ln_get_pluto_earth_dist(JD);
	return So / dist;
}

/*! \fn void ln_get_pluto_rect_helio(double JD, struct ln_rect_posn *position)
* \param JD Julian day.
* \param position pointer to return position
*
* Calculate Plutos rectangular heliocentric coordinates for the
* given Julian day. Coordinates are in AU.
*/
@nogc void ln_get_pluto_rect_helio(double JD, ref ln_rect_posn position) nothrow
{
	ln_helio_posn pluto;

	ln_get_pluto_helio_coords(JD, pluto);
	ln_get_rect_from_helio(pluto, position);
}

}