/**********************************************************************
*   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].                                             *
*                                                                     *
*                                  C++ Release By J-P Moreau, Paris.  *
*                                  (www.jpmoreau.fr)                  *
**********************************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  25
#define  PI    4*atan(1.0)

    double A0,A1,C,C0,C1,C2,D,F8,M,K,S0,T,T0,T1,T2,T3,TMP,W,X3,Y1,Y2;
    double A[NMAX],B[NMAX],T4[NMAX],X[NMAX];   //A(i),B(i), coefficients of periodic force
    int I,J,N,N1;


void main() {

  printf("\n");
  printf(" Mass M = "); scanf("%lf", &M);
  printf(" Stiffness K = "); scanf("%lf", &K);
  printf(" Viscous damping coefficient: "); scanf("%lf", &C);

  printf("\n");
  printf("                N\n");
  printf(" F(t) = A0/2 + Sum [A(P) COS(PWt) + B(P) SIN(PWt)]\n");
  printf("               P=1\n");
  printf("\n");
  printf(" N = "); scanf("%d", &N);
  printf("\n");
  printf(" Excitation frequency = "); scanf("%lf", &F8);
  W = 2.0 * PI * F8;
  printf("\n");
  printf(" A0 = "); scanf("%lf", &A0);
  printf("\n");
  for (I=1; I<=N; I++) {
    printf(" A(%d)= ", I); scanf("%lf", &TMP);
    A[I]=TMP;
  }
  printf("\n");
  for (I=1; I<=N; I++) {
    printf(" B(%d)= ", I); scanf("%lf", &TMP);
    B[I]=TMP;
  }
  printf("\n Time Scanning\n\n");
  T0 = 2.0 * PI / W;
  printf(" Excitation Period = %7.3f\n\n", T0);
  printf(" Starting time = "); scanf("%lf", &T1);
  printf(" Ending time...= "); scanf("%lf", &T2);
  printf(" Time step.....= "); scanf("%lf", &T3);

  N1 = (int) floor(1 + (T2 - T1) / T3);  // Number of response points
  J = 0;
  Y1 = 1e+20;
  Y2 = 0.0;

  // Main loop
  T=T1;
  do {
    J++;
    X[J] = A0 / K / 2.0;
    T4[J] = T;
    for (I=1; I<=N; I++) {
      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.0) {
        printf("STOP - ZERO DIVIDE ERROR!\n\n");
        return;
      }
      else
        X[J] = X[J] + A[I] * (C1 * C0 + C2 * S0) / D + B[I] * (C1 * S0 - C2 * C0) / D;
    }
    X3 = X[J];
    if (Y1 > X3)  Y1 = X3;
    if (Y2 < X3)  Y2 = X3;
    T += T3;
  } while (T <= T2+T3);

  J = 0;
  printf("\n");
  printf("   Time     Displacement \n");
  printf(" ------------------------\n");
  T=T1;
  do {
    J++;
    printf("%8.4f  %11.6f\n", T, X[J]);
    T += T3;
  } while (T <= T2+T3);

  printf("\n");

}

// end of file 1dof02.cpp