#ifndef _utl_ComplexLgAmpPhase_h_
#define _utl_ComplexLgAmpPhase_h_

#include <cmath>
#include <complex>
#include <iostream>
#include <utl/AugerException.h>

/**
  \file
  A class to store complex numbers which are internally represented by log10(amplitude) and phase
  (and not real and imaginary part). This is needed for some operations like interpolation and pow().
  Only developed as far as needed for application inside radio, feel free to extend it ... ;-)

  \author Tim Huege
  \ingroup math
*/


namespace utl {

  /**
    \class ComplexLgAmpPhase ComplexLgAmpPhase.h "utl/ComplexLgAmpPhase.h"
    \brief
    A class to store complex numbers which are internally represented by log10(amplitude) and phase
    (and not real and imaginary part). This is needed for some operations like interpolation and pow().

    \author Tim Huege
  */

  class ComplexLgAmpPhase {

  public:
    /// Standard constructor initializing to zero
    ComplexLgAmpPhase() : fLgAmplitude(0), fPhase(0) { }

    /// Constructor using log10(amplitude) and phase (in radians)
    ComplexLgAmpPhase(const double lgamp, const double phase = 0)
      : fLgAmplitude(lgamp), fPhase(phase) { }

    /// Get log10(amplitude)
    double GetLgAmplitude() const
    { return fLgAmplitude; }

    /// Get the amplitude
    double GetAmplitude() const
    { return std::pow(10.0,fLgAmplitude); }

    /// Get the phase in radians, can be bigger than 2*pi
    double GetPhase() const
    { return fPhase; }

    /// Get the complex number in standard notation; this clips the phase to 2*pi!
    std::complex<double> GetComplex() const
    { return std::pow(10.0,fLgAmplitude) * std::exp(std::complex<double>(0, 1) * fPhase); }

    /// Get the real part
    double GetReal() const
    { return GetComplex().real(); }

    /// Get the imaginary part
    double GetImag() const
    { return GetComplex().imag(); }

    // careful, it is not possible to implement addition or subtraction without losing the extra phase information!

    ComplexLgAmpPhase
    operator*(const ComplexLgAmpPhase& rhs) const
    {
      double newLgAmplitude = this->fLgAmplitude + rhs.fLgAmplitude;
      double newPhase = this->fPhase + rhs.fPhase; 
      return ComplexLgAmpPhase(newLgAmplitude,newPhase);
    }
    
    ComplexLgAmpPhase
    operator/(const ComplexLgAmpPhase& rhs) const
    {
      double newLgAmplitude = this->fLgAmplitude - rhs.fLgAmplitude;
      double newPhase = this->fPhase - rhs.fPhase;
      return ComplexLgAmpPhase(newLgAmplitude,newPhase);
    }

    ComplexLgAmpPhase&
    operator*=(const ComplexLgAmpPhase& rhs)
    {
      fLgAmplitude += rhs.fLgAmplitude;
      fPhase += rhs.fPhase;
      return *this;
    }

    ComplexLgAmpPhase&
    operator/=(const ComplexLgAmpPhase& rhs)
    {
      fLgAmplitude -= rhs.fLgAmplitude;
      fPhase -= rhs.fPhase;
      return *this;
    }

    // should this operator use Test<CloseTo> instead of operator==?
    bool operator==(const ComplexLgAmpPhase& rhs) const
    { return fLgAmplitude == rhs.fLgAmplitude && fPhase == rhs.fPhase; }

    ComplexLgAmpPhase&
    Pow(double power)
    {
      fLgAmplitude *= power;
      fPhase *= power;
      return *this;
    }

  private:
    double fLgAmplitude;
    double fPhase;

  };

  /// Do linear interpolation of two ComplexLgAmpPhases in log10(amplitude) and phase
  inline
  ComplexLgAmpPhase
  Interpolate(double x1, const ComplexLgAmpPhase& y1,
              double x2, const ComplexLgAmpPhase& y2,
              double x)
  {
    const double a1 = y1.GetLgAmplitude();
    const double p1 = y1.GetPhase();
    const double a2 = y2.GetLgAmplitude();
    const double p2 = y2.GetPhase();
    const double frac = (x - x1) / (x2 - x1);
    const double anew = a1 + frac * (a2 - a1);
    const double pnew = p1 + frac * (p2 - p1);
    return ComplexLgAmpPhase(anew, pnew);
  }

} // utl

#endif