#include"Potential.h"
#include<iostream>
#include<math.h>
#include<fstream>
#include<time.h>
#include<stdio.h>
#include<stdlib.h>
#include<complex>
#include <cmath>


Potential::Potential(double& dx_in) :dx(dx_in) {};
Potential::Potential(const int n, double& dx_in) :dx(dx_in) {
	this->init(n);
}

void Potential::init_box(Cvec xlist, double V0, double length, double x0) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	for (int i = 0; i < size; i++) {
		if (real(xlist[i]) >= x0 - length / 2 && real(xlist[i]) <= x0 + length / 2) {
			vec[i] += V0;
		}
		else {
			vec[i] += 0;
		}
	}	
}
void Potential::init_fermibox(Cvec xlist, double depth, double s, double length, double x0) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	for (int i = 0; i < size; i++) {
		if (xlist[i].real() < x0) vec[i] += depth / (exp(-(xlist[i] - (x0 - length / 2)) / s) + 1.0);
		else vec[i] += depth / (exp((xlist[i] - (x0 + length/2)) / s) + 1.0);

	}
}
void Potential::init_Coulomb(Cvec xlist, double x0, int Z, double a) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	for (int i = 0; i < size; i++) vec[i] += -1.0*Z / sqrt(pow((real(xlist[i]) - x0), 2) + a);
}
void Potential::init_Yukawa(Cvec xlist, double x0, int Z, double a) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	for (int i = 0; i < size; i++) vec[i] += -1.0 * Z / sqrt(pow((real(xlist[i]) - x0), 2) + a) * exp(-abs(xlist[i] - x0) / 0.1);
}
void Potential::init_step(Cvec xlist, double x0, double height) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	for (int i = 0; i < size; i++) {
		if (real(xlist[i]) >= x0) {
			vec[i] = height;
		}
		else {
			vec[i] = 0;
		}
	}
}


//void Potential::init_imag(Cvec xlist, double x0, Creal value, double power) {
//	if (xlist.get_size() != this->size) {
//		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
//		exit(0);
//	}
//	int x_size = xlist.get_size();
//	double xleft = xlist[0].real();
//	double xright = xlist[x_size - 1].real();
//	for (int i = 0; i < x_size; i++) {
//		double x_v = xlist[i].real();
//		if (x_v > xright - x0) {
//			vec[i] = -value * pow((x_v - (xright - x0)), power);
//		}
//		else if (x_v < xleft + x0) {
//			vec[i] = -value * pow((x_v - (xleft + x0)), power);
//		}
//	}
//}
void Potential::init_imag(Cvec xlist, double x0, Creal value, double power) {
	if (xlist.get_size() != this->size) {
		std::cout << "ERROR, length of xlist does not fit length of potential" << std::endl;
		exit(0);
	}
	int x_size = xlist.get_size();
	double xleft = xlist[0].real();
	double xright = xlist[x_size - 1].real();
	double x_now = xleft;
	for (int i = 0; i < x_size; i++) {
		x_now = xlist[i].real();
		vec[i] = -value / pow(cosh((xright - x_now) / x0), power);
		vec[i] += -value / pow(cosh((x_now - xleft) / x0), power);
	}
}

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

	Potential out(size, dx);
	for (int i = 0; i < size; i++) out.vec[i] = vec[i] + in.vec[i];
	return out;
}
Potential Potential::operator+(const Cvec& in) {
	if (this->size != in.size) {
		std::cout << "ERROR, Vectors must have the same length" << std::endl;
		exit(0);
	}

	Potential out(size, dx);
	for (int i = 0; i < size; i++) out.vec[i] = vec[i] + in.vec[i];
	return out;
}
Potential Potential::operator*(const double& in) {
	Potential out(size, dx);
	for (int i = 0; i < size; i++) out[i] = vec[i] * in;
	return out;
}


double Potential::minimum() {
	double res = vec[0].real();
	for (int i = 1; i < size; i++) {
		if (vec[i].real() < res) res = vec[i].real();
	}
	return res;
}