Gravity Shape

gravityShape_v2.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   -- gravityShape_v2.m --
%
%	引力曲線の可視化
%   オブジェクト数を拡張
%	野原友幸 2014.10.02
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% $$m \frac{d \vec{v}}{dt}= \vec{F};$$
%
% $$\vec{F}= K \sum_{i}^{}{m\frac{M_g,i}{|\vec{r_i}|} };$$
%
% $$\frac{\vec{v}^{t}-\vec{v}^{t-1}}{ \bigtriangleup t } =K
% \sum_{i}^{}{\frac{M_g,i}{|\vec{r_i}|} };$$
%
% $$\vec{v}^{t} = \vec{v}^{t-1}+K \sum_{i}^{}{\frac{M_g,i}{|\vec{r_i}|} };$$
%
%%
% 
% $$e^{\pi i} + 1 = 0$$
% 

break
clear all; close all; clc;

%% ワークフロー
%% 変数の定義
dt=0.0001;
K=1;
Mg_1=100;
Mg_2=100;
Ma = 1;

%% 数値計算設定値
MAX_STEP = 300000;
INI_STEP=10;

%% 初期化
Xg_1=0;
Yg_1=0;

Xg_2=0;
Yg_2=-10;

X = zeros(100,1);
Y = zeros(100,1);
VX = zeros(100,1);
VY = zeros(100,1);

X(1:INI_STEP) = 5;
Y(1:INI_STEP) = 0;
VX(1:INI_STEP) = 10;
VY(1:INI_STEP) = 10;


%% 時間発展
figure;
for step = INI_STEP:MAX_STEP
    
%% アクターと質量との距離の計算
Dx = X(step-1) - Xg_1;
Dy = Y(step-1) - Yg_1;
D1 = sqrt( Dx^2 + Dy^2 );
cos1=Dx/D1;
sin1=Dy/D1;
Dx = X(step-1) - Xg_2;
Dy = Y(step-1) - Yg_2;
D2 = sqrt( Dx^2 + Dy^2 );
cos2=Dx/D2;
sin2=Dy/D2;

%% 引力の更新
FX(step) = -1*dt*K * ((Mg_1/D1)* cos1 + (Mg_2/D2)* cos2);
FY(step) = -1*dt*K * ((Mg_1/D1)* sin1 + (Mg_2/D2)* sin2);

% nu=0.00001;
% FX(step) = -1*dt*K * ((Mg_1/D1)* cos1 + (Mg_2/D2)* cos2 )+ VY(step-1)*nu;
% FY(step) = -1*dt*K * ((Mg_1/D1)* sin1 + (Mg_2/D2)* sin2) + VX(step-1)*nu;

%% 速度の更新
VX(step) = VX(step-1) + FX(step);
VY(step) = VY(step-1) + FY(step);

%% 座標位置の更新
X(step)=X(step-1)+dt*VX(step);
Y(step)=Y(step-1)+dt*VY(step);

%% アニメーションの表示
% plot(X(step-2:step),Y(step-2:step), 'ob');
% hold on;
% plot(Xg_1, Yg_1, 'or');
% hold off;
% grid on;
% axis([-10 10 -10 10])
% 
% drawnow

end


%% リサージュ plot
ra=2;
figure('Position',[20 50 320*ra 320*ra]);
plot(VX,VY, '-b');
grid on;
axis([min(VX) max(VX) min(VY) max(VY)]);
xlabel('Vx');
ylabel('Vy');
title('引力曲線シミュレーション');

%% 軌跡 plot
ra=2;
figure('Position',[ra*320+20 50 320*ra 320*ra]);
plot(X,Y, '-b');
hold on;
plot(X(INI_STEP),Y(INI_STEP),'ob');
plot(X(end),Y(end),'ob');
plot(Xg_1, Yg_1, 'or');
plot(Xg_2, Yg_2, 'or');
hold off;
grid on;
axis([min(X) max(X) min(Y) max(Y)]);
xlabel('X');
ylabel('Y');
title('引力曲線シミュレーション');

%% アニメーションの表示
figure;
SKIP = 100;
TAILN = 20;
for step = INI_STEP+SKIP*TAILN:SKIP:MAX_STEP
  
    plot(X(step),Y(step), 'ob', 'MarkerSize',5);
    hold on;
    plot(X(step-SKIP*TAILN:step), Y(step-SKIP*TAILN:step) ,'-b');
    plot(Xg_1, Yg_1, '.r', 'MarkerSize',5*Mg_1/Ma/10);
    plot(Xg_2, Yg_2, '.r', 'MarkerSize',5*Mg_2/Ma/10);
    hold off;
    grid on;
    axis([min(X) max(X) min(Y) max(Y)]);
    xlabel('X');
    ylabel('Y');
    title('引力曲線シミュレーション');

	drawnow;
end

break



break
%% figure open
ra=2;
% figure('Position',[20 50 640*ra 320*ra]);
figure;
m=6; n=1;g=1;
% 変位の表示
subplot(m,n,g);g=g+1; plot(X);  title('X');
grid on;
% 変位の表示
subplot(m,n,g);g=g+1; plot(Y);  title('Y');
grid on;
% 変位の表示
subplot(m,n,g);g=g+1; plot(VX);  title('VX');
grid on;
% 変位の表示
subplot(m,n,g);g=g+1; plot(VY);  title('VY');
grid on;
% 変位の表示
subplot(m,n,g);g=g+1; plot(FX,'r'); title('FX');
grid on;
% FYの表示
subplot(m,n,g);g=g+1; plot(FY,'r'); title('FY');
grid on;
% Fの表示
subplot(m,n,g);g=g+1; plot(F,'r'); title('Y');
grid on;