/*
 *   Copyright (C) 1992-2003 by John A. Magliacane KD2BD
 *   Copyright (C) 2003 by Jonathan Naylor G4KLX
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program 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 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "Sun.h"
#include "Inline.h"

#include <cmath>
using namespace std;

const double R  = 6378.16;
const double C1 = 0.9174640600000014;
const double S1 = 0.397818675;

CSun::CSun() :
m_az(0.0),
m_el(0.0),
m_ra(0.0),
m_dec(0.0),
m_range(0.0),
m_obs()
{
}

CSun::~CSun()
{
}

void CSun::findSun(double dayNum)
{
	/* This function finds the position of the Sun */

	/* Zero vector for initializations */
	vector_t zeroVector = {0, 0, 0, 0};

	/* Solar ECI position vector  */
	vector_t solarVector = zeroVector;

	/* Solar observed azi and ele vector  */
	vector_t solarSet = zeroVector;

	/* Solar right ascension and declination vector */
	vector_t solarRad = zeroVector;

	/* Solar lat, long, alt vector */
	geodetic_t solarLatLonAlt;

	double julUTC = dayNum + 2444238.5;

	calculateSolarPosition(julUTC, solarVector);
	calculateObs(julUTC, solarVector, zeroVector, m_obs, solarSet);

	m_az    = solarSet.x;
	m_el    = solarSet.y;
	m_range = 1.0 + ((solarSet.z - AU) / AU);

	calculateLatLonAlt(julUTC, solarVector, solarLatLonAlt);
	calculateRADec(julUTC, solarVector, zeroVector, m_obs, solarRad);

	m_ra  = solarRad.x;
	m_dec = solarRad.y;
}

