/*
 *  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 <lgirdwood@gmail.com>
 */

module nova.transform;

import std.math;

import nova.utility;
import nova.sidereal_time;
import nova.nutation;
import nova.precession;
import nova.ln_types;

extern (C) {

/*! \fn void ln_get_rect_from_helio(struct ln_helio_posn *object, struct ln_rect_posn *position);
* \param object Object heliocentric coordinates
* \param position Pointer to store new position
*
* Transform an objects heliocentric ecliptical coordinates
* into heliocentric rectangular coordinates.
*/
/* Equ 37.1 Pg 264
*/
@nogc void ln_get_rect_from_helio(const ref ln_helio_posn object,
	ref ln_rect_posn position) nothrow
{
	double sin_e, cos_e;
	double cos_B, sin_B, sin_L, cos_L;

	/* ecliptic J2000 */
	sin_e = 0.397777156;
	cos_e = 0.917482062;

	/* calc common values */
	cos_B = cos(ln_deg_to_rad(object.B));
	cos_L = cos(ln_deg_to_rad(object.L));
	sin_B = sin(ln_deg_to_rad(object.B));
	sin_L = sin(ln_deg_to_rad(object.L));

	/* equ 37.1 */
	position.X = object.R * cos_L * cos_B;
	position.Y = object.R * (sin_L * cos_B * cos_e - sin_B * sin_e);
	position.Z = object.R * (sin_L * cos_B * sin_e + sin_B * cos_e);
}

/*! \fn void ln_get_hrz_from_equ(struct ln_equ_posn *object, struct ln_lnlat_posn *observer, double JD, struct ln_hrz_posn *position)
* \param object Object coordinates.
* \param observer Observer cordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects equatorial coordinates into horizontal coordinates
* for the given julian day and observers position.
*
* 0 deg azimuth = south, 90 deg = west.
*/
/* Equ 12.1,12.2 pg 88
*
* TODO:
* Transform horizontal coordinates to galactic coordinates.
*/

@nogc void ln_get_hrz_from_equ(const ref ln_equ_posn object,
	const ref ln_lnlat_posn observer, double JD, ref ln_hrz_posn position) nothrow
{
	double sidereal;

	/* get mean sidereal time in hours*/
	sidereal = ln_get_mean_sidereal_time(JD);
	ln_get_hrz_from_equ_sidereal_time (object, observer, sidereal, position);
}


@nogc void ln_get_hrz_from_equ_sidereal_time(const ref ln_equ_posn object,
	const ref ln_lnlat_posn observer, double sidereal,
	ref ln_hrz_posn position) nothrow
{
	real H, ra, latitude, declination, A, Ac, As, h, Z, Zs;

	/* change sidereal_time from hours to radians*/
	sidereal *= 2.0 * PI / 24.0;

	/* calculate hour angle of object at observers position */
	ra = ln_deg_to_rad(object.ra);
	H = sidereal + ln_deg_to_rad(observer.lng) - ra;

	/* hence formula 12.5 and 12.6 give */
	/* convert to radians - hour angle, observers latitude, object declination */
	latitude = ln_deg_to_rad(observer.lat);
	declination = ln_deg_to_rad(object.dec);

	/* formula 12.6 *; missuse of A (you have been warned) */
	A = sin(latitude) * sin(declination) + cos(latitude) *
		cos(declination) * cos(H);
	h = asin(A);

	/* convert back to degrees */
	position.alt = ln_rad_to_deg(h);

	/* zenith distance, Telescope Control 6.8a */
	Z = acos(A);

	/* is'n there better way to compute that? */
	Zs = sin(Z);

	/* sane check for zenith distance; don't try to divide by 0 */
	if (fabs(Zs) < 1e-5) {
		if (object.dec > 0.0)
			position.az = 180.0;
		else
			position.az = 0.0;
		if ((object.dec > 0.0 && observer.lat > 0.0)
		   || (object.dec < 0.0 && observer.lat < 0.0))
		  	position.alt = 90.0;
		else
		  	position.alt = -90.0;
		return;
	}

	/* formulas TC 6.8d Taff 1991, pp. 2 and 13 - vector transformations */
	As = (cos(declination) * sin(H)) / Zs;
	Ac = (sin(latitude) * cos(declination) * cos(H) - cos(latitude) *
			sin(declination)) / Zs;

	// don't blom at atan2
	if (Ac == 0.0 && As == 0.0) {
	        if (object.dec > 0)
			position.az = 180.0;
		else
			position.az = 0.0;
		return;
	}
	A = atan2(As, Ac);

	/* convert back to degrees */
	position.az = ln_range_degrees(ln_rad_to_deg(A));
}

/*! \fn void ln_get_equ_from_hrz(struct ln_hrz_posn *object, struct ln_lnlat_posn *observer, double JD, struct ln_equ_posn *position)
* \param object Object coordinates.
* \param observer Observer cordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects horizontal coordinates into equatorial coordinates
* for the given julian day and observers position.
*/
@nogc void ln_get_equ_from_hrz(const ref ln_hrz_posn object,
	const ref ln_lnlat_posn observer, double JD, ref ln_equ_posn position) nothrow
{
	real H, longitude, declination, latitude, A, h, sidereal;

	/* change observer/object position into radians */

	/* object alt/az */
	A = ln_deg_to_rad(object.az);
	h = ln_deg_to_rad(object.alt);

	/* observer long / lat */
	longitude = ln_deg_to_rad(observer.lng);
	latitude = ln_deg_to_rad(observer.lat);

	/* equ on pg89 */
	H = atan2(sin(A), (cos(A) * sin(latitude) + tan(h) * cos(latitude)));
	declination = sin(latitude) * sin(h) - cos(latitude) * cos(h) * cos(A);
	declination = asin(declination);

	/* get ra = sidereal - longitude + H and change sidereal to radians*/
	sidereal = ln_get_apparent_sidereal_time(JD);
	sidereal *= 2.0 * PI / 24.0;

	position.ra = ln_range_degrees(ln_rad_to_deg(sidereal - H + longitude));
	position.dec = ln_rad_to_deg(declination);
}

/*! \fn void ln_get_equ_from_ecl(struct ln_lnlat_posn *object, double JD, struct ln_equ_posn *position)
* \param object Object coordinates.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects ecliptical coordinates into equatorial coordinates
* for the given julian day.
*/
/* Equ 12.3, 12.4 pg 89
*/
@nogc void ln_get_equ_from_ecl(const ref ln_lnlat_posn object, double JD,
	ref ln_equ_posn position) nothrow
{
	double ra, declination, longitude, latitude;
	ln_nutation nutation;

	/* get obliquity of ecliptic and change it to rads */
	ln_get_nutation(JD, nutation);
	nutation.ecliptic = ln_deg_to_rad(nutation.ecliptic);

	/* change object's position into radians */

	/* object */
	longitude = ln_deg_to_rad(object.lng);
	latitude = ln_deg_to_rad(object.lat);

	/* Equ 12.3, 12.4 */
	ra = atan2((sin(longitude) * cos(nutation.ecliptic) -
		tan(latitude) * sin(nutation.ecliptic)), cos(longitude));
	declination = sin(latitude) * cos(nutation.ecliptic) + cos(latitude) *
		sin(nutation.ecliptic) * sin(longitude);
	declination = asin(declination);

	/* store in position */
	position.ra = ln_range_degrees(ln_rad_to_deg(ra));
	position.dec = ln_rad_to_deg(declination);
}

/*! \fn void ln_get_ecl_from_equ(struct ln_equ_posn *object, double JD, struct ln_lnlat_posn *position)
* \param object Object coordinates in B1950. Use ln_get_equ_prec2 to transform from J2000.
* \param JD Julian day
* \param position Pointer to store new position.
*
* Transform an objects equatorial cordinates into ecliptical coordinates
* for the given julian day.
*/
/* Equ 12.1, 12.2 Pg 88
*/
@nogc void ln_get_ecl_from_equ(const ref ln_equ_posn object, double JD,
	ref ln_lnlat_posn position) nothrow
{
	double ra, declination, latitude, longitude;
	ln_nutation nutation;

	/* object position */
	ra = ln_deg_to_rad(object.ra);
	declination = ln_deg_to_rad(object.dec);
	ln_get_nutation(JD, nutation);
	nutation.ecliptic = ln_deg_to_rad(nutation.ecliptic);

	/* Equ 12.1, 12.2 */
	longitude = atan2((sin(ra) * cos(nutation.ecliptic) + tan(declination) *
		sin(nutation.ecliptic)), cos(ra));
	latitude = sin(declination) * cos(nutation.ecliptic) - cos(declination) *
		sin(nutation.ecliptic) * sin(ra);
	latitude = asin(latitude);

	/* store in position */
	position.lat = ln_rad_to_deg(latitude);
	position.lng = ln_range_degrees(ln_rad_to_deg(longitude));
}

/*! \fn void ln_get_ecl_from_rect(struct ln_rect_posn *rect, struct ln_lnlat_posn *posn)
* \param rect Rectangular coordinates.
* \param posn Pointer to store new position.
*
* Transform an objects rectangular coordinates into ecliptical coordinates.
*/
/* Equ 33.2
*/
@nogc void ln_get_ecl_from_rect(const ref ln_rect_posn rect, ref ln_lnlat_posn posn) nothrow
{
	double t;

	t = sqrt(rect.X * rect.X + rect.Y * rect.Y);
	posn.lng = ln_range_degrees(ln_rad_to_deg(atan2(rect.X, rect.Y)));
	posn.lat = ln_rad_to_deg(atan2(t, rect.Z));
}

/*! \fn void ln_get_equ_from_gal(struct ln_gal_posn *gal, struct ln_equ_posn *equ)
* \param gal Galactic coordinates.
* \param equ B1950 equatorial coordinates. Use ln_get_equ_prec2 to transform to J2000.
*
* Transform an object galactic coordinates into B1950 equatorial coordinate.
*/
/* Pg 94 */
@nogc void ln_get_equ_from_gal(const ref ln_gal_posn gal, ref ln_equ_posn equ) nothrow
{
	double RAD_27_4, SIN_27_4, COS_27_4;
	double l_123, cos_l_123;
	double sin_b, cos_b, rad_gal_b;
	double y;

	RAD_27_4 = ln_deg_to_rad(27.4);
	SIN_27_4 = sin(RAD_27_4);
	COS_27_4 = cos(RAD_27_4);

	l_123 = ln_deg_to_rad(gal.l - 123);
	cos_l_123 = cos(l_123);

	rad_gal_b = ln_deg_to_rad(gal.b);

	sin_b = sin(rad_gal_b);
	cos_b = cos(rad_gal_b);

	y = atan2(sin(l_123), cos_l_123 * SIN_27_4 - (sin_b / cos_b) * COS_27_4);
	equ.ra = ln_range_degrees(ln_rad_to_deg(y) + 12.25);
	equ.dec =
		ln_rad_to_deg(asin(sin_b * SIN_27_4 + cos_b * COS_27_4 * cos_l_123));
}

/*! \fn void ln_get_equ2000_from_gal(struct ln_gal_posn *gal, struct ln_equ_posn *equ)
* \param gal Galactic coordinates.
* \param equ J2000 equatorial coordinates.
*
* Transform an object galactic coordinates into equatorial coordinate.
*/
@nogc void ln_get_equ2000_from_gal(const ref ln_gal_posn gal, ref ln_equ_posn equ) nothrow
{
	ln_get_equ_from_gal(gal, equ);
	ln_get_equ_prec2(equ, B1950, JD2000, equ);
}

/*! \fn ln_get_gal_from_equ(struct ln_equ_posn *equ, struct ln_gal_posn *gal)
* \param equ B1950 equatorial coordinates.
* \param gal Galactic coordinates.
*
* Transform an object B1950 equatorial coordinate into galactic coordinates.
*/
/* Pg 94 */
@nogc void ln_get_gal_from_equ(const ref ln_equ_posn equ, ref ln_gal_posn gal) nothrow
{
	double RAD_27_4, SIN_27_4, COS_27_4;
	double ra_192_25, cos_ra_192_25;
	double rad_equ_dec;
	double cos_dec, sin_dec;
	double x;

	RAD_27_4 = ln_deg_to_rad(27.4);
	SIN_27_4 = sin(RAD_27_4);
	COS_27_4 = cos(RAD_27_4);

	ra_192_25 = ln_deg_to_rad(192.25 - equ.ra);
	cos_ra_192_25 = cos(ra_192_25);

	rad_equ_dec = ln_deg_to_rad(equ.dec);

	sin_dec = sin(rad_equ_dec);
	cos_dec = cos(rad_equ_dec);

	x = atan2(sin(ra_192_25),
		cos_ra_192_25 * SIN_27_4 - (sin_dec / cos_dec) * COS_27_4);
	gal.l = ln_range_degrees(303 - ln_rad_to_deg(x));
	gal.b = ln_rad_to_deg(asin(sin_dec * SIN_27_4 + cos_dec *
		COS_27_4 * cos_ra_192_25));
}

/*! \fn void ln_get_gal_from_equ2000(struct ln_equ_posn *equ, struct ln_gal_posn *gal)
* \param equ J2000 equatorial coordinates.
* \param gal Galactic coordinates.
*
* Transform an object J2000 equatorial coordinate into galactic coordinates.
*/
@nogc void ln_get_gal_from_equ2000(const ref ln_equ_posn equ, ref ln_gal_posn gal) nothrow
{
	ln_equ_posn equ_1950;
	ln_get_equ_prec2(equ, JD2000, B1950, equ_1950);
	ln_get_gal_from_equ(equ_1950, gal);
}

/*! \example transforms.c
 *
 * Examples of how to use transformation functions.
 */

}