////////////////////////////////////////////////////////////////////
// ステップ関数型ポテンシャルへの衝突
////////////////////////////////////////////////////////////////////
#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 = 5000, Ld = 10;
//空間分割サイズ
double dz = 1.0E-11;
//重ね合わせの数
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;
//電子波のエネルギー
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 = E0;

//計算時間の幅
const int ts = -150, te = 300;
//時間間隔
double dt = 1.0 * 1.0E-16;
/////////////////////////////////////////////////////////////////
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: true false false false" << std::endl;
	fout_1D << "#xrange:" << -Lz / 2 << " " << Lz / 2 << " " << Lz / 10 << std::endl;
	fout_1D << "#yrange:" << -0.30 << " " << 0.30 << " " << 0.05 << std::endl;
	fout_1D << "#all:" << std::endl;
	fout_1D << "0 -0.5" << std::endl;
	fout_1D << "0 0.0" << std::endl;
	fout_1D << "2500 0.0" << 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);

			for (int jz = 0; jz <= N; jz++) {
				double k = (k0 + dk * double(jz - N / 2));
				double omega = hbar / (2.0*me) * pow(k, 2);
				double E = hbar * omega;

				std::complex<double> k1 = std::complex<double>( sqrt(2.0*me*(E - V1)) / hbar, 0.0);
				std::complex<double> k2;
				if (E > V2) {
					k2 = sqrt(2.0*me*(E - V2)) / hbar;
				}else {
					k2 = I * sqrt(2.0*me*(V2 - E)) / hbar;
				}
				std::complex<double> r, t;
				r = (k1 - k2) / (k1 + k2);
				t = (2.0 * k1) / (k1 + k2);

				if (z < 0) {

					Psi += (exp(I * k1 * z)+r* exp(-I * k1 * z))* exp(-I * omega * t_real) * exp(-1.0 / 2.0 * pow((k - k0) / sigma, 2));

				}else {

					Psi += (t * exp(I * k2 * z))* exp(-I * omega * t_real) * exp(-1.0 / 2.0 * pow((k - k0) / sigma, 2));

				}
			}

			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
*/