////////////////

#include <iostream>
#include <iomanip>
using namespace std;

#include <boost/random/normal_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <time.h>       /* hack for initializing RNG */
#include <vector>
#include <fftw.h>

int main(){
  double R(1);
  double C(1/2./M_PI);
  double dt(.04);
  int nsteps(1024*1024);
  double tmax(nsteps * dt);
  double df(1. / tmax);
  double Vc(0);
  double kT(1e-10);
  double fN(1./dt/2.);
  cout << "RC:, " << R*C
        << ",  corner:, " << 1./R/C/2./M_PI
        << endl;
  cout  << "nsteps:, " << nsteps
        << setprecision(10) << ",  1/dt/df:, " << 1. / dt / df
        << ",  fN:, " << fN
        << endl;
  cout << "kT:, " << kT
        << endl;

  boost::variate_generator<boost::mt19937, boost::normal_distribution<> >
    generator(boost::mt19937(time(0)), boost::normal_distribution<>());

  double sum(0);
  double sumsq(0);
  double Vsum(0);
  double Vsumsq(0);
  vector<fftw_complex> Vlist(nsteps);
  vector<fftw_complex> Vhat(nsteps);
  fftw_plan myplan = fftw_create_plan(nsteps, FFTW_FORWARD, FFTW_ESTIMATE);
  for (int ii = 0; ii < nsteps; ii++){
    double rnd = generator();
    sum += rnd;
    sumsq += rnd*rnd;
    double nyquist = 1. / dt / 2.;      // circular frequency
    double Vin(rnd * sqrt(4. * kT * R * nyquist));
    double dv = Vc - Vin;       // will decay exponentially
#if 0
    Vc = Vin + dv * exp(-dt/R/C);
#else
    Vc = Vc + dv * expm1(-dt/R/C);
#endif
    Vsum += Vc;
    Vsumsq += Vc*Vc;
    Vlist[ii].re = Vc;
    Vlist[ii].im = 0;
  }
  cout << "sum:, " << sum
        << ",  sumsq:, " << sumsq
        << endl;
  double mean = sum / nsteps;
  double variance = sumsq / nsteps - mean*mean;
  cout << "variance:, " << variance
        << ",  stdev:, " << sqrt(variance)
        << endl;

// and now for the voltage:
  double ms_v = Vsumsq / nsteps;
//  double rms_v = sqrt(ms_v);
  double energy = 0.5 * C * ms_v;
  cout << "E: " << energy
       << "  E / (half kT): " << energy / (0.5 * kT)
       << endl;

///////////////////

// call the unnormalized version (from the library):
  fftw_one(myplan, &Vlist[0], &Vhat[0]);

// Then do the normalization.
// It won't pass the Parseval check without this:
  for (unsigned int ii=0; ii < Vhat.size(); ii++) {
    Vhat[ii].re *= dt;
    Vhat[ii].im *= dt;
  }

// At this point Vhat has units of volt-seconds,
// i.e. volts per unit bandwidth.

// Parseval, time domain:
  double parse1(Vsumsq * dt);

  double Vhat_sumsq(0);
  for (unsigned int ii = 0; ii < Vhat.size(); ii++) {
    double absq = Vhat[ii].re*Vhat[ii].re + Vhat[ii].im*Vhat[ii].im;
    Vhat_sumsq += absq;
  }
  double parse2(Vhat_sumsq * df);
  cout << "Parse1:, " << parse1
        << ",  Parse2:, " << parse2
        << ", " << (parse1 - parse2) / parse1
        << endl;

  if (0) for (int ii = 0; ii < 10; ii++) {
    cout
        << setw(20) << df * ii << ", "
        << setw(20) << Vhat[ii].re << ", "
        << setw(20) << Vhat[ii].im << ", "
        << endl;
  }

  int next = 0;
// we want more resolution near the end:
  int bmax = nsteps / 20 / 10;
  while (next < (15*nsteps)/20) {
    int batch = next / 10;    // 8 or 9 steps per octave, 26 per decade
    if (batch > bmax) batch = bmax;
    if (batch < 1) batch = 1;
    if (batch < 100) batch = 100;
    double sumsq(0);
    for (int ii = next; ii < next + batch; ii++) {
      double absq = Vhat[ii].re*Vhat[ii].re + Vhat[ii].im*Vhat[ii].im;
      sumsq += absq;
    }
    double avg(sumsq / (batch));

// The following factor of df can be explained in terms of
// dimensional analysis:  The properly-normalized ordinates
// in frequency space have dimensions of volt-seconds, i.e.
// voltage per unit bandwidth.  When we square that, we get
// volts squared per frequency squared.  We want simply
// volts squared per plain old unsquared frequency.
    avg *= df;

// To simplify comparison,
// normalize to the expected V^2 per unit bandwidth:
    double simpler (avg / (4. * kT * R));
    cout
        << setw(20) << df * (next + batch/2.) << ", "
        << setw(20) << simpler << ", "
        << setw(20) << batch << ", "
        << endl;
    next += batch;
  }
}