%%
%繰りこみヌメロフ法による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.
%スクリプト内コメントの式番号はこの論文の式番号
%
%製作者：堀越宗一
%製作日：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=2.5*g/mole;

V = @(x) 500*(x.^2-1).^2;%単位はcm^-1

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

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

xmin=-3.0;%単位はA0
xmax=3.0;%単位はA0
dx=0.0001;%単位はA0

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

Emin=0;%単位はcm^-1
Emax=1000;%単位はcm^-1
dE=0.005;%単位は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(2,2,[1,3])
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
hold off;
xlim([xmin xmax])
ylim([0 1000])
xlabel 'x[A0]'
ylabel 'Energy[cm^{-1}]'
title '計算結果'
grid on;

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に対する接合点のエラー値の変化（丸は決定された固有値の位置）'
