【物理应用】Matlab 带电粒子在放射状电场和均匀磁场下的混沌运动模拟

paths='D:\0 Matlab\20 Matlab 物理应用\带电粒子在放射状电场和均匀磁场下的混沌运动模拟\';% path for saving

% initial values
ix=1;
iy=0;
iz=0;
hozdis=[0 1];
verdis=[-5 5]
vxymag=[0 2];
vzrang=[-1 1];
alpha=[0 360];
step=10;
stepz=10;
stepvxy=10;
stepvz=10;
stepang=10;
timerg=[0 10000];

listm((step+1)^4,8)=0;
k=1;

for nz=0:stepz %z
    zstp=(verdis(2)-verdis(1))/stepz;
    iz=verdis(1)+zstp*nz;
    
    for nvxy=0:stepvxy %vxy
        vxystp=(vxymag(2)-vxymag(1))/stepvxy;
        vxy=vxymag(1)+nvxy*vxystp;
        
        for nvz=0:stepvz %vz
            vzstp=(vzrang(2)-vzrang(1))/stepvz;
            vz=vzrang(1)+nvz*vzstp;
            
            for xyangl=0:stepang %angle
                angstp=(alpha(2)-alpha(1))/stepang;
                vx=vxy*cos(pi*(alpha(1)+angstp*xyangl)/180);
                vy=vxy*sin(pi*(alpha(1)+angstp*xyangl)/180);
                
                [t,x]=ode45('test',timerg,[ix,iy,iz,vx,vy,vz]);
                
                %judging mechanism
                flag1=0;
                flag2=0;
                lastloc=size([t,x]);
                lastdis=sqrt((x(lastloc(1),1))^2+(x(lastloc(1),2))^2+(x(lastloc(1),3))^2);
                
                if abs(max(abs(x(:,3))))<=200
                    flag1=1;
                end
                if abs(max(abs(x(:,3))))<=100
                    flag1=2;
                end
                if abs(max(abs(x(:,3))))<=60
                    flag1=3;
                end
                if abs(max(abs(x(:,3))))<=30
                    flag1=4;
                end
                if abs(max(abs(x(:,3))))<=10
                    flag1=5;
                end
                if abs(max(abs(x(:,3))))<=3
                    flag1=6;
                end
                if abs(max(abs(x(:,3))))<=2
                    flag1=7;
                end
                if abs(max(abs(x(:,3))))<=1
                    flag1=8;
                end
                if abs(max(abs(x(:,3))))<=0.5
                    flag1=9;
                end
                if lastdis<=1
                    flag2=1;
                end
                
                % output list
                listm(k,1)=ix;
                listm(k,2)=iy;
                listm(k,3)=iz;
                listm(k,4)=vxy;
                listm(k,5)=alpha(1)+angstp*xyangl;
                listm(k,6)=vz;
                listm(k,7)=flag1;
                listm(k,8)=flag2;
                k=k+1;
                
                xylmt=abs(max([max(abs(x(:,1))),max(abs(x(:,2)))]));
                
                %plot
                subplot(2,2,1);
                plot3(0,0,0,'B+',x(:,1),x(:,2),x(:,3),'R.');
                xlim([-xylmt xylmt]);
                ylim([-xylmt xylmt]);
                xlabel('X');
                ylabel('Y');
                zlabel('Z');
                grid;
                subplot(2,2,2);
                plot(0,0,'B+',x(:,1),x(:,2),'R.');
                xlim([-xylmt xylmt]);
                ylim([-xylmt xylmt]);
                xlabel('X');
                ylabel('Y');
                grid;
                subplot(2,2,3);
                plot(0,0,'B+',x(:,2),x(:,3),'R.');
                xlim([-xylmt xylmt]);
                xlabel('Y');
                ylabel('Z');
                grid;
                subplot(2,2,4);
                plot(0,0,'B+',x(:,1),x(:,3),'R.');
                xlim([-xylmt xylmt]);
                xlabel('X');
                ylabel('Z');
                grid;
                
                
                %save file
                fln=[paths,num2str(ix),'_',num2str(iy),'_',num2str(iz),'_',num2str(vxy),'_',num2str(alpha(1)+xyangl*angstp),'_',num2str(vz),'_',num2str(flag1),'_',num2str(flag2)];
%                 saveas(gcf,[fln,'.jpg']);
                
                %pause(0.5);
                clear t;
                clear x;
            end %angle
        end %vz
    end %vxy
end %z

完整代码添加QQ1575304183