clear all;
%%
%アルカリ原子超微細構造のゼーマンシフトの計算
%大阪市立大学レーザー量子物理学教室
%ver2 2021/8/30

%Wigner 3j-6jはVladimir Sovkov氏のコードを用いています。以下から関数一式をダウンロードしてください。
%https://jp.mathworks.com/matlabcentral/fileexchange/74069-wigner-3j-6j-9j
%Wigner3j.m, Wigner6j.m, ArranA.m, if3jc.m, sper.mが必要です。

%%
%入力

DataSave=0;%保存したい場合は1,しない場合は0

StartB=0;%開始磁場[Gauss]
FinalB=500;%終了磁場[Gauss]
StepB=1;%ステップ[Gauss]

%6Li,2S1/2状態パラメータ
%J=1/2;%電子の全角運動量
%I=1;%核スピン
%gJ=2.0023010;
%gI=-0.0004476540;
%Ahfs=152.1368407;%[MHz]

%6Li,2P1/2状態パラメータ
%J=1/2;%電子の全角運動量
%I=1;%核スピン
%gJ=0.6668;
%gI=-0.0004476540;
%Ahfs=17.375;%[MHz]

%6Li,2P3/2状態パラメータ
%J=3/2;%電子の全角運動量
%I=1;%核スピン
%gJ=1.335;
%gI=-0.0004476540;
%Ahfs=-1.155;%[MHz]
%Bhfs=-0.010;%[MHz]

%133CsD2状態パラメータ
J=3/2;%電子の全角運動量
I=7/2;%核スピン
gJ=1.334;
gI=-0.00039885395;
Ahfs=50.28827;%[MHz]
Bhfs=-0.4934;%[MHz]
Chfs=0.560*10^(-3);%[MHz]

uB=1.39962449361;%[MHz/Gauss]
B=StartB:StepB:FinalB;%磁場[Gauss]

%%
%全ての内部状態の準備

Nstate=(2*J+1)*(2*I+1);

F=abs(J-I):abs(J+I);
mF=-(J+I):(J+I);

Nfig=ceil(sqrt(length(mF)+1));%サブプロットのための数

F_list=[];
mF_list=[];

for k=1:length(F)
    F_list=[F_list,repmat(F(k),1,2*F(k)+1)];
    mF_list=[mF_list,-F(k):F(k)];
end
FmF=[F_list',mF_list'];%全ての|F,mF>の組み合わせを準備

[mJ,mI]=meshgrid(-J:J,-I:I);
mJmI=[mJ(:),mI(:)];%全ての|mJ,mI>の組み合わせを準備

%%
%関数の準備

Xk=@(F,k) (-1)^(I+J+F)*Wigner6j(F,J,I,k,I,J)/Wigner3j(J,k,J,-J,0,J)/Wigner3j(I,k,I,-I,0,I);

if ( I>1 && J>1 )
    H_hfs=@(F1,mF1,F2,mF2) (Ahfs*I*J*Xk(F1,1)+1/4*Bhfs*Xk(F1,2)+Chfs*Xk(F1,3))*(F1==F2)*(mF1==mF2);
    H_L=@(F,mF) Ahfs*I*J*Xk(F,1)+1/4*Bhfs*Xk(F,2)+Chfs*Xk(F,3);
elseif ( I>1/2 && J>1/2 )
    H_hfs=@(F1,mF1,F2,mF2) (Ahfs*I*J*Xk(F1,1)+1/4*Bhfs*Xk(F1,2))*(F1==F2)*(mF1==mF2);
    H_L=@(F,mF) Ahfs*I*J*Xk(F,1)+1/4*Bhfs*Xk(F,2);
else
    H_hfs=@(F1,mF1,F2,mF2) (Ahfs*I*J*Xk(F1,1))*(F1==F2)*(mF1==mF2);
    H_L=@(F,mF) Ahfs*I*J*Xk(F,1);
end

H_mag=@(F1,mF1,F2,mF2) uB*(-1)^(F1-mF1+J+I+1)*sqrt((2*F1+1)*(2*F2+1))*Wigner3j(F1,1,F2,-mF1,0,mF2)*(gJ*(-1)^(F2)*sqrt(J*(J+1)*(2*J+1))*Wigner6j(J,1,J,F2,I,F1)+gI*(-1)^(F1)*sqrt(I*(I+1)*(2*I+1))*Wigner6j(I,1,I,F2,J,F1))*(mF1==mF2);

H_H=@(mJ,mI,B) Ahfs*mJ*mI+uB*(gJ*mJ+gI*mI)*B;

%%
%計算箇所

E=[];
E_L=[];
E_H=[];
FmFmJmI=[];
fig1=figure;

