damping.m

用matlab编了一个damping程序

function damping(r,R,h,l,m)
a=0;
q=acos((r-h)/r);
t=1;
w=3*pi/2;p=0;
k=((2*q)*sqrt(9.8))/((5*pi)*sqrt(r));u=1;
d=0; x=sqrt(((r-(r*cos(q)))*(r-(r*cos(q))))+(r*sin(q)*(r*sin(q))));x;u=1;
while p<100
    while a<=4*q
    while a<q
        axes('position',[0 .1 1 .6]);
        circle(r*cos(w-a),r*sin(w-a),R);
        line(0,0,r*cos(w-a),r*sin(w-a));
        axis([-20 20 -10 10]);plot(0,0,'vr');
        s=0:(2*pi*sqrt(r/9.8))/20:(((a/4*q)*pi*sqrt(r/9.8))+d);
        i=size(s);
        while u<=i(1,2)
            b(1,u)=x*sin(sqrt(9.8/r)*s(1,u))*(exp(-(l/(2*m))*s(1,u)));
            u=u+1;
        end
        axes('position',[.1 .6 .8 .25]);plot(s,b);
        M(t)=getframe;clf;
        t=t+1;m=t;
        a=a+k;
    end
    q=q*exp(-(l/(2*m))*(pi/2)*sqrt(r/9.8));u=1;b=0;
    while q<=a && a<=3*q
        axes('position',[0 .1 1 .6]);
        circle(r*cos(w-2*q+a),r*sin(w-2*q+a),R);
        line(0,0,r*cos(w-2*q+a),r*sin(w-2*q+a));
        axis([-20 20 -10 10]); plot(0,0,'vr');
        s=0:(2*pi*sqrt(r/9.8))/20:(((a/4*q)*pi*sqrt(r/9.8))+d);
         i=size(s);
        while u<=i(1,2)
            b(1,u)=x*sin(sqrt(9.8/r)*s(1,u))*(exp(-(l/(2*m))*s(1,u)));
            u=u+1;
        end
        axes('position',[.1 .6 .8 .25]);plot(s,b);
        M(t)=getframe;
        clf
        t=t+1;m=t;
        a=a+k;
    end
    q=q*exp(-(l/(2*m))*(pi/2)*sqrt(r/9.8));u=1;b=0;
    while 3*q<a && a<=4*q
        axes('position',[0 .1 1 .6]);
        circle(r*cos(w+4*q-a),r*sin(w+4*q-a),R);
        line(0,0,r*cos(w+4*q-a),r*sin(w+4*q-a));
        axis([-20 20 -10 10]); plot(0,0,'vr');
        s=0:(2*pi*sqrt(r/9.8))/20:(((a/4*q)*pi*sqrt(r/9.8))+d);
         i=size(s);
        while u<=i(1,2)
            b(1,u)=x*sin(sqrt(9.8/r)*s(1,u))*(exp(-(l/(2*m))*s(1,u)));
            u=u+1;
        end
        axes('position',[.1 .6 .8 .25]);plot(s,b);
        M(t)=getframe;
        clf
        t=t+1;m=t;
        a=a+k;
    end
    q=q*exp(-(l/(2*m))*(pi/2)*sqrt(r/9.8));u=1;b=0;
    end
    a=0;d=d+(2*pi*sqrt(r/9.8));
    t=m+1;
    p=p+1;
end
end