#include"Wavefunction.h"
#include<iostream>
#define _USE_MATH_DEFINES
#include<math.h>
#include<fstream>
#include<time.h>
#include<stdio.h>
#include<stdlib.h>
#include<complex>
#include <cmath>

using namespace std;
Wavefunction::Wavefunction(double& dx_in) :dx(dx_in) {};
Wavefunction::Wavefunction(const int n, double& dx_in) :dx(dx_in) {
	this->init(n);
}

Creal Wavefunction::scalar(Wavefunction& rhs) {
	if (size != rhs.size) {
		cout << "ERROR, Vectors must have the same length" << endl;
		exit(0);
	}
	Creal out(0, 0);
	for (int i = 0; i < size; i++) {
		out += conj(vec[i]) * rhs[i];
	}
	out *= dx;
	return out;
}

double Wavefunction::norm() {
	return scalar(*this).real();
}
void Wavefunction::normalize() {
	double n = norm();
	*this /= sqrt(n);
}

double& Wavefunction::get_dx() const{
	return dx;
}

void Wavefunction::set_energy(double Energy) {
	E0 = Energy;
}

double Wavefunction::get_energy() const {
	return E0;
}

Wavefunction& Wavefunction::operator=(const double& in) {
	for (int i = 0; i < size; i++) this->vec[i] = Creal(in, 0.0);
	return *this;
}
Wavefunction& Wavefunction::operator=(const Wavefunction& in) {
	if (this->size != in.size) this->init(in.size);
	for (int i = 0; i < size; i++) this->vec[i] = in.vec[i];
	this->E0 = in.E0;
	return *this;
}
Wavefunction Wavefunction::operator-(const Wavefunction& in) {
	if (size != in.size) {
		std::cout << "ERROR, Vectors must have the same length" << std::endl;
		exit(0);
	}
	Wavefunction out(size,dx);
	for (int i = 0; i < size; i++) out[i] = (*this)[i] - in[i];
	return out;
}
Wavefunction Wavefunction::operator/(const double& in) {
	if (in == 0) {
		std::cout << "ERROR: Divide by Zero" << std::endl;
		exit(0);
	}
	Wavefunction out(size, dx);
	for (int i = 0; i < size; i++) out[i] = vec[i] / in;
	return out;
}
Wavefunction Wavefunction::operator*(const double& in) {
	Wavefunction out(size,dx);
	for (int i = 0; i < size; i++) out[i] = vec[i] * in;
	return out;
}
Wavefunction Wavefunction::operator*(const Creal& in) {
	Wavefunction out(size,dx);
	for (int i = 0; i < size; i++) out[i] = vec[i] * in;
	return out;
}
Wavefunction Wavefunction::operator*(const Cvec& in) {
	if (size != in.size) {
		std::cout << "ERROR, Vectors must have the same length" << std::endl;
		exit(0);
	}
	Wavefunction out(size, dx);
	for (int i = 0; i < size; i++) out[i] = (*this)[i] * in[i];
	return out;
}

bool Wavefunction::operator<(const Wavefunction& in) const {
	return (get_energy() < in.get_energy());
}
bool Wavefunction::operator>(const Wavefunction& in) const {
	return (get_energy() > in.get_energy());
}

void Wavefunction::init_random() {
	srand(time(NULL));
	for (int i = 0; i < size; i++) vec[i] = rand();
}
void Wavefunction::init_gauss(Cvec xlist, double x0, double sigma, double k) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of wavefunction" << std::endl;
			exit(0);
	}
	double C = 1 / sqrt(M_PI * powf(sigma, 2));
	Creal i(0, 1);
	for (int j = 0; j < size; j++) vec[j] = sqrt(C) * exp(-(xlist[j] - x0) * (xlist[j] - x0) / (2 * sigma * sigma)) * exp(i * k * xlist[j]);
}

void Wavefunction::init_gauss_time(Cvec xlist, double x0, double sigma, double k, double t) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of wavefunction" << std::endl;
		exit(0);
	}
	
	Creal i(0, 1);
	for (int j = 0; j < size; j++) {
		vec[j] = exp(Creal(0.0, k * (xlist[j].real() - x0))) * exp(-Creal(0, pow(k, 2) * t / 2.0)) * exp(-pow(xlist[j] - x0 - k * t, 2.0) / (2.0 * pow(sigma, 2.0) * (1.0 + Creal(0, t / pow(sigma, 2.0)))));
	}
	this->normalize();
}

void Wavefunction::init_planeWave(Cvec xlist, Creal Ampl, double k, double phase) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of wavefunction" << std::endl;
		exit(0);
	}
	Creal i(0, 1);
	for (int j = 0; j < size; j++) vec[j] = Ampl * exp(-i * (k * xlist[j] + phase));
}