////////////////////////////////////////////////////////////////////
// 箱型ポテンシャル障壁によるトンネル効果
////////////////////////////////////////////////////////////////////
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cstdio>
#include <iomanip>
#include <stdio.h>
#include <complex>
#include <direct.h>

/////////////////////////////////////////////////////////////////
//物理定数
/////////////////////////////////////////////////////////////////
//光速
const double c = 2.99792458E+8;
//真空の透磁率
const double mu0 = 4.0*M_PI*1.0E-7;
//真空の誘電率
const double epsilon0 = 1.0 / (4.0*M_PI*c*c)*1.0E+7;
//プランク定数
const double h = 6.6260896 * 1.0E-34;
double hbar = h / (2.0*M_PI);
//電子の質量
const double me = 9.10938215 * 1.0E-31;
//電子ボルト
const double eV = 1.60217733 * 1.0E-19;
//複素数
const std::complex<double> I = std::complex<double>(0.0, 1.0);
/////////////////////////////////////////////////////////////////
//物理系の設定
/////////////////////////////////////////////////////////////////
//空間の大きさ
const int Lz = 10000, Ld = 20;
//空間分割サイズ
double dz = 1.0E-13;
//重ね合わせの数
const int N = 5000;
//パルスの幅
double delta_z = 0.5E-8;
double sigma = 2.0*sqrt(2.0*log(2.0)) / (delta_z);
//波数の間隔
double dk = 30.0 / (delta_z * double(N + 1));
//z方向の原点の場所
const int iz0 = Lz / 2 - Lz/10;
//電子波のエネルギー
const double E0 = 10.0 * eV;
//波数の中心
double k0 = sqrt(2.0 * me * E0 / pow(hbar, 2));
//角振動数の中心
double omega0 = hbar / (2.0*me) * pow(k0, 2);

//ポテンシャル
double V1 = 0;
double V2 = 20.0 * eV;
double V3 = 0;

//ポテンシャル障壁の幅
double d = 2000 * dz;

//計算時間の幅
const int ts = 0, te = 100;
//時間間隔
double dt = 2.0*M_PI/omega0 / (te - ts + 2); //1.0 * 1.0E-16;


/////////////////////////////////////////////////////
//Matrix22Complex クラスの定義
/////////////////////////////////////////////////////
class Matrix22Complex {
private:
	std::complex<double> a[2][2];
public:
	void set(int i, int j, std::complex<double> k){
		a[i][j] = k;
	}
	std::complex<double> get(int i, int j){
		return a[i][j];
	}
	Matrix22Complex operator*(Matrix22Complex& c){
		Matrix22Complex temp;
		int i, j, k;
		for (i = 0; i<2; i++) {
			for (j = 0; j<2; j++) {
				temp.a[i][j] = 0;
				for (k = 0; k<2; k++) {
					temp.a[i][j] += a[i][k] * c.a[k][j];
				}
			}
		}
		return temp;
	}
};
/////////////////////////////////////////////////////
//反射係数と透過係数を計算する関数定義
/////////////////////////////////////////////////////
std::complex<double> ReflectionCoefficient(Matrix22Complex& MM) {
	return - ( MM.get(1, 0) / MM.get(1, 1) );
}
std::complex<double> TransmissionCoefficient(Matrix22Complex& MM) {
	return ( MM.get(0, 0) * MM.get(1, 1) - MM.get(0, 1) * MM.get(1, 0) ) / MM.get(1, 1);
}

/////////////////////////////////////////////////////////////////
std::string folder = "output";//作成するフォルダ名
std::string filename = "pulse1.txt";//作成するフォルダ名
const int precision_N = 4;

