/*
    Install the GNU Scientific Library

    To use:
        #include "gsl.hpp"
        using namespace gsl;

    To compile:
        c++  prog.cpp  -lgsl  -lgslcblas  -lm
 */

#ifndef GSL_HPP
#define GSL_HPP

#include <iostream>
#include <string>

namespace gsl {
    struct GSL_Error {
        std::string name;
        GSL_Error(const char* c) : name(c) { }
        GSL_Error(std::string s) : name(s) { }
    };

    inline void error(const char* c)
    {
        std::cerr << c << std::endl;
        throw GSL_Error(c);
    }

} /* end namespace gsl */

#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#include <complex>
#include <vector>

namespace gsl {

    typedef std::vector< std::vector<double> >  Matrix;
    typedef std::vector<Matrix>                 Matrix3;
    typedef std::vector<Matrix3>                Matrix4;

    typedef std::vector< std::vector<std::complex<double> > > CMatrix;

} /* namespace gsl */

#include <gsl/gsl_fft_complex.h>

namespace gsl {

    static void fft(std::vector< std::complex<double> >& z)
    {
        size_t n = z.size();
        if (n == 0 || (n & (n-1)) != 0)
            gsl::error("fft n must be a power of 2");
        double *d = new double[2 * n];
        for (int j = 0; j < n; ++j) {
            d[0+2*j] = z[j].real();
            d[1+2*j] = z[j].imag();
        }        gsl_fft_complex_radix2_forward(d, 1, n);
        for (int j = 0; j < n; ++j)
            z[j] = (d[0+2*j], d[1+2*j]);
        delete [] d;
    }

    static void fft_inverse(std::vector< std::complex<double> >& z)
    {
        size_t n = z.size();
        if (n == 0 || (n & (n-1)) != 0)
            gsl::error("fft_inverse n must be a power of 2");
        double *d = new double[2 * n];
        for (int j = 0; j < n; ++j) {
            d[0+2*j] = z[j].real();
            d[1+2*j] = z[j].imag();
        }        gsl_fft_complex_radix2_inverse(d, 1, n);
        for (int j = 0; j < n; ++j)
            z[j] = (d[0+2*j], d[1+2*j]);
        delete [] d;
    }

    static std::vector<double> power(   // returns one-sided power spectrum
        const std::vector< std::complex<double> >& g)
    {
        size_t n = g.size();
        std::vector<double> P(1 + n/2);
        P[0] = std::norm(g[0]) / double(n);
        for (int j = 1; j < n/2; ++j)
            P[j] = (std::norm(g[j]) + std::norm(g[n-j])) / double(n);
        P[n/2] = std::norm(g[n/2]) / double(n);
        return P;
    }

} /* end namespace gsl */

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

namespace gsl {

    // default Mersenne Twister using mt19937 algorithm
    static gsl_rng *gslcpp_rng = gsl_rng_alloc(gsl_rng_default);

    static double ran_uniform()        // uniform deviate in [0,1)
    {
        return gsl_rng_uniform(gslcpp_rng);
    }

    static double ran_gaussian(double sigma=1.0)    // Gaussian deviate
    {
        return gsl_ran_gaussian(gslcpp_rng, sigma);
    }

} /* namespace gsl */

namespace gsl {

    void RK4_step(
        std::vector<double>& x,
        const double dt,
        std::vector<double> flow(std::vector<double>&) )
    {

    }

} /* namespace gsl */

#endif /* GSL_HPP */