上海電機(jī)學(xué)院文理學(xué)院 上海 201306

摘 要：利用Matlab軟件的griddata函數(shù)對(duì)平均場方法獲取的原子核的勢能面進(jìn)行了插值平滑，通過繪制等高線，得到等勢能面，并以核素74Ge等勢能面的繪制為例，分析了74Ge的形變。

關(guān)鍵詞：原子核;等勢能面;Matlab

中圖分類號(hào)：O571.21

繪制原子核的等勢能面是分析核素形變的重要的手段。從系統(tǒng)的哈密頓量出發(fā)，利用受限的Hartree-Fock方法，[1]可以計(jì)算β自由度（描述原子核的軸對(duì)稱形變）和γ自由度（描述原子核的三軸形變）下的勢能。由于計(jì)算的核勢能是散點(diǎn)圖，還需通過數(shù)值方法對(duì)散點(diǎn)進(jìn)行插值計(jì)算，得到光滑的勢能面，并在β-γ極坐標(biāo)平面上通過繪制等高線得到等勢能面，從而進(jìn)一步判斷原子核的形變情況。

Matlab作為一個(gè)具有較強(qiáng)功能的科學(xué)計(jì)算軟件，能夠直接調(diào)用數(shù)值方法的計(jì)算程序，可以快速高效地實(shí)現(xiàn)各種數(shù)值計(jì)算過程。本文將利用Matlab軟件中的griddata函數(shù)對(duì)離散的原子核勢能進(jìn)行插值平滑，具體說明Matlab軟件在原子核等勢能面繪制中的應(yīng)用。

1 griddata函數(shù)簡介

Matlab的griddata函數(shù)可以對(duì)散點(diǎn)數(shù)據(jù)進(jìn)行兩維和三維插值。以兩維插值為例，其命令語法如下：

zq= griddata（x，y，z，xq，yq，method）

上述命令可以實(shí)現(xiàn)將（x，y，z）的散點(diǎn)數(shù)據(jù)平滑成z=f（x，y）形式的曲面，在（xq，yq）指定的查詢點(diǎn)進(jìn)行插值，并返回插入的值zq，其中xq和yq為兩維查詢點(diǎn)網(wǎng)格，zq則為兩維數(shù)組。method為插值方法，可以選擇linear、nearest、natural、cubic、linear 或 v4，其中前四個(gè)選項(xiàng)是基于三角剖分的插值方法，而 v4則為雙調(diào)和樣條插值方法，默認(rèn)方法為linear。本文中，我們采用cubic選項(xiàng)，即基于三角剖分的三次插值方法。griddata函數(shù)的詳細(xì)語法說明可見Matlab軟件的幫助文件。

下面我們利用griddata函數(shù)來繪制核素74Ge的等勢能面。

2 核素74Ge等勢能面繪制

在原子核的殼模型中，以核素56Ni為不活躍的殼芯，74Ge的價(jià)核子可以在pf5/2g9/2的模型空間中被激發(fā)。利用受限的Hartree-Fock方法，通過JUN45有效相互作用，我們計(jì)算了β和γ自由度下，體系總角動(dòng)量為0時(shí)，74Ge的勢能。相關(guān)的計(jì)算方法參見文獻(xiàn)[2]的討論，這里不再贅述。計(jì)算得到的勢能分布可見圖（a）中的散點(diǎn)。利用griddata函數(shù)可對(duì)散點(diǎn)進(jìn)行插值平滑，再進(jìn)一步地將平滑后的勢能面投影到β-γ極坐標(biāo)平面上，就可以做出等勢能面圖。具體的計(jì)算過程見下文的Matlab程序代碼和對(duì)應(yīng)的注釋：

r=load（'fort.30'）;%讀取計(jì)算的離散的核勢能

subplot（2，1，1）;

plot3（r（：，1），r（：，2），r（：，3），'r.'）;%繪制勢能的三維散點(diǎn)圖

xlabel（'＼＼it x'）; % x=β

