// orbit.h
//
// Copyright (C) 2001 Chris Laurel <claurel@shatters.net>
//
// 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 3
// of the License, or (at your option) any later version.

#ifndef _ORBIT_H_
#define _ORBIT_H_

#include "tools/vecmath.hpp"
#include "iterative_orbits.hpp"
#include <string>
#include <memory>

// The callback type for the external position computation function
typedef void (PositionFunctionType)(double jd,double xyz[3]);
typedef void (OsculatingFunctionType)(double jd0,double jd,double xyz[3]);

class Body;

class Orbit {
public:
	Orbit(void) {}
	virtual ~Orbit(void) {}

	// Compute position for a specified Julian date and return coordinates
	// given in "dynamical equinox and ecliptic J2000"
	// which is the reference frame for VSOP87
	virtual void positionAtTimevInVSOP87Coordinates(double, double, double*) const = 0;
	void positionAtTimevInVSOP87Coordinates(double JD, double* v) const {
		positionAtTimevInVSOP87Coordinates(JD, JD, v);
	}

	// As fastPositionAtTimevInVSOP87Coordinates is const, give a chance to cache some datas
	virtual std::pair<double, double> prepairFastPositionAtTimevInVSOP87Coordinates(double JD0, double deltaJD) {
		return {-deltaJD, deltaJD};
	}

	// If possible, do faster (and less accurate) calculation for orbits
	virtual void fastPositionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const {
		positionAtTimevInVSOP87Coordinates(JD, JD, v);
	}
	// //! updating comet tails is a bit expensive. try not to overdo it.
	// virtual bool getUpdateTails() const;
	// virtual void setUpdateTails(const bool update);

	virtual OsculatingFunctionType * getOsculatingFunction() const {
		return nullptr;
	};

	virtual double getBoundingRadius() const {
		return 0;
	}

	// Is this orbit stable (exactly the same path over time)
	virtual bool isStable(double) const {
		return true;
	}

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double) const {
		return true;
	}

	virtual std::string saveOrbit() const = 0;

private:
	Orbit(const Orbit&);
	const Orbit &operator=(const Orbit&);
};


class EllipticalOrbit : public Orbit {
public:
	EllipticalOrbit(double pericenterDistance,
	                double eccentricity,
	                double inclination,
	                double ascendingNode,
	                double argOfPeriapsis,
	                double meanAnomalyAtEpoch,
	                double period,
	                double epoch, // = 2451545.0,
	                double _parent_rot_obliquity, // = 0.0,
	                double _parent_rot_ascendingnode, // = 0.0
	                double _parent_rot_J2000_longitude,
	                bool useParentPrecession=true);

	// Compute position for a specified Julian date and return coordinates
	// given in "dynamical equinox and ecliptic J2000"
	// which is the reference frame for VSOP87
	// In order to rotate to VSOP87
	// parent_rot_obliquity and parent_rot_ascendingnode must be supplied.
	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	virtual std::pair<double, double> prepairFastPositionAtTimevInVSOP87Coordinates(double JD0, double deltaJD) override;
	virtual void fastPositionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	// Original one
	Vec3d positionAtTime(double) const;
	double getPeriod() const;
	double getBoundingRadius() const;

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double) const {
		return m_UseParentPrecession;
	}

	virtual std::string saveOrbit() const;

private:
	double eccentricAnomaly(double m, double &lastE) const;
	Vec3d positionAtE(double) const;

	//! Last value returned by eccentricAnomaly, used for iterative precision
	mutable double iterativeLastE = 0;
	//! Last value returned by eccentricAnomaly, used for fast precision with batch position
	mutable double batchLastE = 0;

	double pericenterDistance;
	double eccentricity;
	double inclination;
	double ascendingNode;
	double argOfPeriapsis;
	double meanAnomalyAtEpoch;
	double period;
	double epoch;

	// Rotation to VSOP87 coordinate data
	// \todo replace with IAU formulas for planets
	double parent_rot_obliquity;
	double parent_rot_ascendingnode;
	double parent_rot_J2000_longitude;
	double rotate_to_vsop87[9];

	bool m_UseParentPrecession;
};

