%%
%繰りこみヌメロフ法による1次元時間に依存しないシュレーディンガー方程式の数値解法
%
%モース型ポテンシャル
%
%参考：Johnson, Bernard R. "New numerical methods applied to solving the one‐dimensional eigenvalue problem."
%      The Journal of Chemical Physics 67.9 (1977): 4086-4093.
%スクリプト内コメントの式番号はこの論文の式番号
%
%固有値の正解値は、
%Wicke, Brian G., and David O. Harris. "Comparison of three numerical techniques for calculating eigenvalues of an unsymmetrical double minimum oscillator."
%The Journal of Chemical Physics 64.12 (1976): 5236-5242.
%のテーブルⅡの値を使用
%
%製作者：堀越宗一
%製作日：2019/08/09
%最終更新日：2019/08/13

clear all;

%%
%ポテンシャルの準備
%初めからcm^-1とA^0の単位系で計算(光速とプランク定数をcmの単位にしておく)

cm=10^(-2);
K=cm^(-1);
A0=10^(-10);
c=299792458*10^(2);%cm
h=6.62607004*10^(-34)*10^4;%cm
hbar=h/(2*pi);%cm
mole=6.02214076*10^23;
g=1*10^(-3);

mu=5*g/mole;
omegae=1000;
omwgae_xe=8;
B=sqrt(8*pi^2*c*mu*omwgae_xe/h)*10^(-8);
D=31250;
xa=1.5;

V = @(x) D*(1-exp(-B*(x-xa))).^2;%単位はcm^-1

%%
%n=0～15までの固有値の正解値、

En=[498.00000,1482.00000,2450.00000,3402.00000,4338.00000,5258.00000,6162.00000,7050.00000,...
    7922.00000,8778.00000,9618.00000,10442.00000,11250.00000,12042.00000,12818.00000,13578.00000];

%%
%ここからヌメロフ法

%%
%グリッドの設定（ここでは等間隔のグリッド）

xmin=0.0;%単位はA0
xmax=15.0;%単位はA0
dx=0.0001;%単位はA0

x=xmin:dx:xmax;%単位はA0
Numx=length(x);

Emin=0;%単位はcm^-1
Emax=31250;%単位はcm^-1
dE=1;%単位はcm^-1
E=Emin:dE:Emax;%単位はcm^-1
NumE=length(E);

NumNode=zeros(1,NumE);
DE=zeros(1,NumE);

%解きたい時間に依存しないシュレーディンガー方程式：(2)式
Qn=@(E) (2*mu/hbar^2)*(E-V(x))*(h*c)*K^2*A0^2;%cm^-1からA0に単位を変換

parfor i=1:NumE
    
    %繰りこみヌメロフ法で波動関数、接合点、節の数、接合エラー値を計算
    [DE(i),NumNode(i),Psin] = RNNumerov(Qn(E(i)),dx);%xはA0単位の座標、QnもA0単位で与える
    
end

%%
%固有値の高精度計算（各節の数においてDE(E)=0となるエネルギーを決める）

n=0:max(NumNode)-1;
Numn=length(n);
Eigenvalues=zeros(1,Numn);

parfor i=1:Numn
    
    k=find(NumNode==n(i));
    
    errf = @(E) (RNNumerov(Qn(E),dx)).^2;
    Eigenvalues(i)=fminbnd(errf,min(E(k)),max(E(k)));
    
end

%%
%プロット

subplot(1,2,1)
plot(x,V(x),'k-','LineWidth',1)
hold on;
for i=1:Numn
    [~,~,Psin] = RNNumerov(Qn(Eigenvalues(i)),dx);
    Nc=trapz(Psin.^2);
    plot([xmin xmax],[Eigenvalues(i) Eigenvalues(i)],'k--',x,2.5*10^4*Psin/sqrt(Nc)+Eigenvalues(i),'-')
end
for i=1:16
    plot([3 3.5],[En(i) En(i)],'k-','LineWidth',1)
end
hold off;
xlim([xmin xmax])
ylim([0 35000])
xlabel 'x[A0]'
ylabel 'Energy[cm^{-1}]'
title '計算結果（太い黒線は正値）'

subplot(1,2,2)
plot(x,V(x),'k-','LineWidth',1)
hold on;
for i=1:Numn
    [~,~,Psin] = RNNumerov(Qn(Eigenvalues(i)),dx);
    plot([xmin xmax],[Eigenvalues(i) Eigenvalues(i)],'k--',x,5*10^3*Psin+Eigenvalues(i),'-')
end
for i=1:16
    plot([3 3.5],[En(i) En(i)],'k-','LineWidth',1)
end
hold off;
xlim([xmin xmax])
ylim([31000 31300])
xlabel 'x[A0]'
ylabel 'Energy[cm^{-1}]'
title '計算結果（乖離極限）'

%{

subplot(2,2,2)
plot(E,NumNode,'k.-')
hold on;
plot(Eigenvalues,n,'ro')
hold off;
xlabel 'E[cm^{-1}]'
ylabel 'Count of nodes'
xlim([min(E) max(E)])
title '試行固有値Eに対する節の数の変化（丸は決定された固有値の位置）'

subplot(2,2,4)
plot(E,DE,'k.-')
hold on;
plot(Eigenvalues,zeros(1,length(Eigenvalues)),'ro')
plot([min(E) max(E)],[0 0],'k--')
hold off;
xlabel 'E[cm^{-1}]'
ylabel 'DE(E)'
xlim([min(E) max(E)])
ylim([-0.01 0.01])
title '試行固有値Eに対する接合点のエラー値の変化（丸は決定された固有値の位置）'

%}