Contents

原子核の3D

原子核は電子の軌道のように球面調和関数に似た殻構造をもっているらしいのですが、面倒なので1つの球が12個の球に囲まれているとします。 デタラメに中性子と陽子を並べたものです。

close all;clear;
N=30;
NN=ones(N+1,N+1);
[x,y,z]=sphere(N);
a=3^0.5;b=(8/3)^0.5;aa=1;
x=aa*x;
y=aa*y;
z=aa*z;
dx=[2,0,0];
dy=[1,a,0];
dz=[1,1/a,b];
dzz=[1,1/a,-b];%右手・左手
for Z=1:6
    switch Z
        case 1;
            kk=[0,0,0;1,0,0;0,1,0;0,0,1];
            st=['$$^{4} _{2}$$','He'];
            nnn=4-2;
            DZplus=dz;
        case 2;
            kk=[0,0,0;1,0,0;2,0,0;0,1,0;1,1,0;0,2,0;...
                0,0,1;1,0,1;0,1,1;...
                0,0,-1;1,0,-1;0,1,-1];
            st=['$$^{12} _{6}$$','C'];
            nnn=12-6;
            DZplus=dz;
            DZminus=dzz;
        case 3;
            kk=[0,0,0;1,0,0;2,0,0;3,0,0; 0,1,0;1,1,0;2,1,0; 1,-1,0;2,-1,0;3 -1,0; 2,-2,0;3,-2,0;...
                0,0,1;1,0,1;2,0,1;1,-1,1;2,-1,1;2,-2,1;...
                1,0,2;2,0,2;2,-1,2;...
                0,0,-1;1,0,-1;2,0,-1;1,-1,-1;2,-1,-1;2,-2,-1];
            st=['$$^{27} _{13}$$','Al'];
            nnn=27-13;
            DZplus=cumsum([dz;-dzz]);
            DZminus=dzz;
        case 4;
            kk=[0,0,0;1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 0,3,0;1,3,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0;...
                0,0,1;1,0,1;2,0,1;3,0,1; 0,1,1;1,1,1;2,1,1; 0,2,1;1,2,1; 1,-1,1;2,-1,1;3,-1,1;...
                1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2;1,2,2;...
                1,0,3;2,0,3;1,1,3;...
                0,0,-1;1,0,-1;2,0,-1;3,0,-1; 0,1,-1;1,1,-1;2,1,-1; 0,2,-1;1,2,-1; 1,-1,-1;2,-1,-1;3,-1,-1];
            st=['$$^{51} _{23}$$','V'];
            nnn=51-23;
            DZplus=cumsum([dz;-dzz;dz]);
            DZminus=dzz;
        case 5;
            kk=[0,0,0;1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0; 2,-2,0;3,-2,0;4,-2,0;...
                0,0,1;1,0,1;2,0,1;3,0,1;0,1,1;1,1,1;2,1,1;1,-1,1;2,-1,1;3,-1,1; 2,-2,1;3,-2,1;...
                1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2; 2,-1,2;3,-1,2;...
                0,0,-1;1,0,-1;2,0,-1;3,0,-1;0,1,-1;1,1,-1;2,1,-1;1,-1,-1;2,-1,-1;3,-1,-1; 2,-2,-1;3,-2,-1;...
                1,0,-2;2,0,-2;3,0,-2; 1,1,-2;2,1,-2; 2,-1,-2;3,-1,-2];
            st=['$$^{57} _{26}$$','Fe'];
            nnn=57-26;
            DZplus=cumsum([dz;-dzz]);
            DZminus=cumsum([dzz;-dz]);
        case 6;
            kk=[0,0,0;1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0; 2,-2,0;3,-2,0;4,-2,0;...
                0,0,1;1,0,1;2,0,1;3,0,1;0,1,1;1,1,1;2,1,1;1,-1,1;2,-1,1;3,-1,1; 2,-2,1;3,-2,1;...
                1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2; 2,-1,2;3,-1,2;...
                1,0,3;2,0,3;2,-1,3;...
                0,0,-1;1,0,-1;2,0,-1;3,0,-1;0,1,-1;1,1,-1;2,1,-1;1,-1,-1;2,-1,-1;3,-1,-1; 2,-2,-1;3,-2,-1;...
                1,0,-2;2,0,-2;3,0,-2; 1,1,-2;2,1,-2; 2,-1,-2;3,-1,-2];
            % st=sprintf('N=%d',size(kk,1));
            st=['$$^{60} _{28}$$','Ni'];
            nnn=60-28;
            DZplus=cumsum([dz;-dzz;dz]);
            DZminus=cumsum([dzz;-dz]);
    end;
    k=zeros(size(kk,1),3);
    for n=1:size(kk,1)

        k(n,:)=kk(n,1)*dx+kk(n,2)*dy;
        kz=kk(n,3);
        if kz>0
            k(n,:)=k(n,:)+DZplus(kz,:);
        elseif kz<0
            k(n,:)=k(n,:)+DZminus(abs(kz),:);
        end;
    end;

    mk=mean(k,1);
    figure('color',[1,1,1],'InvertHardCopy','off');
    for n=1:size(k,1)
        surface(x+(k(n,1)-mk(1))*NN,...
                y+(k(n,2)-mk(2))*NN,...
                z+(k(n,3)-mk(3))*NN,...
                'tag','p');
        if n==1;hold on;end;
    end;
    daspect([1,1,1]);
    hp=findobj(gcf,'tag','p');
    axis tight;axis off;view([60,10]);
    set(hp,'facecolor',[1,0,0],'edgecolor','none');
    set(hp,'facelighting','phong');
    [~,nn]=sort(rand(1,length(hp)));
    set(hp(nn(1:nnn)),'facecolor',[0,0,1]);
    light('position',100*[1,1,1]);
    light('position',-100*[1,1,1]);
    rotate3d on;
    title(st,'interpreter','latex','FontSize',15);
