#ifndef _utl_FFTDataContainer_h_
#define _utl_FFTDataContainer_h_ 

/**
  \file
  Template class for a data container that offers and takes both
  time series and corresponding frequency data, keeping track of
  which one is the most recent and doing FFTs automatically if
  necessary.
  
  Caution: Do not cache references you have asked from an
  FFTDataContainer. Always refresh a new one to make sure it is
  up to date.
   
  \author Tim Huege
  \ingroup fft
*/

#include <complex>
#include <utl/Trace.h>
#include <utl/ShadowPtr.h>
#include <utl/AugerException.h>

#include <utl/fft++.h>

namespace utl {

  /**
    \class FFTDataContainer FFTDataContainer.h "utl/FFTDataContainer.h"
    \brief
    Template class for a data container that offers and takes both
    time series and corresponding frequency data, keeping track of
    which one is the most recent and doing FFTs automatically if
    necessary. The data are stored in ShadowPtrs to be able to clear
    them without relying on specific clear() methods of the container
    classes.

    \param C Container class that is used to store the data, usually Trace
    \param T Data type that is used to store time-series data, usually double
    \param F Data type that is used to store frequency-spectrum data, usually complex<double>

  \author Tim Huege
  */

  template<template <typename X> class C, typename T, typename F>
  class FFTDataContainer {

  public:

    /// constructor setting up default values
    FFTDataContainer()
    : fUpToDateDomain(Both), fNyquistZone(1u)
    {}

    /// get the Nyquist zone
    unsigned int GetNyquistZone() const
    {
      return fNyquistZone;
    }

    /// set the Nyquist zone
    void SetNyquistZone(unsigned int zone)
    {
      if (zone < 1)
      {
        throw(OutOfBoundException("Requested illegal Nyquist zone of zero."));
      }
      fNyquistZone=zone;
    }

    /// get the frequency corresponding to a bin of the frequency spectrum
    double GetFrequencyOfBin(typename C<F>::SizeType bin) const
    {
      // maybe this has to be done in a more efficient way ...
      if (fUpToDateDomain == Time)
      {
        // we have to update the cached frequency spectrum data, so we need to cast away the constness of the this-object
        // this is ok since the change is "invisible" to the user, who always gets up-to-date data
        FFTDataContainer<C,T,F>* mutablethis = const_cast<FFTDataContainer<C,T,F>*>(this);
        mutablethis->UpdateFrequencySpectrum(); // do the fft
      }
      double binning = fFrequencySpectrum->GetBinning();
      typename C<F>::SizeType highestbin = fFrequencySpectrum->GetSize()-1;
      if (bin>highestbin)
      {
        throw(OutOfBoundException("Requested frequency for a bin that is not in the FrequencySpectrum."));
      }
      double bandwidth = static_cast<double>(highestbin)*binning;
      unsigned int evenedzone = (fNyquistZone/2)*2;
      double startingfrequency=evenedzone*bandwidth;
      int sign=1;
      if (fNyquistZone==evenedzone)
      {
        sign=-1; // we are in an even Nyquist zone, i.e., the frequencies count downwards
      }
      return startingfrequency+sign*static_cast<double>(bin)*binning;
    }

    /// get the bin of a corresponding frequency of the frequency spectrum
    int GetBinOfFrequency(double frequency) const
    {
      // maybe this has to be done in a more efficient way ...
      if (fUpToDateDomain == Time)
      {
        // we have to update the cached frequency spectrum data, so we need to cast away the constness of the this-object
        // this is ok since the change is "invisible" to the user, who always gets up-to-date data
        FFTDataContainer<C,T,F>* mutablethis = const_cast<FFTDataContainer<C,T,F>*>(this);
        mutablethis->UpdateFrequencySpectrum(); // do the fft
      }
      double binning = fFrequencySpectrum->GetBinning();
      typename C<F>::SizeType highestbin = fFrequencySpectrum->GetSize()-1;
      double bandwidth = static_cast<double>(highestbin)*binning;
      if( ((fNyquistZone-1)*bandwidth > frequency) or (fNyquistZone*bandwidth < frequency) ) {
        throw(OutOfBoundException("Requested bin for a frequency that is not in the FrequencySpectrum."));
      }

      unsigned int evenedzone = (fNyquistZone/2)*2;
      double startingfrequency=evenedzone*bandwidth;
      int sign=1;
      if (fNyquistZone==evenedzone)
      {
        sign=-1; // we are in an even Nyquist zone, i.e., the frequencies count downwards
      }
      double timeBinning = fTimeSeries->GetBinning();
      int timeTraceLength = fTimeSeries->GetSize();

      double bin = (frequency-startingfrequency) * sign *  timeBinning * timeTraceLength;
      return static_cast<int>(std::floor(0.5 + bin));
    }

