'***********************************************************************
'*   Response of a 1dof Mass-Spring System with viscous damping to a   *
'*   periodic input force:  M X" + C X' + K X = F(t)                   *
'* ------------------------------------------------------------------- *
'* Main Variables:                                                     *
'*                  M: Mass                                            *
'*                  K: Stiffness                                       *
'*                  C: Viscous Damping Coefficient                     *
'*                  F(t): Periodic Force of Frequency F8               *
'*                                                                     *
'*                         N                                           *
'* Note:    F(t) = A0/2 + Sum [A(P) COS(PWt) + B(P) SIN(PWt)]          *
'*                        P=1                                          *
'* ------------------------------------------------------------------- *
'* SAMPLE RUN:                                                         *
'*                                                                     *
'*  Mass M = 50                                                        *
'*  Stiffness K = 1e6                                                  *
'*  Viscous damping coefficient: 0.1                                   *
'*                                                                     *
'*                   N                                                 *
'*    F(t) = A0/2 + Sum [A(P) COS(PWt) + B(P) SIN(PWt)]                *
'*                  P=1                                                *
'*                                                                     *
'*  N = 1                                                              *
'*                                                                     *
'*  Excitation frequency = 20                                          *
'*                                                                     *
'*  A0 = 0                                                             *
'*                                                                     *
'*  A(1) = 0                                                           *
'*                                                                     *
'*  B(1) = 10000                                                       *
'*                                                                     *
'*  Time Scanning                                                      *
'*                                                                     *
'*  Excitation Period = 0.050                                          *
'*                                                                     *
'*  Starting time = 0                                                  *
'*  Ending time...= 0.05                                               *
'*  Time step.....= 0.0025                                             *
'*                                                                     *
'*   Time     Displacement                                             *
'* ------------------------                                            *
'*  0.0000    -0.000003                                                *
'*  0.0025     0.014682                                                *
'*  0.0050     0.027930                                                *
'*  0.0075     0.038444                                                *
'*  0.0100     0.045195                                                *
'*  0.0125     0.047521                                                *
'*  0.0150     0.045196                                                *
'*  0.0175     0.038447                                                *
'*  0.0200     0.027935                                                *
'*  0.0225     0.014688                                                *
'*  0.0250     0.000003                                                *
'*  0.0275    -0.014682                                                *
'*  0.0300    -0.027930                                                *
'*  0.0325    -0.038444                                                *
'*  0.0350    -0.045195                                                *
'*  0.0375    -0.047521                                                *
'*  0.0400    -0.045196                                                *
'*  0.0425    -0.038447                                                *
'*  0.0450    -0.027935                                                *
'*  0.0475    -0.014688                                                *
'*  0.0500    -0.000003                                                *
'* ------------------------------------------------------------------- *
'* REFERENCE: "Mecanique des vibrations lineaires By M. Lalanne,       *
'*             P. Berthier, J. Der Hagopian, Masson, Paris 1980" [16]. *
'*                                                                     *
'*                          Quick Basic Release By J-P Moreau, Paris.  *
'***********************************************************************
'Program 1DOF-02

DEFDBL A-H, O-Z
DEFINT I-N

DIM M AS DOUBLE
DIM K AS DOUBLE

PI = 3.1415926535#

CLS
PRINT
INPUT " Mass M = ", M
INPUT " Stiffness K = ", K
INPUT " Viscous damping coefficient: ", C

PRINT : PRINT
PRINT "                N"
PRINT " F(t) = A0/2 + Sum [A(P) COS(PWt) + B(P) SIN(PWt)]"
PRINT "               P=1"
PRINT
INPUT " N = ", N
DIM A(N), B(N)  'coefficients of periodic force
PRINT
INPUT " Excitation frequency = ", F8
W = 2 * PI * F8
PRINT
INPUT " A0 = ", A0
PRINT
FOR I = 1 TO N
  PRINT " A("; I; ")= "; : INPUT "", A(I)
NEXT I
PRINT
FOR I = 1 TO N
  PRINT " B("; I; ")= "; : INPUT "", B(I)
NEXT I
PRINT
PRINT " Time Scanning"
PRINT
T0 = 2# * PI / W
PRINT USING " Excitation Period = ###.### Sec."; T0
PRINT
INPUT " Starting time = ", T1
INPUT " Ending time...= ", T2
INPUT " Time step.....= ", T3

N1 = 1 + (T2 - T1) / T3'Number of response points
DIM X(N1), T4(N1)
J = 0
Y1 = 1E+20
Y2 = 0#

'Main loop
FOR T = T1 TO T2 STEP T3
  J = J + 1
  X(J) = A0 / K / 2#
  T4(J) = T
  FOR I = 1 TO N
    A1 = I * W * T
    C0 = COS(A1)
    S0 = SIN(A1)
    C1 = K - M * I * I * W * W
    C2 = C * I * W
    D = C1 * C1 + C2 * C2
    IF D = 0# THEN
      PRINT "STOP - ZERO DIVIDE ERROR!"
      END
    ELSE
      X(J) = X(J) + A(I) * (C1 * C0 + C2 * S0) / D + B(I) * (C1 * S0 - C2 * C0) / D
    END IF
  NEXT I
  X3 = X(J)
  IF Y1 > X3 THEN Y1 = X3
  IF Y2 < X3 THEN Y2 = X3
NEXT T

J = 0
CLS
PRINT
PRINT "   Time     Displacement "
PRINT " ------------------------"
FOR T = T1 TO T2 + T3 STEP T3
  J = J + 1
  PRINT USING " ###.####"; T;
  PRINT USING " ####.######"; X(J)
NEXT T

END

'end of file 1dof02.bas