end;

中性子と陽子を交換しながら自転するCo59...GIF動画

陽子⇔中性子になる頻度はわかりませんが、陽子⇔中性子になりながら自転するCo59を作成してみました。

close all;clear;
seg=36;
ang=linspace(0,2*pi,seg+1);
N=30;
NN=ones(N+1,N+1);
[x,y,z]=sphere(N);
a=3^0.5;b=(8/3)^0.5;aa=1;
GIF=2;% GIF==0だと動画は作成されない。
kk=[1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0; 2,-2,0;3,-2,0;4,-2,0;...
    0,0,1;1,0,1;2,0,1;3,0,1;0,1,1;1,1,1;2,1,1;1,-1,1;2,-1,1;3,-1,1; 2,-2,1;3,-2,1;...
    1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2; 2,-1,2;3,-1,2;...
    1,0,3;2,0,3;2,-1,3;...
    0,0,-1;1,0,-1;2,0,-1;3,0,-1;0,1,-1;1,1,-1;2,1,-1;1,-1,-1;2,-1,-1;3,-1,-1; 2,-2,-1;3,-2,-1;...
    1,0,-2;2,0,-2;3,0,-2; 1,1,-2;2,1,-2; 2,-1,-2;3,-1,-2];
st=['$$^{59} _{27}$$','Co'];
nnn=59-27;

for t=1:seg
    cosA=cos(ang(t));sinA=sin(ang(t));
    X=[cosA,-sinA,0;sinA,cosA,0;0,0,1];
    dx=[2,0,0]*X;
    dy=[1,a,0]*X;
    dz=[1,1/a,b]*X;
    dzz=[1,1/a,-b]*X;
    DZplus=cumsum([dz;-dzz;dz]);
    DZminus=cumsum([dzz;-dz]);

    k=zeros(size(kk,1),3);
    for n=1:size(kk,1)
        k(n,:)=kk(n,1)*dx+kk(n,2)*dy;
        kz=kk(n,3);
        if kz>0
            k(n,:)=k(n,:)+DZplus(kz,:);
        elseif kz<0
            k(n,:)=k(n,:)+DZminus(abs(kz),:);
        end;
    end;
    mk=mean(k,1);
    switch GIF
        case 0;% とりあえずこれ
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
            else
                cla;
            end;
        case 1;% GIF動画作成用
            close(gcf);
            figure('color',[1,1,1],'InvertHardCopy','off');
            fname='Co59.gif';
        case 2;% ホームページ用
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
            end;
            subplot('position',[mod(t-1,6)/6+0.05,(5-floor((t-1)/6))/6-0.05,1/6-0.1,1/6]);
    end;

    for n=1:size(k,1)
        surface(x+(k(n,1)-mk(1))*NN,...
                y+(k(n,2)-mk(2))*NN,...
                z+(k(n,3)-mk(3))*NN,...
                'tag','p');
        if n==1;hold on;end;
    end;
    daspect([1,1,1]);
    hp=findobj(gca,'tag','p');
    axis tight;axis off;
    view([-140,10]);
    set(hp,'facecolor',[1,0,0],'edgecolor','none');
    set(hp,'facelighting','phong');
    [~,nn]=sort(rand(1,length(hp)));
    set(hp(nn(1:nnn)),'facecolor',[0,0,1]);
    set(gca,'xlim',6*[-1,1],'ylim',6*[-1,1],'zlim',[-6,6]);
    light('position',100*[1,1,1]);
    light('position',-100*[1,1,1]);
    if ismember(GIF,[0,1])
        text(2,1,-7,st,'interpreter','latex','FontSize',30);
        set(gcf,'position',[50,50,300,400]);
        set(gca,'position',[0,0,1,1]);
        set(gca,'cameraviewangle',7);
    else
        title(st,'interpreter','latex','FontSize',15);
    end;
    switch GIF
        case 0;
            pause(0.05);
            drawnow;
        case 1;
            F=getframe(gcf);
            [fc,fm]=rgb2ind(F.cdata,256);
            if t==1
                imwrite(fc,fm,fname,'loopcount',inf,'delaytime',0);
            else
                imwrite(fc,fm,fname,'writemode','append','delaytime',0);
            end;
    end;