ylabel（'＼＼it y'）; % y=β *sin γ

zlabel（'{＼＼it E}（MeV）'）; %E為勢能

title（'（a）'）; hold on;

xlin=linspace（min（r（：，1）），max（r（：，1）），800）;

ylin=linspace（min（r（：，2）），max（r（：，2）），800）;

[X，Y]=meshgrid（xlin，ylin）; %將插值計(jì)算的區(qū)域網(wǎng)格化

Z=griddata（r（：，1），r（：，2），r（：，3），X，Y，'cubic'）; %利用griddata函數(shù)插值平滑離散點(diǎn)，并在網(wǎng)格化的查詢點(diǎn)[X，Y]上，返回插入的值Z

mesh（X，Y，Z）; %繪制平滑后的勢能面

grid on;box on;

subplot（2，1，2）;

contourf（X，Y，Z，12）; %將勢能面投影到β-γ極坐標(biāo)平面上，繪制等高線，得到等勢能面

title（'（b）'）;h=colorbar;

set（get（h，'Title'），'string'，'{＼＼it E}（MeV）'）;hold on;

sector（0，0，0，60，0.25）; % 利用自定義的函數(shù)sector 繪制β-γ平面極坐標(biāo)

function sector（x0，y0，angle1，angle2，R） %定義sector函數(shù)

t=linspace（angle1/180*pi，angle2/180*pi，1000）;

x=R*cos（t）;y=R*sin（t）;x（1001）=0;y（1001）=0;

plot（x，y，'k-'）;axis（[x0，x0+R，y0，y0+R]）;

set（gca，'YColor'，'white'）;box off;

x1=R*cos（30/180*pi）;y1=0.01+R*sin（30/180*pi）;

x2=R*cos（60/180*pi）;y2=0.01+R*sin（60/180*pi）;

xlabel（{'＼＼it ＼＼beta'}）;text（x1，y1，'{＼＼it ＼＼gamma}=30^{＼＼circ}'）;

text（x2，y2，'60^{＼＼circ}'）;text（R，0，'0^{＼＼circ}'）;

end

（a）插值平滑后的勢能面

（b）對(duì)應(yīng)的等勢能面

執(zhí)行上述代碼后，得到的計(jì)算結(jié)果如圖所示。在圖（a）中可以看到，通過griddata函數(shù)，利用基于三角剖分的三次插值方法能夠很好地平滑離散點(diǎn)，得到光滑連續(xù)的勢能面。圖（b）中，在β-γ極坐標(biāo)平面上所繪制的等勢能面也能夠準(zhǔn)確地描述核素74Ge的形變。當(dāng)γ處于450到550的范圍內(nèi)時(shí)，存在勢能的極小值，表明74Ge具有一定的三軸形變，這與實(shí)驗(yàn)上的分析是相一致的。[3]

3 結(jié)語

利用Matlab軟件的griddata函數(shù)對(duì)受限的Hartree-Fock方法計(jì)算得到的核素74Ge的離散的勢能面進(jìn)行了基于三角剖分的三次插值平滑，繪制的等勢能面能夠很好地反映出74Ge的結(jié)構(gòu)特征。由于Matlab語言編寫簡單，函數(shù)功能豐富，本文提供的方法不失為繪制原子核等勢能面的一種有效途徑。

參考文獻(xiàn)：

[1]RING P，SCHUCK P.The Nuclear Many-Body Problem[M].New York：Spring-Verlag，1980：266.

[2]金華.76Ge 形變的殼模型分析.上海電機(jī)學(xué)院學(xué)報(bào)[J].2017，20（5）：301-305.

[3]SUN J J，SHI Z，LI X Q，et al.Spectroscopy of 74Ge：From soft to rigid triaxiality[J].Physics Letters B，2014，734：308-313.

作者簡介：金華（1979-），男，上海電機(jī)學(xué)院教師，主要從事物理學(xué)的教學(xué)與研究工作。