MATLAB模拟伽尔顿板实验

MATLAB模拟伽尔顿板实验代码展示：

clear;

b=1;

k=1;

dxdy=0.81;

r=0.1;ttt=0:pi/300:2*pi;xxx=r*exp(i*ttt);

for n=1:80

    for m=10:20

if xor(mod(m,2),mod(n,2))

            xy0(1,k)=dxdy*n;

            xy0(2,k)=dxdy*m;

            ax=xy0(1,k);

            ay=xy0(2,k);

            plot(xxx+-1+i*(-2),'r.'), plot(ax,ay,'r.','markersize',15) ,hold on

            k=k+1;

        end

    end

end

k0=k-1;

axis ([-3 n*dxdy+3 -2 m*dxdy+5])

ll=40/80*dxdy;

L=3000;

r1=n/2*dxdy*ones(1,L)+[1.5*dxdy*rand(1,L/2) -1.5*dxdy*rand(1,L/2)];

r2=(m*dxdy+1)+1.5*dxdy*rand(1,L);

xyxy=[r1;r2];

xy=[r1;r2];

vy=-1*rand(1,L);

vx=[2*rand(1,L/2) -2*rand(1,L/2)];

v=[vx;vy];

ball=plot(r1(1,:),r2(1,:),'.','markersize',15,'color','b')

hold off

pause

xy=xyxy;

m=1;

for k=1:3000

   dt=0.015;

    g=59.8;

vxm(k)=mean(v(1,:));vym(k)=mean(v(2,:));

vxstd(k)=std(v(1,:));vystd(k)=std(v(2,:));

    for numball=1:L

           if xy(2,numball)<=5

                 xy(2,numball)=5;

                 v(2,numball)=0;

                 v(1,numball)=0;

           end

if xy(1,numball)<dxdy | xy(1,numball)>=n*dxdy

   xy(1,numball)=(xy(1,numball)<dxdy)*(2*dxdy-xy(1,numball))+(xy(1,numball))>n*dxdy)*(2*n*dxdy-xy(1,

numball));

                 v(1,numball)=-v(1,numball);

             end

           dtt=dt/3*2;

           dttt=dt/3;

          dtttt=dt/2;

            xy(1,numball)=xy(1,numball)+dt*v(1,numball);

            xy(2,numball)=xy(2,numball)+dt*v(2,numball)-0.5*g*dt^2;

            v(2,numball)=v(2,numball)-g*dt;

            if xy(2,numball)>10*dxdy  

            al=xy(1,numball)-dtt*v(1,numball);

            all=xy(1,numball)-dttt*v(1,numball);

            alll=xy(1,numball)-dtttt*v(1,numball);

            bl=xy(2,numball)-dtt*v(2,numball)-0.5*g*dtt^2;

            bll=xy(2,numball)-dttt*v(2,numball)-0.5*g*dttt^2;

            blll=xy(2,numball)-dtttt*v(2,numball)-0.5*g*dtttt^2;

            cl=[al;bl];

            cll=[all;bll];

            clll=[alll;blll];

             for numxy=1:k0

              dxy=norm(xy(:,numball)-xy0(:,numxy));

              dxyl=norm(cl-xy0(:,numxy));

              dxyll=norm(cll-xy0(:,numxy));

              dxylll=norm(clll-xy0(:,numxy));

              if dxy<=0.38 | dxyl<=0.38 | dxyll<=0.38 | dxylll<=0.38

                 xy(:,numball)=xyxy(:,numball);

                 v(2,numball)=v(2,numball)+g*dt;    

                 a=xy(1,numball);b=v(1,numball);

                 c=xy(2,numball);d=v(2,numball);

                 e=xy0(1,numxy);f=xy0(2,numxy);

                 t = fsolve(@myfun2,0.000,optimset('Display','off'),a,b,c,d,e,f,g);

                 t0=min(t);

                 xy(1,numball)=xy(1,numball)+t0*v(1,numball);

                 xy(2,numball)=xy(2,numball)+t0*v(2,numball)-0.5*g*t0^2;

                  %set(ball,'xdata',xy(1,:),'ydata',xy(2,:));

                 v(2,numball)=v(2,numball)+g*t0;

                 y1=xy(2,numball);y2=xy0(2,numxy);

                 x1=xy(1,numball);x2=xy0(1,numxy);

                 vx=v(1,numball);

                 vy=v(2,numball);

                 axy=atan((y2-y1)/(x2-x1+eps));

                 v(1,numball)=(-vx*cos(2*axy)-vy*sin(2*axy));

                 v(2,numball)=(vy*cos(2*axy)-vx*sin(2*axy));

                 t0=dt-t0;

                 xy(1,numball)=xy(1,numball)+t0*v(1,numball);

                 xy(2,numball)=xy(2,numball)+t0*v(2,numball)-0.5*g*t0^2;

                 v(2,numball)=v(2,numball)-g*t0;

              end

           end 

    end

end

    xyxy=xy;

set(ball,'xdata',xy(1,:),'ydata',xy(2,:));

drawnow;

end

figure(1)

histfit(xy(1,:),31)

figure(2); 

   m=1:k;

   kk=m*dt;

plot(kk,vxm,'r:',kk,vym,'b-')

  figure(3); 

plot(kk,vxstd,'r:',kk,vystd,'b-')

figure(4)

vx1=diff(vxm);vy1=diff(vym);

n=1:(k-1);

kkt=n*dt;

plot(kkt,vx1,'r:',kkt,vy1,'b-')