Attribute VB_Name = "Module2"
'***************************************************************
'*   Program to demonstrate the Acceleration Shock Spectrum    *
'*                                                             *
'*            Visual Basic 4,0 release by J-P Moreau, Paris    *
'*                     (with graphic display).                 *
'*                        (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.                                               *
'***************************************************************

Public df As Double, dt As Double, fbegin As Double, fend As Double
Public dzeta As Double

Dim nfreq As Integer
Dim temp As Double, tbegin As Double, tend As Double

Dim Response(), SP1(), SP2()
Dim ndata As Integer
Dim f As Double, Q As Double, ymini As Double, ymaxi As Double

Sub ShockSp()
  Dim i As Integer
  Dim spmin As Double, spmax As Double
  Form1.AutoRedraw = False
  Form1.Cls
  'open input and output file
  Open "c:\temp\tgfft.dat" For Input As #1
  Open "c:\temp\tgshocksp.lst" For Output As #2
  'read number of input signal points in input file
  Input #1, ndata
  ReDim V(ndata), Response(ndata)
  'read ndata couples T(i), Y(i) in input data file
  For i = 1 To ndata
    Input #1, temp, V(i)
    If i = 1 Then tbegin = temp
    If i = ndata Then tend = temp
  Next i
  Close #1
    
  'calculate sampling increment dt of input signal
  dt = (tend - tbegin) / (ndata - 1)
  'input begin, end frequencies, frequency step df, damping dzeta
  Form2.Show 1      '1=modal dialog
  
  Display 0, 200, MaxY + 1000, "Computing shock spectrum..."
  'calculate number of frequencies nfreq
  nfreq = 1 + Int((fend - fbegin) / df)
  ReDim SP1(nfreq), SP2(nfreq)
  'input damping factor dzeta
  f = fbegin
  'frequency loop
  For i = 1 To nfreq
    'call elementary oscillator response procedure for frequency f
    MassResponse
    'seek minimum and maximum values of response
    ymini = Response(1): ymaxi = ymini
    For j = 2 To ndata
      If Response(j) < ymini Then ymini = Response(j)
      If Response(j) > ymaxi Then ymaxi = Response(j)
    Next j
    'store mini in SP1(i) and maxi in SP2(i)
    SP1(i) = ymini: SP2(i) = ymaxi
    f = f + df
  Next i

  'print shock spectrum (negative and positive) to output file
  f = fbegin
  Print #2, "      Frequency     Shock Spectrum      Shock Spectrum   "
  Print #2, "        (Hz)           negative            positive      "
  Print #2, "  -------------------------------------------------------"
  For i = 1 To nfreq
    Print #2, "        "; Int(f * 100) / 100;
    Print #2, "        "; Int(SP1(i) * 1000000) / 1000000;
    Print #2, "        "; Int(SP2(i) * 1000000) / 1000000
    f = f + df
  Next i

  Close #2
  
  'seek minimum and maximum of SP1, SP2
  For i = 1 To nfreq
    If SP1(i) < spmin Then spmin = SP1(i)
    If SP2(i) > spmax Then spmax = SP2(i)
  Next i
  'draw shock spectrum
  Form1.Cls
  InitWindow 0, 10, fbegin, fend, spmin, spmax
  'positive values
  For i = 0 To nfreq - 1
    V(i) = SP2(i + 1)
  Next i
  CurveXY 10, nfreq, fbegin, fend, 0, 0
  'negative values
  For i = 0 To nfreq - 1
    V(i) = SP1(i + 1)
  Next i
  CurveXY 10, nfreq, fbegin, fend, 0, 0
  'print caption of graph and names of axes
  Legends 0, 10, "Acc. Shock spectrum", "Freq. hz", "Acc. m/s2"

  'print value of FNyquist and df
  Fnyq = 1# / (2# * dt): Fnyq = Int(Fnyq * 10) / 10
  s$ = " FNyquist =" & Str$(Fnyq)
  Display 0, MaxX - 900, 550, s$
  s$ = " DF= " & Str$(df)
  Display 0, MaxX - 900, 800, s$
  Q = 1# / (2# * dzeta)
  s$ = " Q= " & Str$(Int(Q * 100) / 100)
  Display 0, MaxX - 900, 1050, s$
  
  Display 0, 200, MaxY + 1000, "Numerical results in file c:\temp\tgshocksp.lst."

  Form1.Command1.Caption = "Continue"

End Sub  'ShockSp


  '************************************************************
  '* 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:                                                  *
  '*         f............: eigen frequency of oscillator     *
  '*         signal.......: input acceleration signal x"(t)   *
  '*         dt...........: time step of signal (constant)    *
  '*         dzeta........: reduced damping factor            *
  '*         ndata........: number of points of signal        *
  '* OUTPUT:                                                  *
  '*         response.....: table containing acceleration     *
  '*                        response of oscillator (N points) *
  '*                                                          *
  '*                      Basic version by J-P Moreau, Paris  *
  '************************************************************
Sub MassResponse()
    Dim ii As Integer
    Dim omega As Double, Pi As Double
    Dim dQ2 As Double, Q As Double
    Pi = 3.1415926535
    omega = 2# * Pi * f
    If dzeta < 0.000001 Then dzeta = 0.000001
    Q = 1# / (2# * dzeta)
    dQ2 = 2# * Q * Q
    delta = Sqr(2# * dQ2 - 1#)
    p0 = omega * dt / Q
    q2 = Exp(-p0)
    sq2 = Sqr(q2)
    arg = 0.5 * p0 * delta
    cosA = Cos(arg)
    q1 = -2# * sq2 * cosA
    p1 = p0 * sq2 * ((dQ2 - 1#) * Sin(arg) / delta - cosA)
    ynm1 = 0#
    xn = V(1)
    yn = 0#
    Response(1) = yn
    For ii = 2 To ndata
      ynm2 = ynm1: ynm1 = yn: xnm1 = xn
      xn = V(ii)
      yn = p0 * xn + p1 * xnm1 - q1 * ynm1 - q2 * ynm2
      Response(ii) = yn
    Next ii
End Sub

'End of file tgshocksp.bas