main-04

# filename: main-04.m
# writer: won sunggyu
# date: 2025-03-28
# language: octave
# description: newmark-beta method (ode) and transfer function

run("../startup.m");
addpath(genpath("../mylib"));
printf(fmt("{mfilename}\n", "#FF5733"));

% 시스템 매개변수 설정
M = [2, 0; 0, 2] * 0.1;  % 질량 행렬 (2 자유도)
C = [0.5, 0; 0, 0.5] * 0.5;  % 감쇠 행렬
K = [200, -100; -100, 200] * 4;  % 강성 행렬

% 시스템 시뮬레이션
tspan = [0, 10];  % 시간 구간 (0부터 10까지)
x0 = [0; 0];  % 초기 위치
v0 = [0; 0];  % 초기 속도
dt = 0.0001;  % 시간 간격
impulse_span = [1.000, 1.020];  % 임펄스 구간 (0.5초부터 1.0초까지 하프 사인)

# 외력 행렬 생성
F = generate_impulse_force(tspan, dt, impulse_span, x0); # 2xN

# 뉴마크-베타 방법
[t1, x1, v1, a1] = newmark_beta(M, C, K, F, tspan, x0, v0, dt); # 2xN

# FFT
duration = tspan(2) - tspan(1);
df = 1 / duration;
nn = floor(duration / dt);
tt = t1; # 1xN
ff = 0: df: df * (nn - 1); # 1xN

FF = two_sided(F); # 2xN
AA = two_sided(a1); # 2xN
AF = transfer_function(a1, F); # 2xN

% 결과플로팅
figured("Size", [960, 960], "Move", [-500, 0], "Name", "Newmark-β");
axes = subplots(3, 1, "Xlabel", "Time [s]");

plotd(axes(1), t1, x1(1, :), ";DOF 1;");
plotd(axes(1), t1, x1(2, :), ";DOF 2;");
plotd(axes(2), t1, v1(1, :), ";DOF 1;");
plotd(axes(2), t1, v1(2, :), ";DOF 2;");
plotd(axes(3), t1, a1(1, :), ";DOF 1;");
plotd(axes(3), t1, a1(2, :), ";DOF 2;");

set(axes(1), "Ylabel", "Position [m]");
set(axes(2), "Ylabel", "Velocity [m/s]");
set(axes(3), "Ylabel", "Acceleration [m/s^2]");

figured("Size", [960, 960], "Move", [500, 0],"Name", "Newmark-β");
axes = subplots(4, 1, "Xlabel", "Frequency [Hz]", "Xlim", [0, 50]);

plotd(axes(1), ff, abs(FF(1, :)), ";FFm;");
plotd(axes(2), ff, abs(AA(1, :)), ";AAm;");
plotd(axes(3), ff, abs(AF(1, :)), ";AFm;");
plotd(axes(4), ff, arg(AF(1, :)), ";AFp;");

set(axes(1), "Ylabel", "Force [N]");
set(axes(2), "Ylabel", "Acceleration [m/s^2]");
set(axes(3), "Ylabel", "AFm [m/s^2/N]", "Ylim", [0, 200]);
set(axes(4), "Ylabel", "AFp [rad]");