/**************************************************************
*   Program to demonstrate the Acceleration Shock Spectrum    *
*                                                             *
*                         C++ version by J-P Moreau, Paris    *
*                                (www.jpmoreau.fr)            *
* ----------------------------------------------------------- *
* SAMPLE RUN:                                                 *
*                                                             *
* Input data file (tfft.dat):                                 *
*                                                             *
* (The test signal contains 3 frequencies: 50, 250, 500 hz)   *
*                                                             *
* 1024                                                        *
* 0.00000000000000E+0000   0.00000000000000E+0000             *
* 1.95503421309917E-0004   3.53914399999776E+0001             *
* 3.91006842619834E-0004   5.95684899999760E+0001             *
* 5.86510263929974E-0004   6.54621699999552E+0001             *
* 7.82013685239669E-0004   5.24038399999845E+0001             *
* ...                      ...                                *
* 1.98826979472187E-0001   2.77372500000183E-0001             *
* 1.99022482893497E-0001  -2.43361500000174E+0000             *
* 1.99217986314807E-0001  -4.84236799999780E+0000             *
* 1.99413489736116E-0001  -6.02247899999929E+0000             *
* 1.99608993157426E-0001  -5.45615399999951E+0000             *
* 1.99804496578736E-0001  -3.22824200000105E+0000             *
* 2.00000000000045E-0001  -2.96010699999982E-0003             *
*                                                             *
* Output file (tshocksp.lst):                                 *
*                                                             *
*     Frequency     Shock Spectrum      Shock Spectrum        *
*       (Hz)           negative            positive           *
* -------------------------------------------------------     *
*        10.00        -3.037656           3.602233            *
*        20.00        -6.422721           9.166994            *
*        30.00       -16.325201          13.674868            *
*        40.00       -27.368246          27.100563            *
*        50.00       -41.964013          42.437765            *
*        60.00       -27.828617          26.667987            *
*        ...         ...                 ...                  *
*       740.00      -109.784192         112.831489            *
*       750.00      -109.242479         110.265918            *
*       760.00      -108.338595         108.819857            *
*       770.00      -107.086906         108.976530            *
*       780.00      -105.503821         110.103584            *
*       790.00      -103.607643         110.830962            *
*       800.00      -101.418411         109.903625            *
*                                                             *
* Note: the shock spectrum shows 3 peaks around 50, 250 and   *
*       500 hz.                                               *
**************************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  1025

typedef double Vect[SIZE];
typedef double Spectrum[101];

FILE   *f_in, *f_out;
int    i,j,ndata,nfreq;
double df,dt,dzeta,f,fbegin,fend,temp,tbegin,tend,y;
double mini,maxi;
Vect   signal,response;
Spectrum SP1,SP2;    


  /***********************************************************
  * This procedure calculates the acceleration of the sismic *
  * mass of an elementary oscillator (1 degree of freedom),  *
  * the basis of which is submitted to a given acceleration, *
  * x"(t).                                                   *
  *  The input signal x" is digitalized with a constant time *
  * step, Sampling_Incr. The table signal contains N accele- *
  * ration values.                                           *
  * -------------------------------------------------------- *
  * INPUTS:                                                  *
  *         Frequency....: eigen frequency of oscillator     *
  *         signal.......: input acceleration signal x"(t)   *
  *         Sampling_Incr: time step of signal (constant)    *
  *         Dzeta........: reduced damping factor            *
  *         N............: number of points of signal        * 
  * OUTPUT:                                                  *        
  *         x_y..........: table containing acceleration     *
  *                        response of oscillator (N points) *
  *                                                          *
  *            From Pascal version by J-P Moreau/J-P Dumont  *
  ***********************************************************/
  void _1dof_Oscillator_Response
  ( double frequency, Vect signal, Vect x_y,
    double Sampling_Incr,double Dzeta,int N  )   {
  
    double omega,Q,dQ2,delta,
    p0,p1,Pi,q1,q2,sq2,arg,cosA,
    yn,ynm1,ynm2,xn,xnm1;
    int i;
  
    Pi=3.1415926535;
    omega=2*Pi*frequency;
    if (dzeta<1E-6)  dzeta=1E-6;
    Q=1.0/(2.0*Dzeta);
    dQ2=2*Q*Q;
    delta=sqrt(2.0*dQ2-1.0);
    p0=omega*Sampling_Incr/Q;
    q2=exp(-p0);
    sq2=sqrt(q2);
    arg=0.5*p0*delta;
    cosA=cos(arg);
    q1=-2.0*sq2*cosA;
    p1=p0*sq2*((dQ2-1.0)*sin(arg)/delta-cosA);
    ynm1=0.0;
    xn=signal[1];
    yn=0.0;
    x_y[1]=yn;
    for (i=2; i<N+1; i++) {
      ynm2=ynm1;ynm1=yn;xnm1=xn;
      xn=signal[i];
      yn=p0*xn+p1*xnm1-q1*ynm1-q2*ynm2;
      x_y[i]=yn;
    }
  }  //_1dof_Oscillator_Response


void main()  {
  // open input and output file
  f_in =fopen("tfft.dat","r");
  f_out=fopen("tshocksp.lst","w");
  // read number of input signal points in input file
  fscanf(f_in,"%d",&ndata);
  // read ndata couples T(i), Y(i) in input data file
  for (i=1; i<ndata+1; i++) {
    fscanf(f_in,"%lf %lf",&temp,&y);
	signal[i]=y;
    if (i==1)  tbegin=temp;
    if (i==ndata)  tend=temp;
  }
  fclose(f_in);
  // calculate sampling increment dt of input signal
  dt = (tend-tbegin)/(ndata-1);
  // input begin, end frequencies and frequency step df
  fbegin=10.0; fend=800.0; df=10.0;
  // calculate number of frequencies nfreq
  nfreq=(int) (1+floor((fend-fbegin)/df));
  // input damping factor dzeta
  dzeta=0.05;
  f=fbegin;
  // frequency loop
  for (i=1; i<nfreq+1; i++) {
    // call elementary oscillator response procedure for frequency f
    _1dof_Oscillator_Response(f,signal,response,dt,dzeta,ndata);
    // seek minimum and maximum values of response
    mini=response[1]; maxi=mini;
    for (j=2; j<ndata+1; j++) {
      if (response[j]<mini)  mini=response[j];
      if (response[j]>maxi)  maxi=response[j];
    }
    // store mini in SP1[i] and maxi in SP2[i] 
    SP1[i]=mini; SP2[i]=maxi;
    f+=df;
  }

  // print shock spectrum (negative and positive) to output file
  f=fbegin;
  fprintf(f_out,"\n");
  fprintf(f_out,"      Frequency     Shock Spectrum      Shock Spectrum   \n");
  fprintf(f_out,"        (Hz)           negative            positive      \n");
  fprintf(f_out,"  -------------------------------------------------------\n");
  for (i=1; i<nfreq+1; i++) {
    fprintf(f_out,"      %8.2f     %12.6f       %12.6f\n",f,SP1[i],SP2[i]);
    f=f+df;
  }

  fclose(f_out);
  printf("\n Results in file tshocksp.lst.\n");

}

// End of file tshocksp.cpp