!***********************************************************************
!*   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: "Mécanique des vibrations linéaires By M. Lalanne,       *
!*             P. Berthier, J. Der Hagopian, Masson, Paris 1980"       *
!*             [BIBLI 16].                                             *
!*                                                                     *
!*                                  F90 Release By J-P Moreau, Paris.  *
!*                                         (www.jpmoreau.fr)           *
!***********************************************************************
Program ONEDOF02

INTEGER, PARAMETER :: NMAX=50

REAL M, K, A(NMAX), B(NMAX), X(NMAX), T4(NMAX)

PI = 3.1415926535

PRINT *,' '
WRITE(*,10,advance='no'); READ *, M
WRITE(*,20,advance='no'); READ *, K
WRITE(*,30,advance='no'); READ *, C
PRINT *,' '
PRINT *,'                N'
PRINT *,' F(t) = A0/2 + Sum [A(P) COS(PWt) + B(P) SIN(PWt)]'
PRINT *,'               P=1'
PRINT *,' '
WRITE(*,40,advance='no'); READ *, N
PRINT *,' '
WRITE(*,50,advance='no'); READ *, F8
W = 2 * PI * F8
PRINT *,' '
WRITE(*,60,advance='no'); READ *, A0
PRINT *,' '
DO I = 1, N
  WRITE(*,80,advance='no') I; READ *, A(I)
END DO
PRINT *,' '
DO I = 1, N
  WRITE(*,90,advance='no') I; READ *, B(I)
END DO
PRINT *,' '
PRINT *,' Time Scanning'
PRINT *,' '
T0 = 2. * PI / W
WRITE(*,100)  T0
PRINT *,' '

WRITE(*,110,advance='no'); READ *, T1
WRITE(*,120,advance='no'); READ *, T2
WRITE(*,130,advance='no'); READ *, T3

!N1 = 1 + (T2 - T1) / T3  !Number of response points

J = 0
Y1 = 1E+20
Y2 = 0.

!Main loop
DO T = T1, T2, T3
  J = J + 1
  X(J) = A0 / K / 2.
  T4(J) = T
  DO I = 1, 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.EQ.0.) THEN
      PRINT *,' STOP - ZERO DIVIDE ERROR!'
      STOP
    ELSE
      X(J) = X(J) + A(I) * (C1 * C0 + C2 * S0) / D + B(I) * (C1 * S0 - C2 * C0) / D
    END IF
  END DO
  X3 = X(J)
  IF (Y1 > X3)  Y1 = X3
  IF (Y2 < X3)  Y2 = X3
END DO

J = 0
PRINT *,' '
PRINT *,'   Time     Displacement '
PRINT *,' ------------------------'
DO T = T1, T2, T3
  J = J + 1
  WRITE(*,200) T, X(J)
END DO

STOP

10 FORMAT(' Mass M = ')
20 FORMAT(' Stiffness K = ')
30 FORMAT(' Viscous damping coefficient C = ')
40 FORMAT(' N = ')
50 FORMAT(' Excitation frequency = ')
60 FORMAT(' A0 = ')

80 FORMAT(' A(',I2,')= ')
90 FORMAT(' B(',I2,')= ')

100 FORMAT(' Excitation Period = ',F7.3,' Sec.')
110 FORMAT(' Starting time = ')
120 FORMAT(' Ending time...= ')
130 FORMAT(' Time step.....= ')

200 FORMAT(F9.4,F13.6)

END

!end of file 1dof02.f90