#ifndef noise_h
#define noise_h
// include files for gaussian deviates
#include <cmath>
#include "nr3.h"
#include "ran.h"
#include "gamma.h"
#include "deviates.h"

int idum=87654321; // start seed
Ran ran(idum); // random generator for noise
Normaldev gaus(0,1,idum); // gaussian deviate for noise

// noise_psd returns the amplitude of the noise, i.e. the square root of the
// power spectral density, in 1/sqrt(Hz)
double noise_psd(double freq) {
	if (freq>2048) freq=4096-freq; // Nyquist frequency at 2048 Hz
	if (freq<1) return 1e-14;
	double psd= (1e-25*freq+1e-9/pow(freq,8)+5e-20/freq);
	if (psd>1e-14) return 1e-14;
	return psd;
};

Complex noise(double freq) {
// use up to Nyquist.
	Complex a;
	double norm=sqrt(0.5)*noise_psd(freq);
	a.real(norm*gaus.dev());
	a.imag(norm*gaus.dev());
	return a;
};

// filters etc need to be normalized. keep track of the overall factor,
// this avoids needing to multiply all entries
double power(const VecComplex& vec) {
	// normalize a complex vector
	double fac=0;
	for (int i=0; i<vec.size(); i++) fac+=norm(vec[i]); // norm(complex<double>) gives the power
	return fac;
}

#endif