end;
rotate3d on;

Co59がCo60になるところ

Co60はCo59に中性子が当てて作成します。

close all;clear;
seg=36*2;
ang=linspace(0,4*pi,seg+1);
N=30;
NN=ones(N+1,N+1);
[x,y,z]=sphere(N);
a=3^0.5;b=(8/3)^0.5;aa=1;
GIF=2;% GIF==0だと動画は作成されない。
if GIF==1;seg=seg/2;end;
kk=[1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0; 2,-2,0;3,-2,0;4,-2,0;...
    0,0,1;1,0,1;2,0,1;3,0,1;0,1,1;1,1,1;2,1,1;1,-1,1;2,-1,1;3,-1,1; 2,-2,1;3,-2,1;...
    1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2; 2,-1,2;3,-1,2;...
    1,0,3;2,0,3;2,-1,3;...
    0,0,-1;1,0,-1;2,0,-1;3,0,-1;0,1,-1;1,1,-1;2,1,-1;1,-1,-1;2,-1,-1;3,-1,-1; 2,-2,-1;3,-2,-1;...
    1,0,-2;2,0,-2;3,0,-2; 1,1,-2;2,1,-2; 2,-1,-2;3,-1,-2];
st1=['$$^{59} _{27}$$','Co'];
st2=['$$^{60} _{27}$$','Co'];
nnn=59-27;
check=0;
st=st1;
tx=linspace(-25,-5,seg/2);
for t=1:seg
    if check==0
        if t>=seg/2;
            kk=[0,0,0;kk];
            nnn=nnn+1;
            st=st2;
            check=1;
        end;
    end;
    cosA=cos(ang(t));sinA=sin(ang(t));
    X=[cosA,-sinA,0;sinA,cosA,0;0,0,1];
    dx=[2,0,0]*X;
    dy=[1,a,0]*X;
    dz=[1,1/a,b]*X;
    dzz=[1,1/a,-b]*X;
    DZplus=cumsum([dz;-dzz;dz]);
    DZminus=cumsum([dzz;-dz]);

    k=zeros(size(kk,1),3);
    for n=1:size(kk,1)
        k(n,:)=kk(n,1)*dx+kk(n,2)*dy;
        kz=kk(n,3);
        if kz>0
            k(n,:)=k(n,:)+DZplus(kz,:);
        elseif kz<0
            k(n,:)=k(n,:)+DZminus(abs(kz),:);
        end;
    end;
    mk=mean(k,1);
    switch GIF
        case 0;% とりあえずこれ
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
                set(gcf,'position',[50,50,800,400]);
            else
                cla;
            end;
        case 1;% GIF動画作成用
            close(gcf);
            figure('color',[1,1,1],'InvertHardCopy','off');
            set(gcf,'position',[50,50,800,400]);
            fname='Co60.gif';
        case 2;% ホームページ用
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
                set(gcf,'position',[50,50,1000,400]);
            end;
            subplot('position',[mod(t-1,6)/6+0.05,(5-floor((t-1)/6))/6-0.05,1/6-0.1,1/6]);
    end;

    for n=1:size(k,1)
        surface(x+(k(n,1)-mk(1))*NN,...
                y+(k(n,2)-mk(2))*NN,...
                z+(k(n,3)-mk(3))*NN,...
                'tag','p');
        if n==1;hold on;end;
    end;
    if check==0;
        hnt=surface(x+tx(t),y,z,'edgecolor','none','facelighting','phong','facecolor',[0,0,1]);
    end;
    daspect([1,1,1]);
    hp=findobj(gca,'tag','p');
    axis off;
    view([-140,10]);
    set(hp,'facecolor',[1,0,0],'edgecolor','none');
    set(hp,'facelighting','phong');
    [~,nn]=sort(rand(1,length(hp)));
    set(hp(nn(1:nnn)),'facecolor',[0,0,1]);
    set(gca,'xlim',[-12,4],'ylim',[-26,6],'zlim',[-7,5]);
    if ismember(GIF,[0,1]);
        set(gca,'position',[0,0,1,1]);
    end;
    light('position',100*[1,1,1]);
    light('position',-100*[1,1,1]);
    if ismember(GIF,[0,1])
        text(2,1,-7,st,'interpreter','latex','FontSize',30);
        set(gca,'cameraviewangle',5);
    else
        title(st,'interpreter','latex','FontSize',15);
    end;
    switch GIF
        case 0;
            pause(0.05);
            drawnow;
        case 1;
            F=getframe(gcf);
            [fc,fm]=rgb2ind(F.cdata,256);
            if t==1
                imwrite(fc,fm,fname,'loopcount',inf,'delaytime',0);
            else
                imwrite(fc,fm,fname,'writemode','append','delaytime',0);
            end;
    end;
