                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        clc; clear; close all

m = 0.2533;
k = 10;
zeta = 0.05;

% (1)
wn = sqrt(k/m);
wD = wn*sqrt(1 - zeta^2);
disp(['wn = ', num2str(wn), ' rad/sec'])
disp(['wD = ', num2str(wD), ' rad/sec'])

% (2)
dt = 0.01;
t = [0:dt:1];

p0 = 10;
w = pi/0.6;                                                         
C = p0/k*(1 - (w/wn)^2)/((1 - (w/wn)^2)^2 + (2*zeta*(w/wn))^2);
D = p0/k*(-2*zeta*(w/wn))/((1 - (w/wn)^2)^2 + (2*zeta*(w/wn))^2);
u_steady = C*sin(w*t) + D*cos(w*t);
                 
% (3)
c = 2*m*wn*zeta;
u0 = 0; udot0 = 0; 
p = p0*sin(w*t);

gamma = 1/2; beta = 1/4;
[u1, udot1, udotdot1] = newmarkBeta(m, c, k, p, u0, udot0, dt, gamma, beta);
gamma = 1/2; beta = 1/6;
[u2, udot2, udotdot2] = newmarkBeta(m, c, k, p, u0, udot0, dt, gamma, beta);

figure(1)
plot(t, u_steady, t, u1, t, u2)
legend('稳态解', '平均加速度', '线性加速度')
xlabel('t (s)')
ylabel('u (m)')

% 拓展1 - 计算包含瞬态反应的理论解
A = u0 - D;
B = (udot0 - w*C + zeta*wn*A)/wD;
u_total = exp(-zeta*wn*t).*(A*cos(wD*t) + B*sin(wD*t)) + u_steady;

figure(2)
plot(t, u_total, t, u_steady, t, u1, t, u2)
legend('理论解','稳态解', '平均加速度', '线性加速度')
xlabel('t (s)')
ylabel('u (m)')

% 拓展2 - 计算t = [0, 100]s的结构响应
t = [0:dt:100];
p = p0*sin(w*t);

u_steady = C*sin(w*t) + D*cos(w*t);
u_total = exp(-zeta*wn*t).*(A*cos(wD*t) + B*sin(wD*t)) + u_steady;

gamma = 1/2; beta = 1/4;
[u1, udot1, udotdot1] = newmarkBeta(m, c, k, p, u0, udot0, dt, gamma, beta);
gamma = 1/2; beta = 1/6;
[u2, udot2, udotdot2] = newmarkBeta(m, c, k, p, u0, udot0, dt, gamma, beta);

figure(3)
plot(t, u_total, t, u_steady, t, u1, t, u2)
legend('理论解','稳态解', '平均加速度', '线性加速度')
xlabel('t (s)')
ylabel('u (m)')

