#include"Spectral_Field.h"
#include"Constants.h"
#include"Wavefunction.h"
#include<iostream>
#include<complex>
#include<math.h>

Spectral_Field::Spectral_Field(){}
Spectral_Field::~Spectral_Field(){}

void Spectral_Field::set_omega(double E_start, double E_step, double E_end) {
	int size = (int)(E_end - E_start) / E_step;
	omega.init(size);

	for (int i = 0; i < size; i++) {
		omega[i] = (E_start + i * E_step) * const_e / const_hbar;
	}
}

void Spectral_Field::set_Eomega(Cvec Eomega) {
	if (Eomega.get_size() != omega.get_size()) {
		std::cout << "ERROR: Size of omega must be the same as of E(omega)" << std::endl;
		exit(0);
	}
	
	E_omega = Eomega;
}


void Spectral_Field::set_Eomega_gauss(double ampl, double omega_0, double width) {
	E_omega.init(omega.get_size());
	for (int i = 0; i < omega.get_size(); i++) {
		E_omega[i] = ampl * 1/sqrt(width) * exp(-pow(omega[i] - omega_0, 2) / width);
	}
}

void Spectral_Field::set_phi(double phi_0, double phi_1, double phi_2, double omega_0) {
	phi.init(omega.get_size());
	
	phi = phi_0 + (omega - omega_0) * phi_1 + (omega - omega_0) * (omega - omega_0) * phi_2;
	phi.save("phi.dat");
}

void Spectral_Field::calculate_Etime(Cvec tlist) {
	// tlist in [a.u.]
	E_time.init(tlist.get_size());
	double dw = real(omega[1] - omega[0]);
	Wavefunction E_help(E_omega.get_size(), dw);
	Wavefunction Exp_help(E_omega.get_size(), dw);
	
	for (int i = 0; i < E_omega.get_size(); i++) {
		E_help[i] = E_omega[i] * exp(Creal(0, -1) * phi[i]); 
		// das Minus muss hier stehen, da unten in der Berechnung von E_time das Skalarprodukt ausgewertet wird, 
		// dabei wird der erste Part (E_help) complex conjugiert
	}

	// SI units 
	for (int i = 0; i < tlist.get_size(); i++) {
		for (int j = 0; j < E_omega.get_size(); j++) {
			Exp_help[j] = exp(Creal(0, 1) * omega[j] * tlist[i] * au_T);
		}
		E_time[i] = 0.5 * 1 / sqrt(2*M_PI) * E_help.scalar(Exp_help);
		E_time[i] += conj(E_time[i]);
	}

	// back to [a.u.]
	E_time /= au_F;
}