#ifndef TDSE_H
#define TDSE_H
#include<complex>
#include<fstream>
#include"Operator.h"
#include"Wavefunction.h"
#include"Potential.h"
#include"Cvec.h"
#include"Pulse.h"
#include"Ctridiag.h"
#include"TrackManager.h"
#include"Save_Array.h"
#include"Pulsetrain.h"
#include"Spectral_Field.h"
#include"mpimanager.h"
#include<vector>

typedef std::complex<double> Creal;

class TDSE {
public:
	// Basics
	Wavefunction Psi;						// actual wavefunction
	Wavefunction Psi_start;					// start wf for propagation
	Ctridiag H_0;							// Hamiltonian (kin + pot)
	Ctridiag H_kin;							// Only kin Hamiltonian
	Cvec xlist;								// Box for simulation
	Cvec tlist;								// Time of propagation
	string path;							// path for saving
	Wavefunction help_wave;					// help WF
	Ctridiag help_tridi;					// help tridiag matrix
	Cvec help_cvec;							// help vector
	Wavefunction wave_interm;				// intermediate wave function

	// Potentials & Laser-Fields
	Ctridiag V_tri;							// Potential in tridiag form
	Potential V;							// Potential in vector form
	Potential V_imag_vec;					// Imaginary boundary potential (absorbing BC)
	Cvec E_las = 0;							// Laser field (all together, time-dependent)
	Potential V_las;						// Potential of laser fields (time-dependent)
	std::vector<Pulse> Pulse_prop;			// Vector of Pulses used in Propagation
	std::vector<Pulse> Pulse_prop_vp;		// Vector of Pulses used in Propagation (defined via vector potential)
	Pulsetrain Pulsetrain_prop;				// Pulsetrain for propagation (optional)
	Cvec Efield_input = 0;					// Inported Field (optional)
	Cvec E_field_vp = 0;					// Laserfield of pulses defined via vector potential
	Spectral_Field Spectral;				// Laser field in spectral decomposition (optional)
	Cvec gamma = 0;							// Enhancement factor
	double Up = 0;

	// Imaginary time evolution
	std::vector<double> Energylevels;		// List of energylevels of different states with imaginary time propagation
	std::vector<Wavefunction> Psi_states;	// List of WF of different states 
	std::vector<Wavefunction> Psi_states_GS;// List of WF of different states with imaginary time propagation (not changed in propagation)
	int state_max = -1;						// #states calculated in Imag. time prop.

	// Propagation
	Save_Array Psitime;						// Time- and Space-resolved WFs (in Propagation)
	Save_Array Popul;						// Timeresolved population for states calculated in Imag. time prop.
	Save_Array Potential_save;				// Time- and space-resolved potential (all together)
	Save_Array Spec_evolution;				// Time-resolved spectrum
	Operator Am;
	Operator Ap;

	// Spectrum (Window-Operator)
	Cvec Elist = 0;							// List of test energies for spectrum calculation (Window Operator)
	Cvec spec = 0;							// Spectrum dN/dE (Window)
	int write_spec_evol = 0;

	// Track stuff
	int track_var = 0;
	TrackManager tracker;					// Track-Manager for Variables in Propagation
	struct Tracks {
		double time = 1;
		double expectaionvalue_x = 1;
		double expectaionvalue_E = 1;
		double expectaionvalue_Ekin = 1;
		double expectaionvalue_p = 1;
		double expt_x_left = 1;				// Only x < x0 (in propagation)
		double expt_x_right = 1;			// Only x > x0
		double norm_left = 1;
		double norm_right = 1;
	}tracks;

	Ctridiag P_stencil;						// Momentum Operator
	//Numerov Boost:
	Operator Delta;
	Operator M;

	double x_right = 0.0;					// right edge of potential box, needed for gamma(x) and calculating spectra only outside of box

	MpiManager* mpi;


public:

	TDSE(double& dxin);
	~TDSE();

