/**********************************************************************
* Response of a 1dof Mass-Spring System with damping to a sinusoidal  *
* input force                                                         *
* ------------------------------------------------------------------- *
* Main Variables:                                                     *
*                  M: Mass                                            *
*                  K: Stiffness                                       *
*                  C: Viscous Damping Coefficient                     *
*                  E: Structural Damping Coefficient                  *
*                 F0: Input Force                                     *
* ------------------------------------------------------------------- *
* SAMPLE RUN:                                                         *
*                                                                     *
*  Mass M..... = 50                                                   *
*  Stiffness K = 1e6                                                  *  
*  Force F0... = 10000                                                *
*                                                                     *
*  (V)iscous or (S)tructural Damping: S                               *
*                                                                     *
*  Structural damping factor = 0.05                                   *
*                                                                     *
*  Frequency Scanning                                                 *
*                                                                     *
*  Resonance Frequency of undamped system = 22.508                    *
*                                                                     *
*  Starting frequency = 20                                            *
*  Ending frequency...= 25                                            *
*  Frequency step.....= 0.25                                          *
*                                                                     *
* Frequency   Displacement   Phase (deg)                              *
* --------------------------------------                              * 
*  20.00      0.046234       13.366013                                *
*  20.25      0.050756       14.701446                                *
*  20.50      0.056293       16.347699                                *
*  20.75      0.063206       18.423074                                *
*  21.00      0.072037       21.111461                                *
*  21.25      0.083609       24.711299                                *
*  21.50      0.099181       29.729354                                *
*  21.75      0.120525       37.058260                                *
*  22.00      0.149218       48.252814                                *
*  22.25      0.181993       65.500492                                *
*  22.50      0.199980       89.194985                                *
*  22.75      0.183564      -66.609176                                *
*  23.00      0.149839      -48.520651                                *
*  23.25      0.119585      -36.721514                                *
*  23.50      0.097048      -29.028193                                *
*  23.75      0.080680      -23.790753                                *
*  24.00      0.068578      -20.053143                                *
*  24.25      0.059388      -17.273961                                *
*  24.50      0.052222      -15.136029                                *
*  24.75      0.046502      -13.444959                                *
*  25.00      0.041843      -12.076311                                *
*                                                                     *
* ------------------------------------------------------------------- *
* REFERENCE: "Mécanique des vibrations linéaires By M. Lalanne,       *
*             P. Berthier, J. Der Hagopian, Masson, Paris 1980".      *
*                                                                     *
*                                  C++ Release By J-P Moreau, Paris.  *
*                                  (www.jpmoreau.fr)                  *
**********************************************************************/
#include <stdio.h>
#include <math.h>


    double C,D,E,F,F0,F1,F2,F3,F4,M, K, I1,M1,O1,PI,T1,T3,TMP,W;
    int C2;
    char B[2];


void main() {

  PI = 4.0*atan(1.0);

  printf("\n");
  printf(" Mass M..... = "); scanf("%lf", &M);
  printf(" Stiffness K = "); scanf("%lf", &K);
  printf(" Force F0... = "); scanf("%lf", &F0);
  printf("\n");

  do {
    printf(" (V)iscous or (S)tructural Damping: "); scanf("%s", B);
  } while (B[0] != 'V' && B[0] != 'S');

  printf("\n");
  if (B=="V") {
    printf(" Viscous damping coefficient C = "); scanf("%lf", &C);
    C2 = 1;
  }
  else {
    printf(" Structural damping factor = "); scanf("%lf", &E);
    C2 = 2;
  }

  printf("\n Frequency Scanning\n\n");
  F4 = 1.0 / 2.0 / PI * sqrt(K / M);
  printf(" Resonance Frequency of undamped system = %7.3f Hz.\n\n", F4);
  printf(" Starting frequency = "); scanf("%lf", &F1);
  printf(" Ending frequency...= "); scanf("%lf", &F2);
  printf(" Frequency step.....= "); scanf("%lf", &F3);

  F = F1;

  printf(" Frequency   Displacement   Phase (deg) \n");
  printf(" -------------------------------------- \n");

  while (F<=F2) {
    W = 2.0 * PI * F;
    TMP=K-M*W*W;
    if (C2 == 1) {  //Viscous damping
      D  = TMP*TMP + C*W*C*W;
      O1 = TMP / sqrt(D);
      I1 = C * W / sqrt(D);
      M1 = F0 / sqrt(D);
    }
    else {  //Structural damping
      D  = TMP*TMP + E*K*E*K;
      O1 = TMP / sqrt(D);
      I1 = E * K / sqrt(D);
      M1 = F0 / sqrt(D);
    }

    if (M1 == 0.0)
      T1 = 0.0;
    else {
      if (O1 == -1.0)
        T1 = PI;
      else if (O1 == 1.0)
        T1 = 0.0;
      else {
        if (I1 > 0.0)
          //T1 = ACOS(O1)
          T1 = atan(sqrt(1.0 - O1*O1) / O1);
        else
          T1 = 2.0 * PI - atan(sqrt(1.0 - O1*O1) / O1);
      }
    }

    T3 = T1 / PI * 180;

    printf("%7.2f%14.6f  %14.6f\n", F, M1, T3); 
    F += F3;
  }

}

// end of file 1dof01.cpp