Co60の半減期のGIF動画

ガンマナイフ 4CにはCo60が201個あります。半減期は約5年で購入時1TBqですが、5年後には500GBq、10年後には250GBqになります。
球面上の乱数はMATLABの乱数発生器>球内の乱数をそのままパクりました。
clear;close all;
GIF=2;
Nrad=2^7;
num=36;
Nn=Nrad*num;
rng(0,'twister');
rvals=2*rand(Nn,1)-1;
elevation=asin(rvals);
azimuth=2*pi*rand(Nn,1);
[rx,ry,rz]=sph2cart(azimuth,elevation,ones(Nn,1));
rx=reshape(rx,[Nrad,num]);
ry=reshape(ry,[Nrad,num]);
rz=reshape(rz,[Nrad,num]);
rn=rand(Nrad,num);
figure('color',[1,1,1]*0,'InvertHardCopy','off');
set(gcf,'position',[50,50,900,300]);
[qx,qy,qz]=sphere(7);
qxx=qx(:);
qyy=qy(:);
qzz=qz(:);
dist=12;
rd=2;
rd2=rd*2.5;
nnum=1:num;
ang=linspace(0,2*pi,num+1);
for t=1:num
    switch GIF
        case 0;
            cla;
        case 1;
            close(gcf);
            figure('color',[1,1,1]*0,'InvertHardCopy','off');
            set(gcf,'position',[50,50,900,300]);
            fname='HalfLife.gif';
        case 2;
             subplot('position',[mod(t-1,6)/6+0.05,(5-floor((t-1)/6))/6,1/6-0.1,1/6]);
    end;
    tt=mod(nnum+(t-1),num);
    for ttt=1:5
        t0(ttt)=find(tt==ttt);
    end;
    qqx=qx*cos(ang(t))-qy*sin(ang(t));
    qqy=qx*sin(ang(t))+qy*cos(ang(t));
    for hl=-1:1
        surface(qqx+hl*dist,qqy,qz,'facecolor',[0,1,1],'edgecolor',[0,0.9,1],'tag','core');
        hold on;
        pk=rand(1,numel(qx));
        pk=find(pk>0.9);
        if numel(pk)>1
            plot3(qxx(pk)+hl*dist,qyy(pk),qzz(pk),'marker','pentagram','LineStyle','none',...
                'MarkerFacecolor',[1,1,0],'MarkerEdgeColor','none','MarkerSize',8);
        end;
        axis off;
        set(gca,'xlim',[-1,1].*(dist+rd2*1.5),'ylim',[-1,1]*rd2*2,'zlim',[-1,1]*rd2*2);
        if ismember(GIF,[0,1]);
            set(gca,'position',[0,0,1,1]);
            set(gca,'CameraViewAngle',3);
        else
            set(gca,'CameraViewAngle',4);
        end;

        view(0,0);
        daspect([1,1,1]);

        N2=1:(Nrad*0.5^(hl+1));
        for ttt=1:5
            t1=t0(ttt);
            RN=[rn(N2,t1),rn(N2,t1)+rd]+ttt;
            RN(RN>rd2)=rd2;
            px=[rx(N2,t1)*ones(1,2)].*RN;
            py=[ry(N2,t1)*ones(1,2)].*RN;
            pz=[rz(N2,t1)*ones(1,2)].*RN;
            plot3(px'+hl*dist,py',pz','tag','ray');
        end;
        hray=findobj(gca,'tag','ray');
        set(hray(1:2:end),'color',[0,0.5,1]);
        set(hray(2:2:end),'color',[0,0.7,1]);
    end;

    switch GIF
        case 0;
            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;
daspect([1,1,1]);
rotate3d on;