	// In: Å, Out: a.u.
	void init_xlist(double start, double step, double end);

	// In: au, Out: a.u
	void init_xlist_au(double start, double step, double end);

	// In: au, Out: a.u
	void init_xlist_au(Cvec x_list);

	// In: fs, Out: a.u.
	void init_tlist(double start, double step, double end);

	// In: fs, Out: a.u.
	void init_tlist(Cvec t_input);

	// Set Psi random, init H0, V, V_tri
	void init_tdse();

	// Initialize tracker -> tracks variables over propagation time
	void init_tracker(std::string filename);

	// Calculates the variables to be tracked
	void calculate_tracks_variables(Wavefunction psi_prop);
	
	// init Coulomb
	void add_Coulomb(double x0, int Z, double a);

	// init Coulomb
	void add_Yukawa(double x0, int Z, double a);

	// init box potential
	void add_box(double depth, double length, double x0);

	// init box with Fermi edge right
	void add_fermibox(double depth, double s, double length, double x0);

	// init benutzerdefinierte Box
	void set_box_benutzer(double depth, double s, double length, double x0);

	// set pot to zero and init fermi box
	void set_fermibox(double depth, double s, double length, double x0);

	// add potential to Hamiltonian
	void set_Hamiltonian();

	// Numerov-Boost: Set propagation matrices
	void set_AmAp(Creal timestep);

	// Numerov-Boost: calculate energy for given wavefunction
	double calc_Energy(Wavefunction psi);

	// run imaginary time propagation
	void imaginary_time_prop(const int max_state, double time_step);
	void imaginary_time_prop_Vmin(const int max_state, double Vmin);

	// imaginary time propagation - calc states parallel
	void imaginary_time_prop_parallel(int all_states);

	// find the groundstate with nearest energy to input energy
	double find_state(double Energy);

	// find Box potential with corresponding workfunction & fermi level (only one state)
	void find_box_pot(double W, double Ef);

	// for import of stuff (state_max)
	void read_param(string path);

	// import potential and groundstate
	void import_GS(string in_path, string save);

	// set start wavefunction for propagation
	void set_psi_start(std::vector<std::vector<double>> state_number);

	// Add Laserpulse to vector
	void add_pulse(double lambda_SI, double I0_SI, double t0_SI, double phi_SI, double FWHM_SI);

	// Add pulse with definiton via vector potential
	void add_pulse_vp(double lambda_SI, double I0_SI, double t0_SI, double phi_SI, double FWHM_in, int use_SC);

	// Add Pulsetrain to propagation
	void add_pulsetrain(double lambda_SI, double I0_SI, std::vector<double> t0_SI, double phi_SI, double FWHM_SI, double t0_env, double FWHM_env);

	// Inport Field (tlist in a.u., Field in a.u., interpolation)
	void input_field(string filename);

	// Set enhancement factor gamma (Input [Å], output [a.u.])
	void set_gamma(double gamma0, double lambda_nf);

	// Calculate cummulative integral of enhancement factor gamma
	Cvec calculate_int_gamma();

	// run propagation with extern field
	void propagation(Wavefunction psi_start, int write_timesteps, double write_x, int write_t, string input_path);

	// simple test function
	void propagation_test();

	// run propagation with extern field (IN: a.u.)
	void propagation_au(Wavefunction psi_start, int write_timesteps, double write_x, int write_t, string input_path);

	// ouprojection of psi_out of Psi
	void outprojection(Wavefunction psi_out);

	// calculate spetrum with Window_Operator (last entry == 1 -> spectrum for x > x0)
	void calc_spec(Wavefunction psi_start, double Energy_start, double Energy_step, double Energy_end, int order, double gamma, int half_spec);
	void calc_spec(Wavefunction psi_start, Cvec Energy, int order, double gamma, int half_spec);

	void save_param(string input_path);

	// save all 
	void save_all();

};

#endif
