#include<iostream>
#include"Cvec.h"
#include"Ctridiag.h"
#include"Wavefunction.h"
#include"Potential.h"
#include"Constants.h"
#include"Pulse.h"
#include"TDSE.h"
#include"TrackManager.h"


#define _USE_MATH_DEFINES
#include<math.h>
#include<time.h>
#include<stdio.h>
#include<stdlib.h>
#include<fstream>
#include<sstream>
#include<array>
#include<vector>
#include<filesystem>

using namespace std;
namespace fs = std::filesystem;


int main() {


	string path_output = "D:/Projects/TDSE/";

	// ----- Speicherort -----
	int run = 50;
	ostringstream oss;
	oss << path_output << "Output/Test" << run << "/";
	string path1 = oss.str();
	fs::create_directories(path1);

	// Tracken von Variablen über der Zeit? (ja - 1)
	int track = 0;

	// ----- Speicher der Wellenfunktion in jedem Zeitschritt? (ja - 1)
	int write_spec = 1;

	// ----- Import von Feldern? (nein - 0) 
	string input = "./Input/Chirp_test.dat";
	input = "0";

	// Import vom Grundzustand
	ostringstream in;
	in << path_output << "Output/GS/GS5breite50_4/";
	string input_GS = in.str();

	// Vergleich mit TDDFT mit importiertem Potential

	//double dx_au = 0.05;
	//TDSE tdse(dx_au);
	//
	////tdse.init_xlist_au(-1000.0, dx_au, 1000.0);
	////tdse.init_tdse();
	////tdse.path = oss.str();
	////
	////double V0 = 0.75;
	////double width = 50;
	////tdse.add_fermibox(-V0, 2, width, 0);
	////tdse.set_Hamiltonian();
	////tdse.imaginary_time_prop_Vmin(1, -1.5 * tdse.V.minimum()); // use 1.5*V for many box states -> ensure that first state found is groundstate 
	//
	//tdse.import_GS(input_GS, path1);
	//tdse.x_right = 25;
	//
	//vector<vector<double>> state_numbers{};
	//tdse.set_psi_start(state_numbers);
	//
	//tdse.add_pulse(500, 2e14, 0, 0, 5);
	//tdse.set_gamma(1, 250);
	//tdse.propagation(-15.0, 0.0001, 10.0, tdse.Psi_start, write_spec, input);
	//
	//for (auto& psi : tdse.Psi_states) {
	//	tdse.outprojection(psi);
	//}
	////tdse.outprojection(tdse.Psi_start);
	//
	//tdse.calc_spec(tdse.Psi, -30.0, 0.01, 75.0, 2, 0.0025, 1); 
	//
	//cout << "Save data..." << endl;
	//tdse.save_all();

	// ---------- Coulomb-Potential & Imaginärzeitentwicklung & Propagation & Spektrum ----------

	//double dx = 0.005;
	//double dx_au = dx * 1e-10 / au_L;
	//
	//for (double tau = 0; tau <= 0; tau += 0.1) {
	//	ostringstream oss1;
	//	oss1 << oss.str();
	//	oss1 << "tau" << tau << "/";
	//	cout << "tau = " << tau << " fs" << endl;
	//	path1 = oss1.str();
	//	fs::create_directory(path1);
	//	
	//	TDSE tdse(dx_au);
	//	tdse.init_xlist(-1000.0, dx, 2000);
	//
	//	tdse.init_tdse();
	//
	//	tdse.path = oss1.str();
	//
	//	tdse.add_Coulomb(0.0, 1, 2.0);
	//	tdse.set_Hamiltonian();
	//	tdse.imaginary_time_prop(1, 1);
	//
	//	vector<vector<double>> state_numbers{ {1.0, 0.0} };
	//	tdse.set_psi_start(state_numbers);
	//
	//	/*vector<double> t0;
	//	for (double i = -36; i <= 36; i += 1.5) t0.push_back(i);
	//	tdse.add_pulsetrain(50, 1e11, t0, 0.0, 0.5, 0.0, 17);*/
	//
	//	/*tdse.add_pulse(43.3511, 1e8, 0.0, 0.0, 5);
	//	tdse.add_pulse(36.9551, 1e8, 0.0, 0.0, 5);
	//	tdse.add_pulse(500, 1e6, tau, 0.0, 5);*/
	//
	//	/*tdse.add_pulse(40, 1e10, 0.0, 0.0, 0.5);
	//	tdse.add_pulse(800.0, 1e10, tau, 0.0, 5);*/
	//
	//	tdse.add_pulse(800.0, 1e14, 0.0, 0.0, 5);
	//	//tdse.add_pulse(50.0, 1e9, 0.0, 0.0, 5);
	//	
	//	/*double w0 = 11 * const_e / const_hbar;
	//	tdse.Spectral.set_omega(0, 0.001, 25);
	//	tdse.Spectral.set_Eomega_gauss(1e7, w0, 1e30);
	//	tdse.Spectral.set_phi(0, 0, 5*1e-31, w0);*/
	//
	//	if (tau < 0) {
	//		tdse.propagation(tau - 15.0, 0.005, 15.0, tdse.Psi_start, write_spec, input);
	//	}
	//	else if (tau > 0) {
	//		tdse.propagation(-15, 0.005, tau + 15.0, tdse.Psi_start, write_spec, input);
	//	}
	//	else {
	//		tdse.propagation(-10, 0.001, 10.0, tdse.Psi_start, write_spec, input);
	//	}
	//	tdse.outprojection(tdse.Psi_start);
	//	tdse.calc_spec(tdse.Psi, 0.0, 0.05, 80.0, 3, 0.0025, 1);
	//	tdse.save_all();
	//}

	/* ---------- Box-Potential & Imaginärzeitentwicklung & Propagation & Spektrum ---------- */

	double dx_au = 0.01;
	
	 /*for (double lambda_exp = 0.9; lambda_exp <= 4; lambda_exp += 0.1) {
		double lambda = pow(10.0, lambda_exp);
		cout << "lambda = " << lambda << endl;
		ostringstream oss1;
		oss1 << oss.str();
		oss1 << "Tip_lambda_1e" << lambda_exp <<"/";
		path1 = oss1.str();
		fs::create_directory(path1);*/
	
		/*ostringstream oss1;
		oss1 << oss.str();
		oss1 << "Box_30_test/";
		path1 = oss1.str();
		fs::create_directory(path1);*/
	
	
	 /*für Schleifen: Zustände vorher berechnen und übernehmen*/
	TDSE start_tdse(dx_au);
	start_tdse.init_xlist_au(-1000.0, dx_au, 5000.0);
	start_tdse.init_tdse();
	start_tdse.path = oss.str();
	
	//tdse.find_box_pot(4.3, 7);
	
	double V0 = 11.3 / au_Ee;
	double width = 0.980469 / (1e10 * au_L);
	//double width = 150;
	start_tdse.add_box(-V0, width, -width / 2);
	//start_tdse.add_fermibox(-V0, 2, width, -width / 2);
	//start_tdse.set_box_benutzer(-V0, 3, width, -width / 2);
	start_tdse.set_Hamiltonian();
	start_tdse.imaginary_time_prop_Vmin(1, V0); // use 1.5*V for many box states -> ensure that first state found is groundstate
	
	
	//for (double phase = 0*M_PI; phase <= 0*M_PI; phase += 0.05*M_PI) {
	//	cout << "phi = " << round(phase * 100.0) / 100.0 << endl;
	//	ostringstream oss1;
	//	oss1 << oss.str();
	//	oss1 << "phase_" << round(phase * 100.0) / 100.0 << "/";
	
	for (double tau = 0; tau <= 0; tau += 0.25) {
		cout << "tau = " << tau << endl;
		ostringstream oss1;
		oss1 << oss.str();
		oss1 << "tau_" << tau << "/";
		
		path1 = oss1.str();
		fs::create_directory(path1);
	
		double lambda = 250;
		double phase = 0*M_PI;
		//double tau = -0.5;
	
		TDSE tdse(dx_au);
		tdse.init_xlist_au(-1000.0, dx_au, 5000.0);
	
		tdse.init_tdse();
		tdse.path = oss1.str();
	
		tdse.track_var = track;
		
		tdse.V = start_tdse.V;
		tdse.set_Hamiltonian();
	
		tdse.Energylevels = start_tdse.Energylevels;
		tdse.Psi_states = start_tdse.Psi_states;
		tdse.state_max = start_tdse.state_max;
		//tdse.imaginary_time_prop_Vmin(50, 1.2*V0); // use 1.5*V for many box states -> ensure that first state found is groundstate
		
		//double start_state = tdse.find_state(4.29);
		//cout << start_state << " mit E = " << tdse.Energylevels[start_state] << endl;
		vector<vector<double>> state_numbers{};
		/*int num_startstates = 5;
		for (double i = 0.0; i < num_startstates; i++) {
			state_numbers.push_back({1.0 / num_startstates, i});
		}*/
	
		//vector<vector<double>> state_numbers{ {1.0, start_state} };
		tdse.set_psi_start(state_numbers);
	
		//tdse.add_pulse(1500, 1e11, 0, phase, 20);
		tdse.add_pulse(800, 1e12, 0, phase, 20);
	
		//tdse.add_SC_pulse(800, 1e12, 0, phase);
		//tdse.add_SC_pulse(800, 1e10, tau, phase);
	
		tdse.set_gamma(7, lambda);
		tdse.propagation(-50, 0.005, 50.0, tdse.Psi_start, write_spec, input);
	
		for (auto& psi : tdse.Psi_states) {
			tdse.outprojection(psi);
		}
	
		tdse.calc_spec(tdse.Psi, 0.0, 0.05, 50.0, 3, 0.0025, 1);
		cout << "Save data..." << endl;
		tdse.save_all();
	}


	/* ---------- Propagation Gauß gegen Potential ---------- */

	//double dx_au = 0.01;
	//
	//for (double E = 4; E <= 4; E += 0.5) {
	//	ostringstream oss1;
	//	oss1 << oss.str();
	//	oss1 << "E" << E << "eV/"; 
	//	cout << "E = " << E << " eV" << endl;
	//	path1 = oss1.str();
	//	fs::create_directory(path1);
	// 
	//	//double E = 15;
	//	double k = sqrt(2 * E / au_Ee);
	// 
	//	for (double dist_pot = 0.05; dist_pot <= 0.05; dist_pot += 2) {
	//		cout << "Potential:  " << dist_pot << endl;
	//		ostringstream oss2;
	//		oss2 << oss1.str();
	//		oss2 << "Dist" << dist_pot << "/";
	//	 
	//		TDSE tdse(dx_au);
	//		tdse.init_xlist_au(-300.0, dx_au, 500.0);
	//		tdse.init_tdse();
	//		tdse.path = oss2.str();
	//		fs::create_directory(tdse.path);
	//		 
	//		 
	//		tdse.add_fermibox(-dist_pot, 1, 300.0, 150.0);
	//		/*for (double xcoulomb = dist_pot; xcoulomb <= 800; xcoulomb += dist_pot) {
	//			tdse.add_Yukawa(xcoulomb, 2.0, 2.0);
	//		}*/
	//		tdse.set_Hamiltonian();
	//		 
	//		 
	//		double x_start = -20.0;
	//		tdse.Psi_start.init_gauss(tdse.xlist, x_start, 5.0, k);
	//		//tdse.Psi_start.init_gauss_time(tdse.xlist, x_start, 20.0, k, -200.0);
	//		 
	//		 
	//		double time_prop = 2 * abs(x_start) / k;
	//		time_prop = 200;
	//		cout << "Propagation time: " << time_prop << endl;
	//		tdse.propagation_au(0.0, 0.1, time_prop, tdse.Psi_start, write_spec, "0");
	//		 
	//		tdse.save_all();
	//	}
	//}
	return 0;
};
