MC方法模拟氢原子基态电子云（二维）

close all
clear
clc

f=figure();
% f.Position=[0 0 1300 1080];
num=1e+5;
a1=0.0529;    
r=[0.25 1 2];    
d=[1.1 0.14 0.2]; 

for i=1:3
    r0=r(i);
    Dmax=d(i);

    cp=Dmax*rand(1,num);       %     cp choosed point 
    rp=r0*rand(1,num);              %     rp right point
    D=[4*rp.^2/a1^3.*exp(-2*rp/a1) ;
           rp.^2/(8*a1^3).*(2-rp/a1).^2.*exp(-rp/a1);
           4/3/a1^3*(rp/81).^2.*(27-18*rp/a1+2*(rp/a1).^2).^2.*exp(-2*rp/3/a1)];
                                                  %     s1,s2,s3 概率公式
    p=find(cp<=D(i,:));
    tp=rp(p);
    s=length(p);
    sita=2*pi*rand(1,s);
    x=tp.*cos(sita);
    y=tp.*sin(sita);

    text=strcat(num2str(i),'S态氢原子电子云');
    text1=strcat(text,'分布');
    text2=strcat(text,'概率曲线');

    subplot(2,3,i)
    plot(x,y,'.');
    axis([0 r0 0 r0]);
    title(text1);

    subplot(2,3,3+i)
    plot(rp,D(i,:),'.');
    axis([0 r0 0 11]);
    title(text2);
end