class CometOrbit : public Orbit {
public:
	// Prepair orbit computation around JD0 and rectify the deltaJD
	virtual std::pair<double, double> prepairFastPositionAtTimevInVSOP87Coordinates(double JD0, double deltaJD) override;
	// Override position for orbit computation so they are more equally spaced
	virtual void fastPositionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;
	void deltaJDToOrbitJD(double JD, double &deltaJD) const;
	void orbitJDToJD(double &JD) const;
	// //! updating comet tails is a bit expensive. try not to overdo it.
	// virtual bool getUpdateTails() const override{ return updateTails; }
	// virtual void setUpdateTails(const bool update) override{ updateTails=update; }
	double getPeriod() const;
	double getBoundingRadius() const;
	virtual std::string saveOrbit() const;

protected:
	CometOrbit(double pericenter_distance,
	           double eccentricity,
	           double inclination,
	           double ascendingNode,
	           double arg_of_perhelion,
	           double time_at_perihelion,
	           double mean_motion,
	           double parent_rot_obliquity,
	           double parent_rot_ascendingnode,
	           double parent_rot_J2000_longitude);

	const double q;
	const double e;
	const double t0;
	const double n;
	double orbitJDCorrection;
	Vec3d d1;
	Vec3d d2;
	double rotate_to_vsop87[9];
};

class HypCometOrbit : public CometOrbit {
public:
	HypCometOrbit(double pericenter_distance,
	           double eccentricity,
	           double inclination,
	           double ascendingNode,
	           double arg_of_perhelion,
	           double time_at_perihelion,
	           double mean_motion,
	           double parent_rot_obliquity,
	           double parent_rot_ascendingnode,
	           double parent_rot_J2000_longitude);

    void warp(double JD);
   // Compute the orbit for a specified Julian date and return an "application compliant" function
   	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;
   	Vec3d positionAtTime(double) const;
private:
	mutable IterativeHyp orbit; // Mitigation - mutable car positionAtTimevInVSOP87Coordinates ne devrais pas etre const...
};

class EllCometOrbit : public CometOrbit {
public:
	EllCometOrbit(double pericenter_distance,
	           double eccentricity,
	           double inclination,
	           double ascendingNode,
	           double arg_of_perhelion,
	           double time_at_perihelion,
	           double mean_motion,
	           double parent_rot_obliquity,
	           double parent_rot_ascendingnode,
	           double parent_rot_J2000_longitude);


   // Compute the orbit for a specified Julian date and return an "application compliant" function
   	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;
   	Vec3d positionAtTime(double) const;
private:
	mutable IterativeEll orbit; // Mitigation - mutable car positionAtTimevInVSOP87Coordinates ne devrais pas etre const...
};

class ParCometOrbit : public CometOrbit {
public:
	ParCometOrbit(double pericenter_distance,
	           double eccentricity,
	           double inclination,
	           double ascendingNode,
	           double arg_of_perhelion,
	           double time_at_perihelion,
	           double mean_motion,
	           double parent_rot_obliquity,
	           double parent_rot_ascendingnode,
	           double parent_rot_J2000_longitude);


   // Compute the orbit for a specified Julian date and return an "application compliant" function
   	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;
	// If possible, do faster (and less accurate) calculation for orbits
	virtual void fastPositionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;
   	Vec3d positionAtTime(double) const;
};

//! A Special Orbit uses special ephemeris algorithms
class SpecialOrbit : public Orbit {
public:
	SpecialOrbit(std::string ephemerisName);
	virtual ~SpecialOrbit(void) {};

	// Compute position for a specified Julian date and return coordinates
	// given in "dynamical equinox and ecliptic J2000"
	// which is the reference frame for VSOP87
	// In order to rotate to VSOP87
	// parent_rot_obliquity and parent_rot_ascendingnode must be supplied.
	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	virtual OsculatingFunctionType * getOsculatingFunction() const {
		return osculatingFunction;
	}

	// An ephemeris orbit is not exactly the same through time
	virtual bool isStable(double) const {
		return stable;
	}

	virtual bool isValid() {
		return (positionFunction != nullptr);
	}

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double) const {
		return m_UseParentPrecession;
	}
	virtual std::string saveOrbit() const;