    /// read out the time series (read access)
    const C<T>& GetConstTimeSeries() const
    {
//    std::cout << "Calling const GetConstTimeSeries()." << std::endl;
      if (fUpToDateDomain == Frequency)
      {
        // we have to update the cached time series data, so we need to cast away the constness of the this-object
        // this is ok since the change is "invisible" to the user, who always gets up-to-date data
        FFTDataContainer<C,T,F>* mutablethis = const_cast<FFTDataContainer<C,T,F>*>(this);
        mutablethis->UpdateTimeSeries(); // do the fft^-1
      }
      return *fTimeSeries;
    }

    /// read out the time series (write access)
    C<T>& GetTimeSeries()
    {
//    std::cout << "Calling non-const GetTimeSeries()." << std::endl;

      if (fUpToDateDomain == Frequency)
      {
        UpdateTimeSeries(); // do the fft^-1
      }
      fUpToDateDomain = Time;  // write access is possible, so have to assume data will change
      return *fTimeSeries;
    }

    /// read out the frequency spectrum (read access)
    const C<F>& GetConstFrequencySpectrum() const
    {
//    std::cout << "Calling const GetConstFrequencySpectrum()." << std::endl;
      if (fUpToDateDomain == Time)
      {
        // we have to update the cached frequency spectrum data, so we need to cast away the constness of the this-object
        // this is ok since the change is "invisible" to the user, who always gets up-to-date data
        FFTDataContainer<C,T,F>* mutablethis = const_cast<FFTDataContainer<C,T,F>*>(this);
        mutablethis->UpdateFrequencySpectrum(); // do the fft
      }
      return *fFrequencySpectrum;
    }

    /// read out the frequency spectrum (write access)
    C<F>& GetFrequencySpectrum()
    {
//    std::cout << "Calling non-const GetFrequencySpectrum()." << std::endl;

      if (fUpToDateDomain == Time)
      {
        UpdateFrequencySpectrum(); // do the fft
      }
      fUpToDateDomain = Frequency;  // write access is possible, so have to assume data will change
      return *fFrequencySpectrum;
    }

    /// delete contents of the contained traces    
    void Clear()
    {
      InitializedShadowPtr< C<T> > tempTimeSeries;               // create new initialized shadow pointers
      InitializedShadowPtr< C<F> > tempFrequencySpectrum;

      fTimeSeries.Swap(tempTimeSeries);                          // set the member data to the new (empty) objects
      fFrequencySpectrum.Swap(tempFrequencySpectrum);

      fUpToDateDomain = Both;                                    // no data there

      // when the temp shadow pointers now run out of scope, they will be automatically deleted and the memory deallocated    
    }

  private:

    enum eDomain { Time, Frequency, Both };
    eDomain fUpToDateDomain;

    unsigned int fNyquistZone;  // data corresponds to which Nyquist zone

    InitializedShadowPtr< C<T> > fTimeSeries;
    InitializedShadowPtr< C<F> > fFrequencySpectrum;

    friend class ShadowPtr< FFTDataContainer<C,T,F> >;  // necessary so that ShadowPtr can delete members of FFTDataContainer

    /// if the frequency domain data has been changed last, update the time domain data with an fft^-1
    void UpdateTimeSeries();
    /// if the time domain data has been changed last, update the frequency domain data with an fft
    void UpdateFrequencySpectrum();

  };

} // utl

#endif