How to plot by using result data?
    2 vues (au cours des 30 derniers jours)
  
       Afficher commentaires plus anciens
    
Dear Sir or Madam,
              I am new to matlab programming I would like to know how to plot position_x , position_y , velocity x , velocity y by applying result  data. Please show me how to perform the plot as I would like to see them. 
Hopefully and Thank you very much for your help.
n_part=10;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
timestep=50;
v_init=0.5;
% initial of position, velocity, angular velocity
position_x(:,1:n_part)=rand(size(1:n_part));
position_y(:,1:n_part)=rand(size(1:n_part));
theta=2*pi*rand(size(position_x));
velocity_x=v_init*sin(theta);
velocity_y=v_init*cos(theta);
w_x=zeros(size(position_x));
w_y=zeros(size(position_y));
w_half_x=zeros(size(velocity_x));
w_half_y=zeros(size(velocity_y));
w_x_dot=zeros(size(velocity_x));
w_y_dot=zeros(size(velocity_y));
v_half=rand(size(velocity_x));
v_half_x=rand(size(velocity_x));
v_half_y=rand(size(velocity_y));
dt=0.001;
rad(1:n_part)=0.5;
acc_x=zeros(size(position_x));
acc_y=zeros(size(position_y));
acc_i=[0,0];
acc_j=[0,0];
w_i=[0,0];
w_j=[0,0];
Fn=[0,0];Fn_x=[0,0];Fn_y=[0,0];
for t=1:1:4
    for i=1:n_part
        v_half_x(t+1,i)=velocity_x(t,i)+0.5*dt*acc_x(t,i);
        v_half_y(t+1,i)=velocity_y(t,i)+0.5*dt*acc_y(t,i);
        position_x(t+1,i)=position_x(t,i)+v_half_x(t+1,i)*dt;
        position_y(t+1,i)=position_y(t,i)+v_half_y(t+1,i)*dt;
        % update of w_ part
        w_half_x(t+1,i)=w_x(t,i)+dt*0.5*dt*w_x_dot(t,i);
        w_half_y(t+1,i)=w_y(t,i)+dt*0.5*dt*w_y_dot(t,i);
        j=1;
        for j=i+1
            if j<=n_part
            lx=position_x(i)-position_x(j);
            ly=position_y(i)-position_y(j);
            vx=velocity_x(i)-velocity_x(j);
            vy=velocity_y(i)-velocity_y(j);
            root_xy=sqrt(lx.^2+ly.^2);
            R_i=rad(i);
            R_j=rad(j);
            delta=(R_i+R_j)-root_xy;
            if delta>0
                N=[lx;ly];
                n_vector=N/norm(N);
                % relative vector to find tangent vector;
                vij1=[vx;vy];
                w_i=[w_x(i);w_y(i)];
                w_j=[w_x(j);w_y(j)];
                vij2=vij1+(R_i.*w_i+R_j.*w_j).*n_vector;
                v_i_j=-n_vector.*(n_vector.*vij2);
                % tangent vector
                tan=v_i_j/norm(v_i_j);
                Fn=kn.*delta;
                Ft=kt.*v_i_j;
                F_total=Fn.*n_vector+Ft.*tan+m*g;
                Fn_x=Fn_x-F_total(1,:);
                Fn_y=Fn_y+F_total(2,:);
                %Momemtum of particles
                M_i=R_i.*n_vector.*Ft;
                M_j=R_j.*n_vector.*Ft;
                % same radius for movement of inertia of circle
                I=pi*R_i/4;
                w_i=0;
                w_j=0;
                w_i=w_i+M_i'/I;
                w_j=w_i+M_j'/I;
                w_x_dot(t,i)=w_x_dot(t,i)+w_i(1);
                w_x_dot(t,i)=w_x_dot(t,i)+w_j(1);
                w_y_dot(t,i)=w_y_dot(t,i)+w_i(2);
                w_y_dot(t,i)=w_y_dot(t,i)+w_j(2);
                w_x_dot(t+1,:)=w_x_dot(t,:);
                w_y_dot(t+1,:)=w_y_dot(t,:);
                acc_i=0;
                acc_i=acc_i+Fn_x./m;
                acc_j=0;
                acc_j=acc_j+Fn_y./m;
                acc_x(t,i)=acc_x(t,i)+acc_i(1);
                acc_y(t,i)=acc_y(t,i)+acc_i(2);
                acc_x(t,j)=acc_x(t,j)+acc_j(1);
                acc_y(t,j)=acc_y(t,j)+acc_j(2);
                acc_x(t+1,:)=acc_x(t,:);
                acc_y(t+1,:)=acc_y(t,:);
            end
            end
        end
        velocity_x(t+1,i)=v_half_x(t+1,i)+0.5*dt*acc_x(t+1,i);
        velocity_y(t+1,i)=v_half_y(t+1,i)+0.5*dt*acc_y(t+1,i);
        w_x(t+1,i)=w_half_x(t+1,i)+0.5*dt*w_x_dot(t+1,i);
        w_y(t+1,i)=w_half_y(t+1,i)+0.5*dt*w_y_dot(t+1,i);
    end
end
% my intension is using rad, position_x, position_y to create circle as
% particles and using velocity_ x and velocity _y to move them. 
% position_x (t,i); position_y (t,i); velocity_x (t,i); velocity_y (t,i);
% rad=0.5;
0 commentaires
Réponses (1)
Voir également
Catégories
				En savoir plus sur Assembly dans Help Center et File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

