///////////////////////
//
// Simulate RLC circuit,
// including the driving term
// specifically thermal noise driving term
//
// This is a series RLC circuit.
//
// The principal coordinate is the charge,
// i.e. the charge to the left of "here".
// Current "here" is then charge_dot.
// The voltage on the inductor is bigL charge_dot_dot
// The voltage on the resistor is bigR charge_dot
// The voltage on the capacitor is charge / bigC
//  which is the principal output of this program.
//
// Reference for why we use charge:
// "Quantum network theory"
// Phys. Rev. A 29, 1419 (1 March 1984)
// http://dx.doi.org/10.1103/PhysRevA.29.1419
//
// The incoming noise is white noise, with
// a spectrum that is flat from DC to daylight.
// This gets strongly filtered by the RLC circuit,
// so that the voltage on the capacitor is only
// weakly affected by the high-frequency noise.
//
//
// Right now I've got it set up to produce nice round numbers of
// order unity:
//  -- resonant frequency 1Hz
//  -- capacitance = 1 farad
//  -- oscillator quality factor "Q" = 10
//  -- really high temperature so that the noise voltage is 1 V
//    at low frequencies.  This is unreasonable as a temperature,
//    but there are plenty of non-thermal noise sources;  you
//    can always take white noise and amplify it.
//
// You could change the program to produce thermal voltages if
// you wanted, in which case you would need to start
// worring about the following:
//
// To model real thermal physics, the scale of the
// driving term in the simulation has to be:
//        sqrt(4 * kT * R / t__pt)      [1]
// where t__pt is the time interval between points,
// i.e. the timestep in the numerical integration.
// This makes sense on dimensional grounds.
//
// More importantly, it scales the right way if
// you consider what would happen if you did a
// "batch" update consisting of several such impulses,
// then summed over batches.
// The integral dt over batches would have a
// longer dt, but the integrand would be smaller
// by a factor of sqrt(dt), so the net effect of
// the batch scales like sqrt(dt) in the numerator,
// which is the right answer for a random walk.
// Also consider that when you take the FFT, 1/dt
// *is* the bandwidth, so equation [1] is consistent
// with the textbook formula
//        sqrt(4 * kT * R * B)          [2]
//
// Note that the ordinate of Vhat has dimensions
// of volt-seconds (not plain volts).
// Reference: https://www.av8n.com/physics/fourier-refined.htm

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <boost/random.hpp>
#include <boost/random/normal_distribution.hpp>
#include <fftw.h>
#include <inttypes.h>  /* for uint64_t, needed on some systems */
#include "strlib.h"

using namespace std;

uint64_t get_seed() {
  union{
    uint64_t seed;
    char buf[1];
  };
  std::ifstream urandom;
  urandom.open("/dev/urandom");
  urandom.read(buf, sizeof(seed));
  urandom.close();
  return seed;
}

double sumsq(vector<fftw_complex> const vec){
  double rslt(0);
  for (unsigned int ii=0; ii < vec.size(); ii++) {
    rslt += pow(vec[ii].re, 2) + pow(vec[ii].im, 2);
  }
  return rslt;
}

double rms(vector<fftw_complex> const vec){
  return sqrt(sumsq(vec) / vec.size());
}

string const Q("\"");           // unduly ugly;  needs less-ugly synonym

// properly normalized FFT
void fftw_one(fftw_plan& myplan, double const t__pt,
  vector<fftw_complex>& vec,    // should be const, but fftw doesn't support it
  vector<fftw_complex>& vhat    // (not const anyway)
){
  fftw_one(myplan, &vec[0], &vhat[0]);

// do the normalization
  for (unsigned int ii=0; ii < vec.size(); ii++) {
    vhat[ii].re *= t__pt;
    vhat[ii].im *= t__pt;
  }
}

