/* 
   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_FFT_H_
#define _ORSA_FFT_H_

#include "orsa_analysis.h"
#include "orsa_orbit.h"

#include <vector>

#include <fftw.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_min.h>

namespace orsa {
  
  // gsl d=1 minimization stuff
  typedef struct gsl_d1_minimization_parameters {
    fftw_complex * pointer_points_sequence;
    int size;
  } gsl_d1_minimization_parameters;
  
  class FFTDataStructure {
  public:
    double time,amplitude,phase;
    virtual ~FFTDataStructure() {};
  };
  
  class FFTDataStream : public std::vector<FFTDataStructure> {
  public:
    double timestep;
  };
  
  class FFTPowerSpectrumBaseElement {
  public:
    double frequency;
    double power;
  };
  
  class FFTPowerSpectrum : public std::vector<FFTPowerSpectrumBaseElement> {};
  
  class Peak {
  public:
    inline void Set(double f, double a, double p) {
      frequency = f;
      amplitude = a;
      phase     = p;
    } 
    
  public:
    /* 
       // used to sort peaks by increasing frequency 
       inline bool operator < (const Peak &p) const {
       return (frequency < p.frequency);
       }
    */
    // used to sort peaks by decreasing amplitude
    inline bool operator < (const Peak &p) const {
      return (p.amplitude < amplitude);
    }
    
  public:
    double frequency,amplitude,phase;
    
  public:
    fftw_complex phiR, phiL;
    
  public:
    virtual ~Peak() {};
  };
  
  //! signal types
  enum FFTSearch {HK=0,D,PQ,NODE,ANOMALY,ANOMALY_PHASE,A_M,TESTING};
  
  enum FFTSearchAmplitude {A,E,I,SIN_I,TAN_I_2,ONE};
  enum FFTSearchPhase     {OMEGA_NODE,OMEGA_PERICENTER,OMEGA_TILDE,MM,LAMBDA,ZERO};
  
  enum FFTAlgorithm { 
    algo_FFT,   // Search only power spectrum peaks
    algo_FFTB,  // Search only power spectrum peaks, on positive and negative frequencies
    algo_MFT,   // Original method by Laskar
    algo_FMFT1, // MFT with linear freq. corrections (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137)
    algo_FMFT2  // FMFT1 with addition non-linear corrections (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137)
  };
  
  //! FFT analysis 
  class FFT : public Analysis {
    
  public: 
    FFT(OrbitStream&, FFTPowerSpectrum&, std::vector<Peak>&, FFTDataStream&);
    
    void Search(FFTSearch,FFTAlgorithm=algo_FFT);
    void Search(FFTSearchAmplitude,FFTSearchPhase,FFTAlgorithm=algo_FFT);
    
  private:
    void Search_FFT();   // FourierTransformPeaks
    void Search_FFTB();
    void Search_MFT();   // ModifiedFourierTransform
    void Search_FMFT1(); // FrequencyModifiedFourierTransform_1
    void Search_FMFT2(); // FrequencyModifiedFourierTransform_2
    //
    void Search_FMFT_main();
    
  public:   
    //! frequency range parameters
    // double freq_start, freq_stop;
    
  public:
    std::vector<int> candidate_bin;
    
  public:
    double relative_amplitude;
      
  private:
    double default_peak_reset_frequency;
    double default_peak_reset_amplitude;
    double default_peak_reset_phase;
    
  public:
    bool  HiResSpectrum;
    
  private:
    void FillDataStream(FFTSearch);
    void FillDataStream(FFTSearchAmplitude,FFTSearchPhase);
    void ComputeCommonPowerSpectrum();
    void ComputeCommonReconstructedSignal();
    
  public:
    gsl_d1_minimization_parameters par;
    gsl_function F;
    gsl_min_fminimizer *s;
    const gsl_min_fminimizer_type *T;
    
  protected:
    std::vector<double>     psd;
    FFTPowerSpectrum  *fps;

    std::vector<Peak>   *peaks;
    OrbitStream    *os;
    
    FFTDataStream  *reconstructed_data_stream;    
    
  public:
    FFTDataStream  data_stream; // should be used readonly from outside this class...
    
  public:
    unsigned int nfreq;
  };
  
} // namespace orsa

#endif // _ORSA_FFT_H_