private:
	PositionFunctionType *positionFunction;
	OsculatingFunctionType *osculatingFunction;
	bool stable;  // does not osculate noticeably for performance caching orbit visualization
	bool m_UseParentPrecession;
};


/*! A mixed orbit is a composite orbit, typically used when you have a
 *  custom orbit calculation that is only valid over limited span of time.
 *  When a mixed orbit is constructed, it computes elliptical orbits
 *  to approximate the behavior of the primary orbit before and after the
 *  span over which it is valid.
 */
class MixedOrbit : public Orbit {
public:
	MixedOrbit(std::unique_ptr<Orbit> orbit, double period, double t0, double t1, double mass,
	           double _parent_rot_obliquity, // = 0.0,
	           double _parent_rot_ascendingnode, // = 0.0
	           double _parent_rot_J2000_longitude,
	           bool useParentPrecession=true);

	virtual ~MixedOrbit();

	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	virtual bool isStable(double jd) const;

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double jd) const;

	virtual std::string saveOrbit() const;

private:
	std::unique_ptr<Orbit> primary;
	Orbit* afterApprox;
	Orbit* beforeApprox;
	double begin;
	double end;
	double boundingRadius;
};



/*! A binary orbit is a two body gravitational system where each orbits
 *  a common barycenter, which itself orbits a parent.  One body is
 *  designated the primary body, the other the secondary.
 *  Primary body position is calculated from secondary body position,
 *  the barycenter position, and the relative masses.
 *  Example: the Earth/Moon system.
 *  Because primary body is loaded before the secondary at startup
 *  the secondary has to be added after construction.
 *  Note that if we define r =  mass(secondary)/mass(primary)
 *  ratio = r/(1+r)
 */
class BinaryOrbit : public Orbit {
public:
	BinaryOrbit(std::unique_ptr<Orbit> barycenter, double ratio);

	virtual ~BinaryOrbit();

	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	// If possible, do faster (and less accurate) calculation for orbits
	virtual void fastPositionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	virtual bool isStable(double jd) const;

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double jd) const {
		return barycenter->useParentPrecession(jd);
	}

	virtual void setSecondaryOrbit(Orbit *second) {
		secondary = second;
	}

	virtual std::string saveOrbit() const;

private:
	std::unique_ptr<Orbit> barycenter;
	Orbit* secondary;
	double ratio;

};

/*! An orbit still represents a fixed orbit with respect to its reference frame
 * Example: geostationary satellite frozen
 * It is mainly used for observation purposes
 *
 */
class stillOrbit : public Orbit {
public:
	//! creation of an orbit still at the Cartesian coordinate point (x,y,z)
	stillOrbit(double _x, double _y, double _z);

	virtual ~stillOrbit();

	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double jd) const {
		return false;
	}

	std::string saveOrbit() const;

private:
	double x,y,z;
};

class LocationOrbit : public Orbit {
public:
	LocationOrbit(double _lon, double _lat, double _alt, double parentRadius, double parentPeriod, double parentOffset);

	virtual ~LocationOrbit();

	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	std::string saveOrbit() const;

private:
	double lon, lat, alt, JDToRotation;
};

class linearOrbit : public Orbit {
public:
	linearOrbit(double _t_start, double _t_end, double *_posInitial, double *_posFinal );
	virtual ~linearOrbit();

	virtual void positionAtTimevInVSOP87Coordinates(double JD0, double JD, double *v) const;

	// Do the body coordinates precess with the parent?
	virtual bool useParentPrecession(double jd) const {
		return false;
	}

	std::string saveOrbit() const;


private:
	double t_start;
	double t_end;
	float t_duration;
	Vec3d posInitial;
	Vec3d posFinal;
};

class BarycenterOrbit : public Orbit {
public:
	BarycenterOrbit() = delete;
	BarycenterOrbit(std::shared_ptr<Body> bodyA, std::shared_ptr<Body> bodyB, double a, double b);
	virtual ~BarycenterOrbit() { }

	void positionAtTimevInVSOP87Coordinates(double, double, double*) const;

	bool useParentPrecession(double) const {
		return false;
	}

	std::string saveOrbit() const;

private:
	std::shared_ptr<Body> bodyA, bodyB;
	double a,b;
};

#endif // _ORBIT_H_
