Contents
Co-60の作成
Co-60はCo-59に中性子を当てて作成します。 Co原子核は一応回転させてますが、乱数で核子の位置を決定してるので、回転してるのはよくわらない絵になっています。 時々乱数不足？で失敗するかもしれません。
clear;close all; sc=0.5; [sx,sy,sz]=sphere(30); sx=sc*sx;sy=sc*sy;sz=sc*sz; N=1000; hf=figure('color',[1,1,1],'InvertHardCopy','off'); htitle=title('$$^{59}_{27}$$Co','interpreter','latex','FontSize',18,'tag','title'); map=flipud(jet(9)); blue=[0,0,1]; red=[1,0,0]; lim=4*[-1,1]; limx=[-4,14]; xt=linspace(13,2.3,40); GIF=0; for tm=1:1:72 if GIF==1 close(gcf); figure('color',[1,1,1],'InvertHardCopy','off'); htitle=title('$$^{59}_{27}$$Co','interpreter','latex','FontSize',18,'tag','title'); fname='Co59_shell.gif'; end; ang=tm/18*pi; cla; mapno=1; nuc=zeros(60,3); if tm>35 set(htitle,'string','$$^{60}_{27}$$Co'); end; if tm<41 nuc(1,:)=[xt(tm),0,0]; nucno=2; else nucno=1; end; 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); if tm<41 && k==27;continue;end; 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;blue];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+ang,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; if tm<41 h=surface(sx+xt(tm),sy,sz,'facecolor',blue,'edgecolor','none','facelighting','phong'); end; view([0,0]); set(gca,'cameraviewangle',6.5,'position',[0,0,1,0.9]); set(gca,'xlim',limx,'ylim',lim,'zlim',lim); axis equal;axis off; light('position',100*[1,1,1]); light('position',-100*[1,1,1]); if GIF==0 pause(0.05); drawnow; if ismember(tm,[1,18,36,40,72]); ha=gca; hf1=figure('color',[1,1,1],'invertHardCopy','off'); copyobj(ha,hf1); figure(hf); end; else F=getframe(gcf); [fc,fm]=rgb2ind(F.cdata,256); if tm==1 imwrite(fc,fm,fname,'loopcount',inf,'delaytime',0); else imwrite(fc,fm,fname,'writemode','append','delaytime',0); end; end; end; rotate3d on;
Co-60のβ・γ崩壊
Co-60はβマイナス崩壊してNi*-60となります。β線（電子）と反電子ニュートリノを出します。 さらにNi*-60はγ崩壊します。γ線２本出します。 わかりにくいですが、陽子が1f5/2から2p3/2、2p3/2から1f7/2に軌道を変えるときにγ線がでることにしました。
clear;close all; sc=0.5;wsc=sc*4.3;esc=0.2;nsc=0.1; [sx,sy,sz]=sphere(30); wx=wsc*sx;wy=wsc*sy;wz=wsc*sz;% ウィークボゾン ex=esc*sx;ey=esc*sy;ez=esc*sz;% 電子 nx=nsc*sx;ny=nsc*sy;nz=nsc*sz;% ニュートリノ sx=sc*sx;sy=sc*sy;sz=sc*sz; N=1000; hf=figure('color',[1,1,1],'InvertHardCopy','off'); htitle=title('$$^{60}_{27}$$Co','interpreter','latex','FontSize',18,'tag','title'); map=flipud(jet(9)); blue=[0,0,1]; red=[1,0,0]; lim=4*[-1,1]; limx=[-4,14]; xte=linspace(3,13,10); xtn=linspace(3,13,5); xtg=linspace(0,5,100); GIF=0; for tm=1:(36*2) if GIF==1 close(gcf); figure('color',[1,1,1],'InvertHardCopy','off'); htitle=title('$$^{60}_{27}$$Co','interpreter','latex','FontSize',18,'tag','title'); fname='Co60_betadecay.gif'; end; ang=tm/18*pi; cla; mapno=1; nuc=zeros(60,3); nucno=1; 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 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; case 27;% if tm<25;% 1f5/2 N rd=3.3;Nn=N*3; kt=cos(phi).*(5*cos(th).^2-1); col=blue; elseif tm<40;% 1f5/2 P rd=3.3;Nn=N*3; kt=cos(phi).*(5*cos(th).^2-1); col=red; elseif tm<45% 2p3/2 P rd=3.1;Nn=N*3; kt=10000; col=red; else %1f7/2 P rd=2.9;Nn=N*3; kt=sin(phi).*sin(th); col=red; end; 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+ang,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; if tm>46 set(htitle,'string','$$^{60}_{28}$$Ni'); elseif tm>24 set(htitle,'string','$$^{60}_{28}$$Ni*'); end; if tm>24 && tm<28 pwx=3.3; h=surface(wx+pwx,wy,wz,'facecolor',[0.5,0.5,0],'facealpha',0.7,'edgecolor','none','facelighting','phong'); text(pwx,0,-4,'$$W^{-}$$','interpreter','latex','FontSize',18,'tag','text'); end; if tm>27 && tm<36 xx=xte(tm-27);zz=xx*tan(-15/180*pi); h=surface(ex+xx,ey,ez+zz,'facecolor',[0,1,1],'edgecolor','none','facelighting','phong'); text(xx,0,zz-0.5,'$$\beta^{-}(e^{-})$$','interpreter','latex','fontsize',18); end; if tm>27 && tm<33 xx=xtn(tm-27);zz=xx*tan(10/180*pi); h=surface(nx+xx,ny,nz+zz,'facecolor',[1,1,1],'edgecolor','none','facelighting','phong'); text(xx,0,zz+0.5,'$$\bar{\nu}_{e}$$','interpreter','latex','fontsize',18); end; if tm>39 && tm<43 gxx=xtg+(tm-40)*3+3; gyy=0*xtg; gzz=0.5*sin(xtg*10); ang=5/180*pi; gx=gxx*cos(ang)-gzz*sin(ang); gy=gyy; gz=gxx*sin(ang)+gzz*cos(ang); plot3(gx,gy,gz,'color',[0,0.5,1],'LineWidth',2); text(mean(gx),0,mean(gz)+1,'$$\gamma_1$$','interpreter','latex','fontsize',18); end; if tm>45 && tm<49 gxx=xtg+(tm-46)*4+3; gyy=0*xtg; gzz=0.5*sin(xtg*10); ang=-5/180*pi; gx=gxx*cos(ang)-gzz*sin(ang); gy=gyy; gz=gxx*sin(ang)+gzz*cos(ang); plot3(gx,gy,gz,'color',[0,1,1],'LineWidth',2); text(mean(gx),0,mean(gz)+1,'$$\gamma_2$$','interpreter','latex','fontsize',18); end; view([0,0]); set(gca,'cameraviewangle',6.5,'position',[0,0,1,0.9]); set(gca,'xlim',limx,'ylim',lim,'zlim',lim); axis equal;axis off; light('position',100*[1,1,1]); light('position',-100*[1,1,1]); if GIF==0 pause(0.05); drawnow; if ismember(tm,[1,25,30,41,47,60]); ha=gca; hf1=figure('color',[1,1,1],'invertHardCopy','off'); copyobj(ha,hf1); figure(hf); end; else F=getframe(gcf); [fc,fm]=rgb2ind(F.cdata,256); if tm==1 imwrite(fc,fm,fname,'loopcount',inf,'delaytime',0); else imwrite(fc,fm,fname,'writemode','append','delaytime',0); end; end; end; rotate3d on;