int main(){
// Some miscellaneous constants:
  fftw_complex zip;  zip.re = zip.im = 0;
  double const pi(M_PI);
  double const boltzmann(1.3806488e-23); // J/K
  double temperature(300);               // kelvin
  double kT(boltzmann * temperature);

// Some basic properties of the simulation:
// t__cyc is prounounced time *per* cycle
// t__pt  is pronounced time *per* point
  int ncycles(128);
  int p__cyc(128);                 // points per cycle
  int ptot(p__cyc * ncycles);      // points total
  double t__cyc(1);                // time per cycle in seconds
  double ttot(t__cyc * ncycles);   // total time in fft buffer
  double t__pt(t__cyc / p__cyc);   // time per point
  double f__pt(1. / ttot);         // frequency per point
  double ftot(f__pt * ptot);
  double ftot2(1. / t__pt);
  if (0) {
    cerr << "t__pt:  " << t__pt  << endl;
    cerr << "t__cyc: " << t__cyc << endl;
    cerr << "ttot:   " << ttot   << endl;
    cerr << "ptot:   " << ptot   << endl;
    cerr << "f__pt:  " << f__pt  << endl;
    cerr << "ftot:   " << ftot   << endl;
    cerr << "ftot2:  " << ftot2
        << "  check1: " << ftot / ftot2 << endl;
  }

// Some properties of the circuit:
  double natFreq_Hz(1. / t__cyc);     // natural freq in Hz
  double natFreq_rps(2. * pi * natFreq_Hz);     // radians per second
  double bigC(1);        // in Farads
  double bigL(1. / bigC / natFreq_rps / natFreq_rps);
  if (1) {
    // check;  should be inverse of t__cyc:
    double xfreq = 1 / sqrt(bigL * bigC) / 2. / pi;
    cerr << "xfreq: " << xfreq
        << "  check1: " << xfreq * t__cyc
        << endl;
  }
  double Qual(10);              // Q as in quality (not Q as in charge)
  double bigZ(sqrt(bigL / bigC));       // reactive impedance at resonance
  double bigR(bigZ / Qual);
  double vscale(sqrt(4 * kT * bigR / t__pt));
  if (1) {
    // This gives nice pretty round numbers for debugging:
    vscale = 1.;
  }

// Cobble up the drive voltage waveform:

  vector<fftw_complex> drive_voltage(ptot);
  vector<fftw_complex> drive_vhat(ptot);

  boost::mt19937 rng(get_seed());
  boost::normal_distribution<> nd(0.0, 1.0);
  boost::variate_generator<boost::mt19937&, boost::normal_distribution<> >
        gaussian(rng, nd);

  for (int ii=0; ii < ptot; ii++) {
    drive_voltage[ii].re = gaussian() * vscale;
    drive_voltage[ii].im = 0;
  }
// RMS should be very close to 1.
  if (0) cerr << "RMS: " << rms(drive_voltage) << endl;

  fftw_plan myplan = fftw_create_plan(ptot, FFTW_FORWARD, FFTW_ESTIMATE);

  fftw_one(myplan, t__pt, drive_voltage, drive_vhat);

// Print out drive voltage.
  if (0) for (int ii=0; ii < ptot; ii++) {
    cout << setw(10) << ii
         << setw(15) << drive_vhat[ii].re
         << endl;
  }

// check that Parseval's theorem is still a theorem:
  cerr << "time norm^2: " << sumsq(drive_voltage) * t__pt
        << " volt^2 seconds" << endl;

  cerr << "freq norm^2: " << sumsq(drive_vhat) * f__pt
        << " volt^2 seconds" << endl;

// At this point we have the drive voltage set up.

// do the numerical integration
  double charge(0);
  double charge_dot(0);
  double charge_dot_dot(0);
  vector<fftw_complex> rslt_v(ptot);
  vector<fftw_complex> rslt_vhat(ptot);
  for (int ii=0; ii < ptot; ii++) {
    double foo = drive_voltage[ii].re - bigR * charge_dot - charge / bigC;
    charge_dot_dot = foo / bigL;
    charge_dot += charge_dot_dot * t__pt;
    charge += charge_dot * t__pt;
    rslt_v[ii].re = charge / bigC;
    rslt_v[ii].im = 0;
  }

// Write results of numerical integration to a file.
  ofstream ouch;
  string fname("harmonic-noise-v.csv");
  ouch.open(fname.c_str());
  if (!ouch.good()) {
    cerr << "Cannot open output file '" << fname << "'" << endl;
    cerr << Strerror() << endl;
    exit(1);
  }

  ouch << setw(15) << Q+"time"+Q << ","
        << setw(15) << Q+"Voltage"+Q << ","
        << endl;

  for (int ii=0; ii < ptot; ii++) {
    ouch << setw(15) << ii * t__pt << ","
         << setw(15) << rslt_v[ii].re << ","
         << endl;
  }
  ouch.close();

// Do the big FFT.
  fftw_one(myplan, t__pt, rslt_v, rslt_vhat);

// Write FFT results to a file.
  string fname_vhat("harmonic-noise-vhat.csv");
  ouch.open(fname_vhat.c_str());
  if (!ouch.good()) {
    cerr << "Cannot open output file '" << fname_vhat << "'" << endl;
    cerr << Strerror() << endl;
    exit(1);
  }

  ouch << setw(15) << Q+"freq"+Q << ","
        << setw(15) << Q+"drive_Vhat^2"+Q << ","
        << setw(15) << Q+"rslt_Vhat^2"+Q << ","
        << endl;

  for (int ii=0; ii < ptot/2; ii++) {
    double drive_vhat2 = pow(drive_vhat[ii].re, 2) + pow(drive_vhat[ii].im, 2);
    double rslt_vhat2 = pow(rslt_vhat[ii].re, 2) + pow(rslt_vhat[ii].im, 2);
    ouch << setw(15) << ii * f__pt << ","
         << setw(15) << drive_vhat2 << ","
         << setw(15) << rslt_vhat2 << ","
         << endl;
  }
  ouch.close();
  fftw_destroy_plan(myplan);
}


// kT has units of energy, i.e. kg m^2 / s^2
// V^2 / R has units of energy per unit time
// R has units if voltage^2 time / energy
// kT R / dt has units of voltage^2