/**************************************************************
*   Program to demonstrate the Acceleration Shock Spectrum    *
*                                                             *
*                         C++ version by J-P Moreau, Paris    *
*                                (www.jpmoreau.fr)            *
*                                                             *
*           (use with tgfft.mak, graph2d.cpp, basis_r.cpp(*)  *
*               and vmblock.cpp(*)).                          *
*                                                             *
* (*) to be found in MATRICES directory.                      *
* ----------------------------------------------------------- *
* 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>
#include <basis.h>     // ABS, MACH_EPS, REAL, WriteEnd(), WriteHead()   
#include <vmblock.h>   // dynamic allocations
#include <windows.h>

HDC  hdc;
RECT rect;
HPEN hpen, hpenOld;

//Functions used here of Graph2D.cpp
void Initwindow(int,double,double,double,double);
void MinMax(int,float *,double *,double *);
void TracerXY(int,float *,double,double);
void Legendes (int,char *,char *,char *,int,int,char *,int,int,char *);

//"home made" graphic commands for hdc environment used by above functions
void DrawPixel(int ix,int iy) {
  //sorry, no other available command found
  Rectangle(hdc,rect.left+ix,rect.top+iy,
	            rect.left+ix+2,rect.top+iy+1);
}

void Swap(int *i1,int *i2) {
  int it;
  it=*i1; *i1=*i2; *i2=it;
}

void DrawLine(int ix1,int iy1,int ix2,int iy2) {
  int i,il,ix,iy;
  double dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(double) (ix2-ix1)/il; 
  dy=(double) (iy2-iy1)/il;
  ix=ix1; iy=iy1;
  for (i=1;i<il+1; i++) {
	DrawPixel(ix,iy);
	ix=(int) (ix1+i*dx); 
	iy=(int) (iy1+i*dy);
  }
}

void OutText(int ix,int iy,char *s) {
  RECT rct;
  rct.left=ix; rct.right=ix+10*strlen(s);
  rct.top=iy-5; rct.bottom=iy+15;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}

#define  SIZE  1025

typedef double Vect[SIZE];
typedef float  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 Shock_spectrum()  {
  extern MaxX, MaxY;
  double spmin,spmax,temp;
  //client drawing zone in pixels
  MaxX=(rect.right-rect.left);
  MaxY=(rect.bottom-rect.top);
  // 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]=(float) mini; SP2[i]=(float) 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;
  }
  //restore last frequency
  f-=df;
  //shift index
  for (i=0; i<nfreq; i++) {
    SP1[i]=SP1[i+1];
    SP2[i]=SP2[i+1];
  }
  //draw positive and negative shock spectrum
  MinMax(nfreq,SP1,&spmin,&temp);
  MinMax(nfreq,SP2,&temp,&spmax);
  Initwindow(10,0,f,spmin,spmax);
  TracerXY(nfreq,SP1,0,f);
  TracerXY(nfreq,SP2,0,f);
  fclose(f_out);
  Legendes (10,"   SHOCK SPECTRUM   ","freq. hz","Acc. m/s2",
	        MaxX-360,50,"Numerical results in file gshocksp.lst",
			MaxX-250,70,"Input file: tfft.dat");
}

//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
   
    PAINTSTRUCT ps;

    switch (message) {
        case WM_PAINT:
	    hdc = BeginPaint(hwnd, &ps);
	    GetClientRect(hwnd, &rect);
            hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255));
            hpenOld = SelectObject(hdc,hpen);

            Shock_spectrum();
               
            SelectObject(hdc,hpenOld);
            DeleteObject(hpen);
            EndPaint(hwnd, &ps);
            return 0;
        case WM_DESTROY:
            PostQuitMessage(0);
            return 0;
    }
    return DefWindowProc(hwnd, message, wParam, lParam);
}

// main program
int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance,
    LPSTR lpszCmdParam, int nCmdShow)
{
    static char szAppName[] = "Gshocksp";
    HWND        hwnd;
    MSG         msg;
    WNDCLASS    wndclass;

    if (!hPrevInstance) {
        wndclass.style         = CS_HREDRAW | CS_VREDRAW;
        wndclass.lpfnWndProc   = WndProc;
        wndclass.cbClsExtra    = 0;
        wndclass.cbWndExtra    = 0;
        wndclass.hInstance     = hInstance;
        wndclass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
        wndclass.hCursor       = LoadCursor(NULL, IDC_ARROW);
        wndclass.hbrBackground = GetStockObject(WHITE_BRUSH);
        wndclass.lpszMenuName  = NULL;
        wndclass.lpszClassName = szAppName;

        RegisterClass(&wndclass);
    }

    hwnd = CreateWindow(szAppName,  // window class name
        " ACC. SHOCK SPECTRUM ",    // window caption
        WS_OVERLAPPEDWINDOW,        // window style                       
        CW_USEDEFAULT,              // initial x position
        CW_USEDEFAULT,              // initial y position
        CW_USEDEFAULT,              // initial x size
        CW_USEDEFAULT,              // initial y size
        NULL,                       // parent window handle
        NULL,                       // window menu handle
        hInstance,                  // program instance handle
        NULL);                      // creation parameters

   ShowWindow(hwnd, nCmdShow);
   UpdateWindow(hwnd);

   while (GetMessage(&msg, NULL, 0, 0)) {
       TranslateMessage(&msg);
       DispatchMessage(&msg);
   }
   return msg.wParam;
}

//end of file gshocksp.cpp