for k=1:length(mF)
    
    %|F,mF>基底で任意の磁場での厳密計算
    subFmF=FmF(find(FmF(:,2)==mF(k)),:);
    NsubFmF=size(subFmF,1);
    
    M_hfs=zeros(NsubFmF);
    M_mag=zeros(NsubFmF);
    for i=1:NsubFmF
        for j=1:NsubFmF
            F1=subFmF(i,1);
            mF1=subFmF(i,2);
            F2=subFmF(j,1);
            mF2=subFmF(j,2);
            
            M_hfs(i,j)=H_hfs(F1,mF1,F2,mF2);
            M_mag(i,j)=H_mag(F1,mF1,F2,mF2);
            
        end
        
    end
    
    e=zeros(NsubFmF,length(B));
    for i=1:length(B)
        e(:,i) =sort(eig(M_hfs+B(i)*M_mag));
    end
    e=e';
    
    %低磁場の|F,mF>基底での近似計算
    e_L=zeros(NsubFmF,length(B));
    for i=1:NsubFmF
        F1=subFmF(i,1);
        mF1=subFmF(i,2);
        if F1>0%F=0ではない場合
            gF=@(F) gJ*(F.*(F+1)+J*(J+1)-I*(I+1))./(2*F.*(F+1))+gI*(F.*(F+1)-J*(J+1)+I*(I+1))./(2*F.*(F+1));
        else%F=0の場合
            gF=@(F) 0;
        end
        e_L(i,:) =H_L(F1,mF1)+uB*gF(F1)*mF1*B;
    end
    e_L=e_L';
    
    %高磁場の|mJ,mI>基底での近似計算
    submJmI=mJmI(find((mJmI(:,1)+mJmI(:,2))==mF(k)),:);
    NsubmJmI=size(submJmI,1);
    e_H=zeros(NsubmJmI,length(B));
    for i=1:NsubmJmI
        mJ1=submJmI(i,1);
        mI1=submJmI(i,2);
        e_H(i,:) =H_H(mJ1,mI1,B);
    end
    e_H=e_H';
    
    %低磁場の|F,mF>基底と高磁場の|mJ,mI>基底の対応付け
    [~,k_M]=sort(e(1,:));
    [~,k_L]=sort(e_L(1,:));
    [~,k_H]=sort(e_H(length(B),:));
    
    e=e(:,k_M);
    e_L=e_L(:,k_L);
    e_H=e_H(:,k_H);
    
    subFmFmJmI=[subFmF(k_L,:),submJmI(k_H,:)];
    
    %同じmF内での計算結果のプロット
    subplot(Nfig,Nfig,k)
    plot(B,e,'-');
    hold off
    xlabel('Magnetic field [Gauss]')
    ylabel('Energy [MHz]')
    title(strcat('m_F=',num2str(mF(k))))
    
    %計算結果を結合
    E=[E,e];
    E_L=[E_L,e_L];
    E_H=[E_H,e_H];
    FmFmJmI=[FmFmJmI;subFmFmJmI];
    
end

%%
%ゼロ磁場でエネルギーの低い順番に並べ替え

[~,k]=sort(E_L(2,:));
E=E(:,k);
E_L=E_L(:,k);
E_H=E_H(:,k);
FmFmJmI=FmFmJmI(k,:);

%%
%計算結果をプロット

%subplot(Nfig,Nfig,length(mF)+1)

w=round(length(B)/5);
fig2=figure;
subplot(1,3,[1,2])
plot(B,E,'k-');
hold on
plot(B(1:w),E_L(1:w,:),'g--');
plot(B(length(B)-w:length(B)),E_H(length(B)-w:length(B),:),'r--');
hold off;
xlim([StartB,FinalB])
ylim([min(E(:)),max(E(:))])
xlabel('Magnetic field [Gauss]')
ylabel('Energy [MHz]')
for k=1:Nstate
    text(FinalB,E(length(B),k),strcat('|F=',num2str(FmFmJmI(k,1)),', m_F=',num2str(FmFmJmI(k,2)),'> \rightarrow |m_J=',num2str(FmFmJmI(k,3)),', m_I=',num2str(FmFmJmI(k,4)),'>'));
end

%%
%計算結果を保存

out=[B',E];

if DataSave==1
    
    data=cell(size(out,1)+1,size(out,2));
    
    data(1,1)=cellstr('B[Gauss]');
    for k=1:Nstate
        data(1,1+k)=cellstr(strcat('|F=',num2str(FmFmJmI(k,1)),',mF=',num2str(FmFmJmI(k,2)),',mJ=',num2str(FmFmJmI(k,3)),',mI=',num2str(FmFmJmI(k,4)),'>'));
    end
    data(2:size(data,1),1:size(data,2))=num2cell(out);
    
    filter = {'*.csv'};
    [file, path] = uiputfile(filter);
    writecell(data,strcat(path,file));
    
end