int main() {

	//フォルダ生成
	_mkdir(folder.c_str());
	std::string filepath = folder + "/" + filename;

	//出力ストリームによるファイルオープン
	std::ofstream fout_1D;
	fout_1D.open(filepath.c_str());

	fout_1D << "#x:位置" << std::endl;
	fout_1D << "#y:確率振幅" << std::endl;
	fout_1D << "#legend: ポテンシャル障壁 実部 虚部 絶対値" << std::endl;
	fout_1D << "#showLines: true true true true" << std::endl;
	fout_1D << "#showMarkers: false false false false" << std::endl;
	fout_1D << "#fills: false false false false" << std::endl;
	fout_1D << "#xrange:" << -iz0 << " " << Lz-iz0 << " " << Lz / 10 << std::endl;
	fout_1D << "#yrange:" << -2.50 << " " << 2.5 << " " << 0.5 << std::endl;

	fout_1D << "#all:" << std::endl;
	fout_1D << -iz0 <<" -2.3" << std::endl;
	fout_1D << "0 -2.3" << std::endl;
	fout_1D << "0 2.3" << std::endl;
	fout_1D << d /dz << " 2.3" << std::endl;
	fout_1D << d / dz << " -2.3" << std::endl;
	fout_1D << Lz - iz0 << " -2.3" << std::endl;

	//各時刻における計算を行う
	for (int tn = ts; tn <= te; tn++) {
		double t_real = dt * double(tn);
		std::cout << tn << std::endl;
		fout_1D << "#coma:" << tn << std::endl;

		for (int iz = 0; iz <= Lz; iz = iz + Ld) {
			double z = dz * double(iz - iz0);
			std::complex<double> Psi = std::complex<double>(0.0, 0.0);

			double k = k0;
			double omega = hbar / ( 2.0 * me ) * pow(k, 2);
			double E = hbar * omega;

			std::complex<double> k1 = sqrt(std::complex<double>(pow(k, 2) - (2.0 * me * V1) / pow(hbar, 2), 0.0));
			std::complex<double> k2 = sqrt(std::complex<double>(pow(k, 2) - (2.0 * me * V2) / pow(hbar, 2), 0.0));
			std::complex<double> k3 = sqrt(std::complex<double>(pow(k, 2) - (2.0 * me * V3) / pow(hbar, 2), 0.0));

			Matrix22Complex M21, T2, M32, M31;
			M21.set(0, 0, (std::complex<double>(1.0, 0.0) + k1 / k2) / 2.0);
			M21.set(0, 1, (std::complex<double>(1.0, 0.0) - k1 / k2) / 2.0);
			M21.set(1, 0, (std::complex<double>(1.0, 0.0) - k1 / k2) / 2.0);
			M21.set(1, 1, (std::complex<double>(1.0, 0.0) + k1 / k2) / 2.0);
			M32.set(0, 0, (std::complex<double>(1.0, 0.0) + k2 / k3) / 2.0);
			M32.set(0, 1, (std::complex<double>(1.0, 0.0) - k2 / k3) / 2.0);
			M32.set(1, 0, (std::complex<double>(1.0, 0.0) - k2 / k3) / 2.0);
			M32.set(1, 1, (std::complex<double>(1.0, 0.0) + k2 / k3) / 2.0);
			T2.set(0, 0, exp(I * k2 * d));
			T2.set(0, 1, 0);
			T2.set(1, 0, 0);
			T2.set(1, 1, exp(-I * k2 * d));

			M31 = M32 * T2 * M21;

			std::complex<double> r = ReflectionCoefficient(M31);
			std::complex<double> t = TransmissionCoefficient(M31);

			if (z < 0) {

				Psi += (exp(I * k1 * z) + r * exp(-I * k1 * z))* exp(-I * omega * t_real);

			}else if ( z < d ) {

				std::complex<double> Ap = M21.get(0, 0) * 1.0 + M21.get(0, 1) * r;
				std::complex<double> Am = M21.get(1, 0) * 1.0 + M21.get(1, 1) * r;

				Psi += ( Ap * exp(I * k2 * z) + Am * exp(-I * k2 * z))* exp(-I * omega * t_real);

			}else {

				Psi += (t * exp(I * k3 * (z - d)))* exp(-I * omega * t_real);

			}

//			Psi = Psi / double(N);
			fout_1D << std::setprecision(precision_N);
			fout_1D << z / dz << " " << Psi.real() << " " << Psi.imag() << " " << abs(Psi) << std::endl;
		}
		fout_1D << std::endl;
	}

	fout_1D.close();
}

/*
gccコンパイル
g++ main.cpp
*/