end;
rotate3d on;

Co60のβマイナス崩壊

基本的にコードの使いまわしです。

close all;clear;
seg=36*3;
ang=linspace(0,4*pi,seg+1);
N=30;
NN=ones(N+1,N+1);
[x,y,z]=sphere(N);
a=3^0.5;b=(8/3)^0.5;aa=0.5;ab=0.3;ac=4.3;
GIF=2;% GIF==0だと動画は作成されない。
if GIF==2;seg=seg/3;end;
kk=[0,0,0;1,0,0;2,0,0;3,0,0;4,0,0; 0,1,0;1,1,0;2,1,0;3,1,0; 0,2,0;1,2,0;2,2,0; 1,-1,0;2,-1,0;3,-1,0;4,-1,0; 2,-2,0;3,-2,0;4,-2,0;...
    0,0,1;1,0,1;2,0,1;3,0,1;0,1,1;1,1,1;2,1,1;1,-1,1;2,-1,1;3,-1,1; 2,-2,1;3,-2,1;...
    1,0,2;2,0,2;3,0,2; 1,1,2;2,1,2; 2,-1,2;3,-1,2;...
    1,0,3;2,0,3;2,-1,3;...
    0,0,-1;1,0,-1;2,0,-1;3,0,-1;0,1,-1;1,1,-1;2,1,-1;1,-1,-1;2,-1,-1;3,-1,-1; 2,-2,-1;3,-2,-1;...
    1,0,-2;2,0,-2;3,0,-2; 1,1,-2;2,1,-2; 2,-1,-2;3,-1,-2];
