#include"Pulse.h"
#include"Constants.h"
#include<iostream>
#include<math.h>


Pulse::Pulse(){}
Pulse::~Pulse(){}


Pulse::Pulse(double lambda_SI, double I0_SI, double t0_SI, double phi_SI, double FWHM_SI) {
	this->lambda = lambda_SI * 1e-9 / au_L;
	this->omega = 2 * M_PI * au_c / this->lambda;
	this->E0 = sqrt(2.0 * I0_SI * 1e4 / (const_epsilon_0 * const_c)) / au_F;
	this->t0 = t0_SI / (au_T * 1e15);
	this->FWHM = FWHM_SI / (au_T * 1e15);
	this->tau = this->FWHM / (2 * sqrt(log(2.0)));
	this->phi = phi_SI;
	this->Up = pow(au_e, 2) * pow(this->E0, 2) / (4 * au_m * pow(this->omega, 2)) * au_Ee;
}

double Pulse::get_Up() {
	return Up;
}

double Pulse::getField(double t_au) {
	return E0 * cos(omega * (t_au - t0) + phi) * exp(-pow(t_au - t0, 2.0) / (2.0 * pow(tau, 2.0)));
}

Cvec Pulse::getField_vp(Cvec t_au) {
	Cvec Efield(t_au.get_size());
	Cvec vectorpotential(t_au.get_size());
	double dt = t_au[2].real() - t_au[1].real();

	for (int i = 1; i < t_au.get_size(); i++) {
		vectorpotential[i] = E0 / omega * cos(omega * (t_au[i] - t0) + phi) * exp(-pow(t_au[i] - t0, 2.0) / (2.0 * pow(tau, 2.0)));
	}

	for (int i = 2; i < t_au.get_size() - 1; i++) {
		Efield[i] = -(vectorpotential[i + 1] - vectorpotential[i - 1]) / (2 * dt);
	}
	Efield[1] = -(vectorpotential[2] - vectorpotential[1]) / dt;
	Efield[t_au.get_size() - 1] = -(vectorpotential[t_au.get_size() - 1] - vectorpotential[t_au.get_size() - 2]) / dt;
	
	return Efield;
}

