% ビームの動的変形（断面積一定）

L = 10; A = 2;
n = 6; h = L/n;
E = 200; % 弾性率（ヤング率）
c = 10;  % 粘性率
rho = 0.5; % 密度
alpha = 2000; % 制約安定化法
tf = 20.0; % シミュレーションの終端時刻

e0 = (4/6)*ones(n+1,1);
e0(1) = (2/6); e0(n+1) = (2/6);
e1 = (1/6)*ones(n+1,1);
M = (rho*A*h)*spdiags([e1 e0 e1], -1:1, n+1, n+1); % 三重対角行列
fprintf("慣性行列 M\n");
full(M)

e0 = 2*ones(n+1,1);
e0(1) = 1; e0(n+1) = 1;
e1 = (-1)*ones(n+1,1);
K = (E*A/h)*spdiags([e1 e0 e1], -1:1, n+1, n+1); % 三重対角行列
fprintf("剛性行列 K\n");
full(K)

B = (c*A/h)*spdiags([e1 e0 e1], -1:1, n+1, n+1); % 三重対角行列
fprintf("粘性行列 K\n");
full(B)

a = zeros(n+1,1);
a(1) = 1;
fprintf("ベクトル a\n");
a

beam_dynamic_equation = @(t,q) beam_dynamic_equation_param (t,q, n,M,B,K,a,alpha);

interval = [0,tf];
qinit = [ zeros(n+1,1); zeros(n+1,1) ];
[time,q] = ode45(beam_dynamic_equation, interval, qinit);

clf;
plot(time, beam_dynamic_external_force(time));
fprintf("外力\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

clf;
hold on;
for k=1:n+1
    plot(time,q(:,k));
end
hold off;
fprintf("節点の変位\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

clf;
hold on;
plot(time,q(:,n+1));
plot(time,q(:,2*(n+1)));
hold off;
fprintf("右端の変位と速度\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;
