spectral.h

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#ifndef __spectral__
#define __spectral__

#include <cmath>

/// Namespace containing entire library.

namespace spectral
{

  /// Real number type definition.
  /// The entire library uses this type as the real number type.  
  /// By default, it is set to double, but it can also be
  /// set to float or long double by modifying this typedef and recompiling.
  
  typedef double real; 
  
  /// Complex number type definition.
  /// Emulates most of the STL complex<real> behavior but with public real and imaginary components, 
  /// entirely inline for performance.  Included are functions and operators for the complex type.
  /// \sa real
  
  class complex 
  {
  public:
    complex(real x=0.0,real y=0.0) : re(x),im(y) {};
    real re;
    real im;
  };
  
  /// The square of the absolute value \f$ x^2 \f$.

  inline real norm(real x) 
  { 
    return(x*x); 
  } 
  
  /// The square of the magnitude \f$ |z|^2 \f$.

  inline real norm(complex z) 
  { 
    return(z.re*z.re+z.im*z.im); 
  }
  
  /// The complex magnitude operator \f$ |z| \f$.

  inline real abs(complex z) 
  { 
    return(std::sqrt(norm(z))); 
  }
  
  /// The complex conjugate operator \f$ \overline{z} \f$.

  inline complex conj(complex z) 
  { 
    return(complex(z.re,-z.im)); 
  }
  
  /// The complex multiplication operator.

  inline complex operator*(complex a,complex b) 
  { 
    return(complex(a.re*b.re-a.im*b.im,a.re*b.im+a.im*b.re)); 
  }

  /// The complex addition operator.
    
  inline complex operator+(complex a,complex b) 
  { 
    return(complex(a.re+b.re,a.im+b.im)); 
  }
    
  /// The complex subtraction operator.

  inline complex operator-(complex a,complex b) 
  { 
    return(complex(a.re-b.re,a.im-b.im)); 
  }
    
  /// The complex-real multiplication operator.

  inline complex operator*(complex a,real x) 
  { 
    return(complex(a.re*x,a.im*x)); 
  }
  
  /// The real-complex multiplication operator.

  inline complex operator*(real x,complex a) 
  { 
    return(complex(a.re*x,a.im*x)); 
  } 
    
  /// The complex-real division operator.

  inline complex operator/(complex a,real b) 
  { 
    return(complex(a.re/b,a.im/b)); 
  }
    
  /// The complex division operator.

  inline complex operator/(complex a,complex b) 
  { 
    return(a*conj(b)/norm(b)); 
  }
    
  /// The complex addition-assignment operator.

  inline complex operator+=(complex &a,complex b) 
  { 
    a=a+b; 
    return a; 
  }
    
  /// The complex subtraction-assignment operator.

  inline complex operator-=(complex &a,complex b) 
  { 
    a=a-b; 
    return a; 
  }
    
  /// The complex-real multiplication-assignment operator.

  inline complex operator*=(complex &a,real b) 
  { 
    a=a*b; 
    return a; 
  }

  /// The complex multiplication-assignment operator.
  
  inline complex operator*=(complex &a,complex b) 
  { 
    a=a*b; 
    return a; 
  }
  
  /// The unary complex negation operator $-z$.

  inline complex operator-(complex z) 
  { 
    return(complex(-z.re,-z.im)); 
  }
  
  /// The complex exponetial function $e^z$.

  inline complex exp(complex z) 
  { 
    real r=std::exp(z.re); 
    return(complex(r*std::cos(z.im),r*std::sin(z.im))); 
  }
  
  /// Generic binary maximum function template.

  template <typename T> inline T max(T a,T b) 
  { 
    return((a>b)?(a):(b)); 
  }
    
  /// Generic binary minimum function template.

  template <typename T> inline T min(T a,T b) 
  { 
    return((a<b)?(a):(b)); 
  }
    
  /// Wavenumber translation function.
  /// Returns real-valued number of wavenumber corresponding to array index for Fourier coefficients.
  /// The Fourier coefficients of a sequence of length n are stored as 
  /// \f[ k=\{0,1,\ldots ,\lfloor n/2 \rfloor -1,- \lfloor n/2 \rfloor ,\ldots ,-1\} \f]
  /// corresponding to the index space 
  /// \f[ i=0,\ldots ,n-1. \f]
  /// \param n length of spectral coefficient sequence
  /// \param i wavenumber index   
  /// \return real-valued wavenumber corresponding to index
  
  inline real wavenumber(int n,int i) 
  { 
    return((i<n/2)?(real)i:real(i-n)); 
  }
  
  /// A long double representation of \f$ \pi \f$.

  const real pi   =  3.141592653589793238462643383279L;

  /// A long double representation of \f$ \sqrt{2} \f$.

  const real root2   =  1.41421356237309504880168872420969L;

  /// A long double representation of \f$ \sqrt{3} \f$.

  const real root3   =  1.732050807568877293527446341505L;
  
  /// Exception for illegal parameter error.  Thrown from methods and object constructors
  /// when an illegal parameter that is out of range is passed.

  class bad_parameter {};
}
#endif

// Local Variables:
// mode:C++
// End: