Contents

球内の乱数

Matlabのドキュメントによると球内の乱数は以下のコードで作成されます。

clear;close all;
rng(0,'twister');
N=1000;
rvals=2*rand(N,1)-1;
elevation=asin(rvals);
azimuth=2*pi*rand(N,1);
radii=rand(N,1).^(1/3);
[x,y,z]=sph2cart(azimuth,elevation,radii);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;

5層の球内の乱数

上記のコードを発展させます。

clear;close all;
rng(0,'twister');
N=1000;
map=jet(5);
for layer=1:5
    rvals=2*rand(N,1)-1;
    elevation=asin(rvals);
    azimuth=2*pi*rand(N,1);
    radii=rand(N,1).^(1/3)+layer;
    [x,y,z]=sph2cart(azimuth,elevation,radii);
    if layer==1;figure('color',[1,1,1]);end;
    h=plot3(x,y,z,'b.');hold on;
    set(h,'color',map(layer,:));
end;
axis equal;axis tight;
rotate3d on;

5層の球内の乱数の内部表示

手前側の1/8を削ってみました。

clear;close all;
rng(0,'twister');
N=1000;
map=jet(5);
for layer=1:5
    rvals=2*rand(N,1)-1;
    elevation=asin(rvals);
    azimuth=2*pi*rand(N,1);
    t=find(elevation<=0 | (azimuth<=pi | azimuth>=1.5*pi));
    radii=rand(length(t),1).^(1/3)+layer;
    [x,y,z]=sph2cart(azimuth(t),elevation(t),radii);
    if layer==1;figure('color',[1,1,1]);end;
    h=plot3(x,y,z,'b.');hold on;
    set(h,'color',map(layer,:));
end;
axis equal;axis tight;
rotate3d on;
axis equal;axis tight;

5層の球内の乱数の各層を正規分布として内部表示

rand関数の代わりにrandn関数を使います。 マイナスの立方根を計算しないで済むように幅±1を6σとします。

clear;close all;
rng(0,'twister');
N=1000;
map=jet(5);
for layer=1:5
    rvals=2*rand(N,1)-1;
    elevation=asin(rvals);
    azimuth=2*pi*rand(N,1);
    t=find(elevation<=0 | (azimuth<=pi | azimuth>=1.5*pi));
    radii=(randn(length(t),1)/6+1).^(1/3)+layer;
    [x,y,z]=sph2cart(azimuth(t),elevation(t),radii);
    if layer==1;figure('color',[1,1,1]);end;
    h=plot3(x,y,z,'b.');hold on;
    set(h,'color',map(layer,:));
end;
axis equal;axis tight;
rotate3d on;
axis equal;axis tight;

球面調和関数の結果を乱数で表示

von Neumannの棄却法という方法を使いました。 中心からの距離は半径1+正規分布/6σとしました。

s軌道

clear;close all;
rng(0,'twister');
N=1000;rd=1;
theta=2*rand(N,1)-1;
theta=asin(theta);
phi=2*pi*rand(N,1);
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^0_0) \ \ s \ \ 1$$','interpreter','latex','FontSize',18);

p軌道

Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
t=rand(Nn,1);
th=theta+pi/2;%補正用
k=cos(th);
t=find(t<abs(k));
theta=theta(t(1:N));
phi=2*pi*rand(N,1);
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^0_1) \ \ p_z \ \ \cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(phi).*sin(th);
t=find(t<abs(k));
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^1_1) \ \ p_x \ \ \cos\phi\sin\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=sin(phi).*sin(th);
t=find(t<abs(k));
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{1}_1) \ \ p_y \ \ \sin\phi\sin\theta$$','interpreter','latex','FontSize',18);

d軌道

Nn=N*4;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(th).^2;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{0}_2) \ \ d_{3z^2-r^2} \ \ 3\cos^2\theta-1$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=2*cos(phi).*sin(th).*cos(th);
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{1}_2) \ \ d_{xz} \ \ \cos\phi\sin\theta\cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=2*sin(phi).*sin(th).*cos(th);
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{1}_2) \ \ d_{yz} \ \ \sin\phi\sin\theta\cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(2*phi).*sin(th).^2;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{2}_2) \ \ d_{x^2-y^2} \ \ \cos2\phi\sin^2\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=sin(2*phi).*sin(th).^2;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{2}_2) \ \ d_{xy} \ \ \sin2\phi\sin^2\theta$$','interpreter','latex','FontSize',18);

f軌道

Nn=N*4;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=(-3*cos(th)+5*cos(th).^3)/2;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{0}_3) \ \ f_{5z^3-3zr^2} \ \ 5\cos^3\theta-3\cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*6;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(phi).*(5*cos(th).^2-1)/4;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{1}_3) \ \ f_{xr^2-5xz^2} \ \ \cos\phi(5\cos^2\theta-1)$$','interpreter','latex','FontSize',18);
Nn=N*6;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=sin(phi).*(5*cos(th).^2-1)/4;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{1}_3) \ \ f_{yr^2-5yz^2} \ \ \sin\phi(5\cos^2\theta-1)$$','interpreter','latex','FontSize',18);
Nn=N*6;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(2*phi).*sin(th).^2.*cos(th)/0.4;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{2}_3) \ \ f_{x^2z-y^2z} \ \ \cos2\phi\sin^2\theta\cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*6;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=sin(2*phi).*sin(th).^2.*cos(th)/0.4;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{2}_3) \ \ f_{xyz} \ \ \sin2\phi\sin^2\theta\cos\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=cos(3*phi).*sin(th).^3;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Re(Y^{3}_3) \ \ f_{x^3-3xy^2} \ \ \cos3\phi\sin^3\theta$$','interpreter','latex','FontSize',18);
Nn=N*3;
theta=2*rand(Nn,1)-1;
theta=asin(theta);
th=theta+pi/2;%補正用
phi=2*pi*rand(Nn,1);
t=rand(Nn,1);
k=sin(3*phi).*sin(th).^3;
t=find(t<abs(k));% 注意!
theta=theta(t(1:N));
phi=phi(t(1:N));
r=(randn(N,1)/6+rd).^(1/3);
[x,y,z]=sph2cart(phi,theta,r);
figure('color',[1,1,1]);
plot3(x,y,z,'b.');
axis equal;axis tight;
rotate3d on;
title('$$ Im(Y^{3}_3) \ \ f_{y^3-3x^2y} \ \ \sin3\phi\sin^3\theta$$','interpreter','latex','FontSize',18);

乱数分布から適度な距離を持った任意の数の点を選択

いろんな方法があると思いますが、statistics toolboxのclusterdata関数を使ってみました。 中心からの距離が、正規分布で平均1、標準偏差6σの半径とします。 概ね球面上に乱数で1000個の点を設け、11個のクラスターに分割し、その中心を表示させます。

clear;close all;
sc=0.3;
[sx,sy,sz]=sphere(30);
N=1000;rd=1;
rx=rd*sx;ry=rd*sy;rz=rd*sz;
sx=sc*sx;sy=sc*sy;sz=sc*sz;
r=randn(N*3,1)/6;
r(find(abs(r)>=1))=[];r=r(1:N)+rd;
figure('color',[1,1,1],'InvertHardCopy','off');
Nn=N*2;
theta=asin(2*rand(Nn,1)-1);
phi=2*pi*rand(Nn,1);
theta=theta(1:N);
phi=phi(1:N);
[x,y,z]=sph2cart(phi,theta,r);
surface(rx,ry,rz,'facecolor',[1,0,0],'facealpha',0.1,'edgecolor','none');hold on;
seg=11;
C=clusterdata([x,y,z],'linkage','ward','savememory','on','maxclust',seg);
map=jet(seg);
for n=1:seg
    c=find(C==n);
    cx=x(c);cy=y(c);cz=z(c);
    cc=[cx,cy,cz]-ones(length(c),1)*mean([cx,cy,cz],1);
    [~,cc]=min(sum(cc.^2,2));
    h=plot3(cx,cy,cz,'.');
    set(h,'color',map(n,:));
    h=surface(sx+cx(cc),sy+cy(cc),sz+cz(cc),'facecolor',map(n,:),'edgecolor','none');
end;
axis equal;axis tight;
rotate3d on;
light('position',100*[1,1,1]);
light('position',-100*[1,1,1]);

原子核の殻模型

原子の殻模型(電子の分布)のように、原子核にも殻模型があるとされます。 原子核の殻模型の場合は原子の殻模型よりも複雑で、各核子の軌道角運動量(公転)とスピン角運動量(自転)の向きを考慮したものになっています。 核子別、すなわち陽子・中性子別の殻模型になっています。 埋まり方は、1s1/2、1p3/2、1p1/2、1d5/2、2s1/2、1d3/2、1f7/2、2p3/2、1f/5/2、2p1/2、1g9/2、1g7/2、2d5/2の順になっています。 動径分布なるものがあるのですが、無視して、適当な半径で3σの正規分布になると仮定します。

He-4は1s1/2に中性子2、陽子2個です。

clear;close all;
sc=0.5;
[sx,sy,sz]=sphere(30);
N=1000;rd=1;
rx=rd*sx;ry=rd*sy;rz=rd*sz;
sx=sc*sx;sy=sc*sy;sz=sc*sz;
r=randn(N*3,1)/3;
r(find(abs(r)>=1))=[];r=r(1:N)+rd;
figure('color',[1,1,1],'InvertHardCopy','off');
Nn=N*2;
theta=asin(2*rand(Nn,1)-1);
phi=2*pi*rand(Nn,1);
theta=theta(1:N);
phi=phi(1:N);
[x,y,z]=sph2cart(phi,theta,r);
surface(rx,ry,rz,'facecolor',[1,0,0],'facealpha',0.1,'edgecolor','none');hold on;
seg=4;
C=clusterdata([x,y,z],'linkage','ward','savememory','on','maxclust',seg);
map=jet(seg);
for n=1:seg
    c=find(C==n);
    cx=x(c);cy=y(c);cz=z(c);
    cc=[cx,cy,cz]-ones(length(c),1)*mean([cx,cy,cz],1);
    [~,cc]=min(sum(cc.^2,2));
    h=surface(sx+cx(cc),sy+cy(cc),sz+cz(cc),'facecolor',[0,0,1],'edgecolor','none','facelighting','phong');
    if n>2;set(h,'facecolor',[1,0,0]);end;
end;
axis equal;axis tight;axis off;
rotate3d on;
light('position',100*[1,1,1]);
light('position',-100*[1,1,1]);
title('$$^4_2$$He','interpreter','latex','FontSize',18);

C-12は1s1/2に中性子2、陽子2、1p3/2に中性子4、陽子4個です。 核子同士が重ならないように乱数で作成した点で核子周囲の点はナシにしました。 乱数の無駄遣いのような気もします。

