'***************************************************************
'*   Program to demonstrate the Acceleration Shock Spectrum    *
'*                                                             *
'*                       Basic 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.                                               *
'***************************************************************
defint i-n
defdbl a-h,o-z

'i,j,ndata,nfreq : INTEGER
'df,dt,dzeta,f,fbegin,fend,temp,tbegin,tend : DOUBLE
'ymini,ymaxi : DOUBLE

dim  signal(1024),response(1024),SP1(100),SP2(100)

  cls
  'open input and output file
  open "tfft.dat" for input as #1
  open "tshocksp.lst" for output as #2
  'read number of input signal points in input file
  input #1, ndata
  'read ndata couples T(i), Y(i) in input data file
  for i=1 to ndata
    input #1, temp, signal(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 and frequency step df
  fbegin=10# : fend=800# : df=10#
  'calculate number of frequencies nfreq
  nfreq=1+INT((fend-fbegin)/df)
  'input damping factor dzeta
  dzeta=0.05#
  f=fbegin
  'frequency loop
  for i=1 to nfreq
    'call elementary oscillator response procedure for frequency f
    gosub 1000
    '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,using "      #####.##     #####.######       #####.######"; f; SP1(i); SP2(i)
    f=f+df
  next i

  close #2
  print
  print " Results in file tshocksp.lst."

END 'of main program


  '************************************************************
  '* 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  *
  '************************************************************
1000  'Subroutine _1dof_Oscillator_Response
    'ii:integer
    Pi=3.1415926535
    Omega=2#*Pi*f
    IF dzeta<1E-6 THEN dzeta=1E-6
    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=signal(1)
    yn=0#
    response(1)=yn
    for ii=2 to ndata
      ynm2=ynm1 : ynm1=yn : xnm1=xn
      xn=signal(ii)
      yn=p0*xn+p1*xnm1-q1*ynm1-q2*ynm2
      response(ii)=yn
    next ii
return

'End of file tshocksp.bas