st1=['$$^{60} _{27}$$','Co'];
st2=['$$^{60} _{28}$$','Ni','$$ ^{*}$$'];
st3=['$$^{60} _{28}$$','Ni'];
nnn=60-27;
check=0;
st=st1;
tx=linspace(-5,-25,seg/3);
ty=linspace(0,-15,seg/3);
tz=linspace(0,7,seg/3);
for t=1:seg
    switch check
        case 0;
            if t>seg/3
                nnn=nnn-1;
                % st=st2;
                check=1;
            end;
        case 1;
            if t>seg/3*2
                st=st3;
                check=2;
            end;
    end;

    cosA=cos(ang(t));sinA=sin(ang(t));
    X=[cosA,-sinA,0;sinA,cosA,0;0,0,1];
    dx=[2,0,0]*X;
    dy=[1,a,0]*X;
    dz=[1,1/a,b]*X;
    dzz=[1,1/a,-b]*X;
    DZplus=cumsum([dz;-dzz;dz]);
    DZminus=cumsum([dzz;-dz]);

    k=zeros(size(kk,1),3);
    for n=1:size(kk,1)
        k(n,:)=kk(n,1)*dx+kk(n,2)*dy;
        kz=kk(n,3);
        if kz>0
            k(n,:)=k(n,:)+DZplus(kz,:);
        elseif kz<0
            k(n,:)=k(n,:)+DZminus(abs(kz),:);
        end;
    end;
    mk=mean(k,1);
    switch GIF
        case 0;% とりあえずこれ
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
                set(gcf,'position',[50,50,800,400]);
            else
                cla;
            end;
        case 1;% GIF動画作成用
            close(gcf);
            figure('color',[1,1,1],'InvertHardCopy','off');
            set(gcf,'position',[50,50,800,400]);
            fname='betadecay.gif';
        case 2;% ホームページ用
            if t==1
                figure('color',[1,1,1],'InvertHardCopy','off');
                set(gcf,'position',[50,50,1000,400]);
            end;
            subplot('position',[mod(t-1,6)/6+0.05,(5-floor((t-1)/6))/6-0.05,1/6-0.1,1/6]);
    end;

    for n=1:size(k,1)
        surface(x+(k(n,1)-mk(1))*NN,...
                y+(k(n,2)-mk(2))*NN,...
                z+(k(n,3)-mk(3))*NN,...
                'tag','p');
        if n==1;hold on;end;
    end;
    if check==0
        if t+3>seg/3
            surface(ac*x-3,ac*y,ac*z,'edgecolor','none','facelighting','phong','facecolor',[0.5,0.5,0],'facealpha',0.7);
            st=st2;
            if ismember(GIF,[0,1]);
                text(-3,0,-ac-1,'$$W^{-}$$','interpreter','latex','fontsize',30,'tag','hW');
            end;
        end;
    elseif check==1;
        ts=tx(t-seg/3);
        tsx=ts*0.7;
        tsy=ty(t-seg/3);
        tsz=tz(t-seg/3);
        hnt=surface(aa*x+ts,aa*y,aa*z,'edgecolor','none','facelighting','phong','facecolor',[0,1,1]);
        hnn=surface(ab*x+tsx,ab*y+tsy,ab*z+tsz,'edgecolor','none','facelighting','phong','facecolor',[1,1,1]*0.5);
    elseif check==2
        if t<seg/3*2+3
             gt=linspace(0,-50,200);
             [gx,gy]=meshgrid(gt,-1:0.25:1);
             for tt=1:2
                 gz=0.2*sin(gx*200*pi);
                 if tt==1

                     gz=gz+ones(size(gx,1),1)*linspace(0,-5,200);
                     col=[0,1,1];
                 else
                     gz=gz+ones(size(gx,1),1)*linspace(0,5,200);
                     col=[0,0,1];
                 end;
                 gh=surf(gy,gx,gz);hold on;
                 set(gh,'EdgeColor',col);
                 set(gh,'facecolor','none');
             end;

             if ismember(GIF,[0,1]);
                 text(0,-25,2.5+1,'$$\gamma_{1}$$','interpreter','latex','FontSize',30);
                 text(0,-25,-2.5-1,'$$\gamma_{2}$$','interpreter','latex','FontSize',30);
             end;
        end;
    end;
    daspect([1,1,1]);
    hp=findobj(gca,'tag','p');
    axis off;
    view([-140,10]);
    set(hp,'facecolor',[1,0,0],'edgecolor','none');
    set(hp,'facelighting','phong');
    [~,nn]=sort(rand(1,length(hp)));
    set(hp(nn(1:nnn)),'facecolor',[0,0,1]);
    set(gca,'xlim',[-12,4],'ylim',[-26,6],'zlim',[-7,5]);

    light('position',100*[1,1,1]);
    light('position',-100*[1,1,1]);
    if ismember(GIF,[0,1])
        set(gca,'position',[0,0,1,1]);
        text(2,1,-7,st,'interpreter','latex','FontSize',30);
        set(gca,'cameraviewangle',5);
        if check==1;
            text(ts,0,-1,'$$e^{-}$$','interpreter','latex','fontsize',30);
            text(tsx,tsy,tsz+1,'$$\bar{\nu _{e}}$$','interpreter','latex','fontsize',30);
        end;
    else
        title(st,'interpreter','latex','FontSize',15);
    end;
    switch GIF
        case 0;
            pause(0.05);
            drawnow;
        case 1;
            F=getframe(gcf);
            [fc,fm]=rgb2ind(F.cdata,256);
            if t==1
                imwrite(fc,fm,fname,'loopcount',inf,'delaytime',0);
            else
                imwrite(fc,fm,fname,'writemode','append','delaytime',0);
            end;
    end;
end;
rotate3d on;