clear;close all;
sc=0.5;
[sx,sy,sz]=sphere(30);
sx=sc*sx;sy=sc*sy;sz=sc*sz;
N=1000;
figure('color',[1,1,1],'InvertHardCopy','off');
map=[1,0,0;0,0,1];
nuc=zeros(12,3);
nucno=1;
blue=[0,0,1];
red=[1,0,0];
for k=1:4
    switch k
        case 1;rd=1;Nn=N*3;
        case 2;rd=1.5;Nn=N*3;
        case 3;rd=1.5;Nn=N*3;
        case 4;rd=1.5;Nn=N*3;
    end;

    rds=rd/sc;
    rx=rds*sx;ry=rds*sy;rz=rds*sz;

    if ismember(k,[1,2])
        surface(rx,ry,rz,'facecolor',map(k,:),'facealpha',0.1,'edgecolor','none');hold on;
    end;
    theta=asin(2*rand(Nn,1)-1);
    th=theta+pi/2;% 補正用
    phi=2*pi*rand(Nn,1);
    t=rand(Nn,1);
    switch k
        case 1;% 1s
            kt=10000;
            col=[blue;blue;red;red];
        case 2;% 1pz
            kt=cos(th);
            col=[blue;blue;red;red];
        case 3;% 1px
            kt=cos(phi).*sin(th);
            col=[blue;red];
        case 4;% 1py
            kt=sin(phi).*sin(th);
            col=[blue;red];
    end;
    seg=size(col,1);
    t=find(t<abs(kt));
    theta=theta(t);
    phi=phi(t);
    r=randn(round(length(t)*1.2),1)/3;
    r(find(abs(r)>=1))=[];
    r=r(1:length(t))+rd;
    [x,y,z]=sph2cart(phi,theta,r);
    for kk=1:size(nuc,1)
        dist=((x-nuc(kk,1)).^2+(y-nuc(kk,2)).^2+(z-nuc(kk,3)).^2).^0.5;
        ndist=find(dist>2*sc);
        x=x(ndist);
        y=y(ndist);
        z=z(ndist);
    end;
    nt=length(x);
    if nt<seg*100;return;end;
    C=clusterdata([x,y,z],'linkage','ward','savememory','on','maxclust',seg);

    for n=1:seg
        c=find(C==n);
        cx=x(c);cy=y(c);cz=z(c);
        cc=[cx,cy,cz]-ones(length(c),1)*mean([cx,cy,cz],1);
        [~,cc]=min(sum(cc.^2,2));
        nuc(nucno,:)=[cx(cc),cy(cc),cz(cc)];
        nucno=nucno+1;
        h=surface(sx+cx(cc),sy+cy(cc),sz+cz(cc),'facecolor',col(n,:),'edgecolor','none','facelighting','phong');
    end;
end;
axis equal;axis tight;axis off;
rotate3d on;
light('position',100*[1,1,1]);
light('position',-100*[1,1,1]);
title('$$^{12}_{6}$$C','interpreter','latex','FontSize',18);

Al-27は1s1/2に中性子2、陽子2、1p3/2に中性子4、陽子4、1p1/2に中性子2、陽子2、1d5/2に中性子6、陽子5個です。 コードの追加・上書きになります。半径のrdは適当です。

clear;close all;
sc=0.5;
[sx,sy,sz]=sphere(30);
sx=sc*sx;sy=sc*sy;sz=sc*sz;
N=1000;
figure('color',[1,1,1],'InvertHardCopy','off');
map=flipud(jet(4));
mapno=1;
nuc=zeros(27,3);
nucno=1;
blue=[0,0,1];
red=[1,0,0];
for k=1:11
    switch k
        case 1;rd=1;Nn=N*3;
        case {2,3,4};rd=1.5;Nn=N*3;
        case {5,6};rd=1.7;Nn=N*3;
        case {7,8,9,10,11};rd=2.2;Nn=N*3;
    end;

    rds=rd/sc;
    rx=rds*sx;ry=rds*sy;rz=rds*sz;

    if ismember(k,[1,2,5,7])
        surface(rx,ry,rz,'facecolor',map(mapno,:),'facealpha',0.1,'edgecolor','none');hold on;
        mapno=mapno+1;
    end;
    theta=asin(2*rand(Nn,1)-1);
    th=theta+pi/2;% 補正用
    phi=2*pi*rand(Nn,1);
    t=rand(Nn,1);
    switch k
        case 1;% 1s1/2
            kt=10000;
            col=[blue;blue;red;red];
        case 2;% 1pz3/2
            kt=cos(th);
            col=[blue;blue;red;red];
        case {3,5};% 1px3/2 1/2
            kt=cos(phi).*sin(th);
            col=[blue;red];
        case {4,6};% 1py3/2 1/2
            kt=sin(phi).*sin(th);
            col=[blue;red];
        case 7;% 1d3z^2-r^2
            kt=cos(th).^2;
            col=[blue;blue;red;red];
        case 8;% 1dxz
            kt=cos(phi).*sin(th).*cos(th);
            col=[blue;red];
        case 9;% 1dyz
            kt=sin(phi).*sin(th).*cos(th);
            col=[blue;red];
        case 10;% 1dx^2-y^2
            kt=cos(2*phi).*sin(th).^2;
            col=[blue;red];
        case 11;% 1dxy
            kt=sin(2*phi).*sin(th).^2;
            col=blue;

    end;
    seg=size(col,1);
    t=find(t<abs(kt)/max(abs(kt(:))));
    theta=theta(t);
    phi=phi(t);
    r=randn(round(length(t)*1.2),1)/3;
    r(find(abs(r)>=1))=[];
    r=r(1:length(t))+rd;
    [x,y,z]=sph2cart(phi,theta,r);
    for kk=1:size(nuc,1)
        dist=((x-nuc(kk,1)).^2+(y-nuc(kk,2)).^2+(z-nuc(kk,3)).^2).^0.5;
        ndist=find(dist>2*sc);
        x=x(ndist);
        y=y(ndist);
        z=z(ndist);
    end;
    nt=length(x);
    if nt<seg*100;return;end;
    C=clusterdata([x,y,z],'linkage','ward','savememory','on','maxclust',seg);

    for n=1:seg
        c=find(C==n);
        cx=x(c);cy=y(c);cz=z(c);
        cc=[cx,cy,cz]-ones(length(c),1)*mean([cx,cy,cz],1);
        [~,cc]=min(sum(cc.^2,2));
        nuc(nucno,:)=[cx(cc),cy(cc),cz(cc)];
        nucno=nucno+1;
        h=surface(sx+cx(cc),sy+cy(cc),sz+cz(cc),'facecolor',col(n,:),'edgecolor','none','facelighting','phong');
    end;
