/* 
   ORSA - Orbit Reconstruction, Simulation and Analysis
   Copyright (C) 2002-2004 Pasquale Tricarico
   
   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.
   
   As a special exception, Pasquale Tricarico gives permission to
   link this program with Qt commercial edition, and distribute the
   resulting executable, without including the source code for the Qt
   commercial edition in the source distribution.
   
   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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#ifndef _ORSA_ORBIT_GSL_H_
#define _ORSA_ORBIT_GSL_H_

#include <string>

#include "orsa_orbit.h"
#include "orsa_universe.h"

namespace orsa {
  
  enum CovarianceMatrixElements { Osculating, // a,e,i,node,peri,M
				  Equinoctal  // a,k,h,q,p,lambda
  };
  
  class OrbitWithCovarianceMatrixGSL : public OrbitWithEpoch {
  public:
    OrbitWithCovarianceMatrixGSL();
    OrbitWithCovarianceMatrixGSL(Orbit&);
    OrbitWithCovarianceMatrixGSL(Orbit&,double**,CovarianceMatrixElements);
  public:
    void SetCovarianceMatrix(double**,CovarianceMatrixElements);
    void GetCovarianceMatrix(double**,CovarianceMatrixElements&) const;
  private:
    double covm[6][6];
    CovarianceMatrixElements cov_base;
  public:
    // void generate(std::vector<OrbitWithEpoch>&,const int, const int=12345) const;
    void GenerateUsingCholeskyDecomposition(std::vector<OrbitWithEpoch>&,const int, const int=12345) const;
    void GenerateUsingPrincipalAxisTransformation(std::vector<OrbitWithEpoch>&,const int, const int=12345) const;
  public:
    bool have_covariance_matrix() const { return bool_have_covariance_matrix; }
  private:
    bool bool_have_covariance_matrix;
  };
  
  // class PreliminaryOrbit : public OrbitWithEpoch {
  class PreliminaryOrbit : public OrbitWithCovarianceMatrixGSL {
  public:
    double GetRMS() const { return rms; };
    
  public:	
    void ComputeRMS(const std::vector<Observation>&);
    
  public:
    // used to sort 
    inline bool operator < (const PreliminaryOrbit &orbit) const {
      return (rms < orbit.rms);
    }
    
  private:
    double rms;
  };
  
  //! Minimal Orbit Intersection Distance between two orbits
  double MOID(const Orbit&, const Orbit&, Vector&, Vector&);
  
  //! Minimal Orbit Intersection Distance between two orbits with different reference bodies  
  double MOID2RB(const Vector&, const Vector&, const Orbit&, const Orbit&, Vector&, Vector&);
  
  //! General interface to the Orbit computation from a set of observations.
  OrbitWithCovarianceMatrixGSL Compute(const std::vector<Observation> &);
  
  void Compute_Gauss(const std::vector<Observation>&, std::vector<PreliminaryOrbit>&);
  void Compute_TestMethod(const std::vector<Observation>&, std::vector<PreliminaryOrbit>&);
  
  //
  
  class Asteroid {
  public:
    int n;    
    OrbitWithCovarianceMatrixGSL orb;
    std::string name;
    double mag;
  };
  
  class AsteroidDatabase : public std::vector<Asteroid> {
    
  };
  
  class ObservationCandidate { 
  public:
    Asteroid        a;
    Observation     o;
    Angle           delta;
  };
  
  // Close Approaches
  
  class CloseApproach {
  public:
    std::string name;
    UniverseTypeAwareTime epoch;
    double distance;
    double relative_velocity;
    // target plane xy and z components
    double tp_xy, tp_z;
  };
  
  // void ComputeCloseApproaches(const Frame &, const unsigned int, const UniverseTypeAwareTime &, const UniverseTypeAwareTime &, std::vector<CloseApproach> &, const double, const double);
  // OLD!
  // void ComputeCloseApproaches(const Frame &, const unsigned int, const UniverseTypeAwareTime &, const UniverseTypeAwareTime &, std::vector<CloseApproach> &, const double, const UniverseTypeAwareTimeStep &);
  
  // void SearchCloseApproaches(const Evolution *, const unsigned int, std::vector<CloseApproach> &, const double, const bool only_planets=true);
  
  void SearchCloseApproaches(const Evolution *, const unsigned int, const unsigned int, std::vector<CloseApproach> &, const double, const double = FromUnits(1,SECOND));
  
  // test
  void OrbitDifferentialCorrectionsLeastSquares(OrbitWithCovarianceMatrixGSL&, const std::vector<Observation>&);
  
}

#endif // _ORSA_ORBIT_GSL_H_