void CSun::calculateSolarPosition(double time, vector_t& solarVector) const
{
	double mjd  = time - 2415020.0;
	double year = 1900 + mjd / 365.25;
	double    T = (mjd + deltaET(year) / secday) / 36525.0;

	double M = RAD(modulus(358.47583 + modulus(35999.04975 * T, 360.0) - (0.000150 + 0.0000033 * T) * T * T, 360.0));

	double L = RAD(modulus(279.69668 + modulus(36000.76892 * T, 360.0) + 0.0003025 * T * T, 360.0));
	double e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T;

	double C   = RAD((1.919460 - (0.004789 + 0.000014 * T) * T) * ::sin(M) + (0.020094 - 0.000100 * T) * ::sin(2.0 * M) + 0.000293 * ::sin(3.0 * M));
	double O   = RAD(modulus(259.18 - 1934.142 * T, 360.0));
	double Lsa = modulus(L + C- RAD(0.00569 - 0.00479 * ::sin(O)), 2.0 * M_PI);
	double nu  = modulus(M + C, 2.0 * M_PI);
	double R   = 1.0000002 * (1.0 - e * e) / (1.0 + e * ::cos(nu));

	double eps = RAD(23.452294 - (0.0130125 + (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * ::cos(O));
	R = AU * R;

	solarVector.x = R * ::cos(Lsa);
	solarVector.y = R * ::sin(Lsa) * ::cos(eps);
	solarVector.z = R * ::sin(Lsa) * ::sin(eps);
	solarVector.w = R;
}

void CSun::calculateObs(double time, vector_t& pos, vector_t& vel, geodetic_t& geodetic, vector_t& obsSet) const
{
	/* The procedures calculateObs and calculateRADec calculate           */
	/* the *topocentric* coordinates of the object with ECI position,     */
	/* {pos}, and velocity, {vel}, from location {geodetic} at {time}.    */
	/* The {obs_set} returned for calculateObs consists of azimuth,       */
	/* elevation, range, and range rate (in that order) with units of     */
	/* radians, radians, kilometers, and kilometers/second, respectively. */
	/* The WGS '72 geoid is used and the effect of atmospheric refraction */
	/* (under standard temperature and pressure) is incorporated into the */
	/* elevation calculation; the effect of atmospheric refraction on     */
	/* range and range rate has not yet been quantified.                  */

	/* The {obs_set} for calculateRADec consists of right ascension and   */
	/* declination (in that order) in radians.  Again, calculations are   */
	/* based on *topocentric* position using the WGS '72 geoid and        */
	/* incorporating atmospheric refraction.                              */

	vector_t obsPos, obsVel;
	calculateUserPosVel(time, geodetic, obsPos, obsVel);

	vector_t range;
	range.x = pos.x - obsPos.x;
	range.y = pos.y - obsPos.y;
	range.z = pos.z - obsPos.z;

	vector_t rgVel;
	rgVel.x = vel.x - obsVel.x;
	rgVel.y = vel.y - obsVel.y;
	rgVel.z = vel.z - obsVel.z;

	magnitude(range);

	double sinLat = ::sin(geodetic.lat);
	double cosLat = ::cos(geodetic.lat);
	double sinTheta = ::sin(geodetic.theta);
	double cosTheta = ::cos(geodetic.theta);

	double topS = sinLat * cosTheta * range.x + sinLat * sinTheta * range.y - cosLat * range.z;
	double topE = -sinTheta * range.x + cosTheta * range.y;
	double topZ = cosLat * cosTheta * range.x + cosLat * sinTheta * range.y + sinLat * range.z;
	double azim = ::atan(-topE / topS);		/* Azimuth */

	if (topS > 0.0)
                azim += M_PI;

	if (azim < 0.0)
                azim += 2.0 * M_PI;

	double el = ::asin(topZ / range.w);

	obsSet.x = azim;				/* Azimuth (radians)   */
	obsSet.y = el;					/* Elevation (radians) */
	obsSet.z = range.w;				/* Range (kilometers)  */

	/* Range Rate (kilometers/second) */

	obsSet.w = dot(range, rgVel) / range.w;

	/* Corrections for atmospheric refraction */
	/* Reference:  Astronomical Algorithms by Jean Meeus, pp. 101-104    */
	/* Correction is meaningless when apparent elevation is below horizon */

	if (obsSet.y < 0.0)
                obsSet.y = el;		/* Reset to true elevation */
}

void CSun::calculateLatLonAlt(double time, vector_t& pos,  geodetic_t& geodetic) const
{
	/* Procedure calculateLatLonAlt will calculate the geodetic   */
	/* position of an object given its ECI position pos and time. */
	/* It is intended to be used to determine the ground track of */
	/* a satellite.  The calculations  assume the earth to be an  */
	/* oblate spheroid as defined in WGS '72.                     */

	/* Reference:  The 1992 Astronomical Almanac, page K12. */

	geodetic.theta = ::atan2(pos.y, pos.x);				/* radians */
	geodetic.lon   = modulus(geodetic.theta - thetaGJD(time), 2.0 * M_PI);/* radians */

	double  r = ::sqrt(pos.x * pos.x + pos.y * pos.y);
	double e2 = f * (2.0 - f);

	geodetic.lat = ::atan2(pos.z, r);					/* radians */

	double phi, c;
	do {
		phi = geodetic.lat;

		c = 1.0 / ::sqrt(1.0 - e2 * ::pow(::sin(phi), 2.0));

		geodetic.lat = ::atan2(pos.z + xkmper * c * e2 * ::sin(phi), r);
	} while (::fabs(geodetic.lat - phi) >= 1E-10);

	geodetic.alt = r / ::cos(geodetic.lat) - xkmper * c;		/* kilometers */

	if (geodetic.lat > (M_PI / 2.0))
		geodetic.lat -= 2.0 * M_PI;
}

void CSun::calculateRADec(double time, vector_t& pos, vector_t& vel, geodetic_t& geodetic, vector_t& obsSet) const
{
	/* Reference:  Methods of Orbit Determination by  */
	/*             Pedro Ramon Escobal, pp. 401-402   */

	calculateObs(time, pos, vel, geodetic, obsSet);

	double az = obsSet.x;
	double el = obsSet.y;
	double phi = geodetic.lat;
	double theta = modulus(thetaGJD(time) + geodetic.lon, 2.0 * M_PI);
	double sinTheta = ::sin(theta);
	double cosTheta = ::cos(theta);
	double sinPhi = ::sin(phi);
	double cosPhi = ::cos(phi);
	double Lxh = -::cos(az) * ::cos(el);
	double Lyh = ::sin(az) * ::cos(el);
	double Lzh = ::sin(el);
	double Sx = sinPhi * cosTheta;
	double Ex = -sinTheta;
	double Zx = cosTheta * cosPhi;
	double Sy = sinPhi * sinTheta;
	double Ey = cosTheta;
	double Zy = sinTheta * cosPhi;
	double Sz = -cosPhi;
	double Ez = 0.0;
	double Zz = sinPhi;
	double Lx = Sx * Lxh + Ex * Lyh + Zx *Lzh;
	double Ly = Sy * Lxh + Ey * Lyh + Zy * Lzh;
	double Lz = Sz * Lxh + Ez * Lyh + Zz * Lzh;

        obsSet.y = ::asin(Lz);			/* Declination (radians) */
	double cosDelta = ::sqrt(1.0 - Lz * Lz);
	double sinAlpha = Ly / cosDelta;
	double cosAlpha = Lx / cosDelta;
	obsSet.x = ::atan2(sinAlpha, cosAlpha);	/* Right Ascension (radians) */
        obsSet.x = modulus(obsSet.x, 2.0 * M_PI);
}

void CSun::calculateUserPosVel(double time, geodetic_t& geodetic, vector_t& obsPos, vector_t& obsVel) const
{
	/* Calculate_User_PosVel() passes the user's geodetic position
	   and the time of interest and returns the ECI position and
	   velocity of the observer.  The velocity calculation assumes
	   the geodetic position is stationary relative to the earth's
	   surface. */

	/* Reference:  The 1992 Astronomical Almanac, page K11. */

	geodetic.theta = modulus(thetaGJD(time) + geodetic.lon, 2.0 * M_PI);/* LMST */
	double c = 1.0 / ::sqrt(1.0 + f *( f - 2.0) * ::pow(::sin(geodetic.lat), 2.0));
	double sq = ::pow(1.0 - f, 2.0) * c;
	double achcp = (xkmper * c + geodetic.alt) * ::cos(geodetic.lat);

	obsPos.x = achcp * ::cos(geodetic.theta);	/* kilometers */
	obsPos.y = achcp * ::sin(geodetic.theta);
	obsPos.z = (xkmper * sq + geodetic.alt) * ::sin(geodetic.lat);

        obsVel.x = -mfactor * obsPos.y;		/* kilometers/second */
        obsVel.y = mfactor * obsPos.x;
	obsVel.z = 0.0;

	magnitude(obsPos);
	magnitude(obsVel);
}

void CSun::magnitude(vector_t& v) const
{
	/* Calculates scalar magnitude of a vector_t argument */
	v.w = ::sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

double CSun::deltaET(double year) const
{
	/* The function deltaET has been added to allow calculations on    */
	/* the position of the sun.  It provides the difference between UT */
	/* (approximately the same as UTC) and ET (now referred to as TDT).*/
	/* This function is based on a least squares fit of data from 1950 */
	/* to 1991 and will need to be updated periodically. */

	/* Values determined using data from 1950-1991 in the 1990
	   Astronomical Almanac.  See DELTA_ET.WQ1 for details. */

	return 26.465 + 0.747622 * (year - 1950.0) + 1.886913 * ::sin(2.0 * M_PI * (year - 1975.0) / 33.0);
}

double CSun::thetaGJD(double jd) const
{
	/* Reference:  The 1992 Astronomical Almanac, page B6. */

	double UT = frac(jd + 0.5);
	jd -= UT;
	double TU = (jd - 2451545.0) / 36525.0;
	double GMST = 24110.54841 + TU * (8640184.812866 + TU * (0.093104 - TU * 6.2E-6));
	GMST = modulus(GMST + secday * omegaE * UT, secday);

	return 2.0 * M_PI * GMST / secday;
}

double CSun::modulus(double arg1, double arg2) const
{
	// Returns arg1 mod arg2

	double ret = arg1;
	int i = int(ret / arg2);
	ret -= double(i) * arg2;

	if (ret < 0.0)
		ret += arg2;

	return ret;
}

double CSun::dot(vector_t& v1, vector_t& v2) const
{
	/* Returns the dot product of two vectors */
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

double CSun::frac(double arg) const
{
	/* Returns fractional part of double argument */
	return arg - ::floor(arg);
}

double CSun::getAzimuth() const
{
	return DEG(m_az);
}

double CSun::getElevation() const
{
	return DEG(m_el);
}

double CSun::getRA() const
{
	return DEG(m_ra);
}

double CSun::getDec() const
{
	return DEG(m_dec);
}

double CSun::getRange() const
{
	return m_range;
}

void CSun::setLocation(double latitude, double longitude, double height)
{
	m_obs.lat = RAD(latitude);
	m_obs.lon = RAD(longitude);
	m_obs.alt = height / 1000.0;
	m_obs.theta = 0.0;
}

double CSun::getDayNum(const wxDateTime& dateTime) const
{
	dateTime.ToGMT(true);

	int  secs = dateTime.GetTicks();
	int msecs = dateTime.GetMillisecond();

	return ((double(secs) + 0.001 * double(msecs)) / 86400.0) - 3651.0;
}