end;
axis equal;axis tight;axis off;
rotate3d on;
light('position',100*[1,1,1]);
light('position',-100*[1,1,1]);
title('$$^{27}_{13}$$Al','interpreter','latex','FontSize',18);

Co-60は1s1/2に中性子2、陽子2、1p3/2に中性子4、陽子4、1p1/2に中性子2、陽子2、1d5/2に中性子6、陽子6、 2s1/2に中性子2、陽子2、1d3/2に中性子4、陽子4、1f7/2に中性子8、陽子7、2p3/2に中性子4、1f/2に中性子1個です。 コードの追加・上書きになります。半径のrdは適当です。

clear;close all;
sc=0.5;
[sx,sy,sz]=sphere(30);
sx=sc*sx;sy=sc*sy;sz=sc*sz;
N=1000;
figure('color',[1,1,1],'InvertHardCopy','off');
map=flipud(jet(9));
mapno=1;
nuc=zeros(60,3);
nucno=1;
blue=[0,0,1];
red=[1,0,0];
for k=1:27
    switch k
        case 1;rd=1;Nn=N*3;
        case {2,3,4};rd=1.5;Nn=N*3;
        case {5,6};rd=1.7;Nn=N*3;
        case {7,8,9,10,11};rd=2.2;Nn=N*3;
        case 12;rd=2.5;Nn=N*3;
        case {13,14,15,16};rd=2.7;Nn=N*3;
        case {17,18,19,20,21,22,23};rd=2.9;Nn=N*3;
        case {24,25,26};rd=3.1;Nn=N*3;
        case 27;rd=3.3;Nn=N*3;
    end;

    rds=rd/sc;
    rx=rds*sx;ry=rds*sy;rz=rds*sz;

    if ismember(k,[1,2,5,7,12,13,17,24,27])
        surface(rx,ry,rz,'facecolor',map(mapno,:),'facealpha',0.1,'edgecolor','none');hold on;
        mapno=mapno+1;
    end;
    theta=asin(2*rand(Nn,1)-1);
    th=theta+pi/2;% 補正用
    phi=2*pi*rand(Nn,1);
    t=rand(Nn,1);
    switch k
        case 1;% 1s1/2
            kt=10000;
            col=[blue;blue;red;red];
        case 2;% 1pz3/2
            kt=cos(th);
            col=[blue;blue;red;red];
        case {3,5};% 1px3/2 1/2
            kt=cos(phi).*sin(th);
            col=[blue;red];
        case {4,6};% 1py3/2 1/2
            kt=sin(phi).*sin(th);
            col=[blue;red];
        case 7;% 1d3z^2-r^2 5/2
            kt=cos(th).^2;
            col=[blue;blue;red;red];
        case {8,13};% 1dxz 5/2 3/2
            kt=cos(phi).*sin(th).*cos(th);
            col=[blue;red];
        case {9,14};% 1dyz 5/2 3/2
            kt=sin(phi).*sin(th).*cos(th);
            col=[blue;red];
        case {10,15};% 1dx^2-y^2 5/2 3/2
            kt=cos(2*phi).*sin(th).^2;
            col=[blue;red];
        case {11,16};% 1dxy 5/2 3/2
            kt=sin(2*phi).*sin(th).^2;
            col=[blue;red];
        case 12;% 2s1/2
            kt=10000;
            col=[blue;blue;red;red];
        case 17;% 1f5z^3-3zr^2 7/2
            kt=5*cos(th).^3-3*cos(th);
            col=[blue;blue;red;red];
        case {18,27};% 1fxr^2-5xz^2
            kt=cos(phi).*(5*cos(th).^2-1);
            col=[blue;red];if k==27;col=blue;end;
        case 19;% 1f5yr^2-ryz^2
            kt=sin(phi).*(5*cos(th).^2-1);
            col=[blue;red];
        case 20;% 1f5x2z-y2z
            kt=cos(2*phi).*sin(th).^2.*cos(th);
            col=[blue;red];
        case 21;% 15fxyz
            kt=sin(2*phi).*sin(th).^2.*cos(th);
            col=[blue;red];
        case 22;% 15fx^3-3xy^2
            kt=cos(3*phi).*sin(th).^3;
            col=[blue;red];
        case 23;% 15fy^3-3x^2y
            kt=sin(3*phi).*sin(th).^3;
            col=blue;
        case 24;% 2pz3/2
             kt=10000;
            col=[blue;blue];
        case 25;% 2px3/2 1/2
            kt=cos(phi).*sin(th);
            col=blue;
        case 26;% 2py3/2 1/2
            kt=sin(phi).*sin(th);
            col=blue;

    end;
    seg=size(col,1);
    t=find(t<abs(kt)/max(abs(kt(:))));
    theta=theta(t);
    phi=phi(t);
    r=randn(round(length(t)*1.2),1)/3;
    r(find(abs(r)>=1))=[];
    r=r(1:length(t))+rd;
    [x,y,z]=sph2cart(phi,theta,r);
    for kk=1:size(nuc,1)
        dist=((x-nuc(kk,1)).^2+(y-nuc(kk,2)).^2+(z-nuc(kk,3)).^2).^0.5;
        ndist=find(dist>2*sc);
        x=x(ndist);
        y=y(ndist);
        z=z(ndist);
    end;
    nt=length(x);
    if nt<seg*100;return;end;
    C=clusterdata([x,y,z],'linkage','ward','savememory','on','maxclust',seg);

    for n=1:seg
        c=find(C==n);
        cx=x(c);cy=y(c);cz=z(c);
        cc=[cx,cy,cz]-ones(length(c),1)*mean([cx,cy,cz],1);
        [~,cc]=min(sum(cc.^2,2));
        nuc(nucno,:)=[cx(cc),cy(cc),cz(cc)];
        nucno=nucno+1;
        h=surface(sx+cx(cc),sy+cy(cc),sz+cz(cc),'facecolor',col(n,:),'edgecolor','none','facelighting','phong');
    end;
end;
axis equal;axis tight;axis off;
rotate3d on;
light('position',100*[1,1,1]);
light('position',-100*[1,1,1]);
title('$$^{60}_{27}$$Co','interpreter','latex','FontSize',18);