DECLARE SUB OSCIL (FREQU!, VIT!(), ACC!(), T!, DZETA!, NDATA%)
DECLARE SUB EXTREM (ACC!(), YMIN!, YMAX!, NDATA%)
DECLARE SUB READATA (VIT!(), TS!, NDATA%)
DECLARE SUB SPEC (VIT!(), TS!, FMIN!, FMAX!, DZETA!, NDATA%, NFREQ%, S!())
'*************************************************************
'* Calculate the acceleration shock spectrum of a speed V(t) *
'* --------------------------------------------------------- *
'* SAMPLE RUN:                                               *
'* Input speed: in file speed.dat                            *
'* Output shock spectrum in file shocksp1.txt                *
'*                                                           *
'*                 Quick Basic Version By J-P Moreau, Paris. *
'*                             (www.jpmoreau.fr)             *
'*************************************************************
DEFINT I-N
NMAX = 2048
DIM VIT(NMAX), S(NMAX)
'Dim NDATA, NFREQ as integer
'Dim DF, dzeta, F, FMAX, FMIN, TS As Single

    OPEN "speed.dat" FOR INPUT AS #1
    OPEN "shocksp1.txt" FOR OUTPUT AS #2

    PRINT #2,
    PRINT #2, " Acceleration Shock Spectrum of a speed(T) at Basis"
    PRINT #2,

    DZETA = .05   'damping coefficient

' READ INPUT SPEED
    READATA VIT(), TS, NDATA

    CLOSE #1

    FMAX = .5 / TS
    PRINT #2, " NYQUIST FREQUENCY: "; FMAX

    FMIN = 10!: FMAX = 1000!: NFREQ = 991

'Calculate shock spectrum S
    SPEC VIT(), TS, FMIN, FMAX, DZETA, NDATA, NFREQ, S()

'Write spectrum to output file
    DF = (FMAX - FMIN) / (NFREQ - 1)
    PRINT #2,
    PRINT #2,
    PRINT #2, " THE ACCELERATION SHOCK SPECTRUM  IS COMPUTED FROM"
    PRINT #2, " THE FILTERED BASIS SPEED IF IT EXISTS, ELSE FROM"
    PRINT #2, " THE UNFILTERED BASIS SPEED."
    PRINT #2,
    PRINT #2, " FREQUENCY STEP OF SHOCK SPECTRUM: "; DF
    PRINT #2,
    PRINT #2,
    PRINT #2, "   FREQUENCY HZ  SPECTRUM +     SPECTRUM -  IN M/S2"
    PRINT #2, "   ------------  -----------    -------------------"

    FOR I = 1 TO NFREQ
	F = FMIN + (I - 1) * DF
	PRINT #2, TAB(4); F;
	PRINT #2, TAB(18); S(2 * I - 1);
	PRINT #2, TAB(36); S(2 * I)
    NEXT I
    PRINT #2,
    CLOSE #2

'Closing section
    CLS
    PRINT
    PRINT " Results in file shockspl.txt."
    PRINT

END 'of main program

SUB EXTREM (ACC(), YMIN, YMAX, NDATA)
    '-----------------------------------------------------------
    ' SEEK MINIMUM AND MAXIMUM OF A TABLE ACC(I)
    '
    ' Inputs:
    '            ACC(I)  Table with NDATA values of acceleration
    '            NDATA   Number of points of ACC(I)
    ' Outputs:
    '            YMIN    minimum value of table ACC
    '            YMAX    maximum value of table ACC
    '-----------------------------------------------------------
    YMIN = 0: YMAX = 0
    FOR I = 1 TO NDATA
	IF ACC(I) > YMAX THEN YMAX = ACC(I)
	IF ACC(I) < YMIN THEN YMIN = ACC(I)
    NEXT I

END SUB 'EXTREM

SUB OSCIL (FREQU, VIT(), ACC(), T, DZETA, NDATA)
    '*************************************************************
    '* This subroutine calculates the acceleration of the sismic *
    '* mass of an elementary oscillator (1 degree of freedom),   *
    '* the basis of which is submitted to a given speed, VIT(t). *
    '* The input signal VIT is digitalized with a constant time  *
    '* step, TS. The table VIT(I) contains NDATA speed values.   *
    '* --------------------------------------------------------- *
    '* INPUTS:                                                   *
    '*         FREQU........: eigen frequency of oscillator      *
    '*         VIT..........: input speed signal of basis VIT(t) *
    '*         TS...........: time step of signal (constant)     *
    '*         DZETA........: reduced damping factor             *
    '*         NDATA........: number of points of signal         *
    '* OUTPUT:                                                   *
    '*         ACC..........: table containing acceleration res- *
    '*                        ponse of oscillator (NDATA points) *
    '*                                                           *
    '*             VB 2008 Express Version By J-P Moreau, Paris. *
    '*************************************************************
    'Dim ARG, COSARG, DELTA, EXPA, GOMEGA, OMEGA, OMEGA2, SINARG As single
    'Dim Q0, Q1, Q2, QSI, R0, R1, R2, R3, XPNP1, YPNP1, YPPNP1 As single
    'Dim PI, XPN, YPN, YPPN As single
    PI = 4! * ATN(1!)
    OMEGA = 2 * PI * FREQU
    OMEGA2 = OMEGA * OMEGA
    DELTA = SQR(1 - DZETA * DZETA)
    EXPA = EXP(-OMEGA * DZETA * T)
    GOMEGA = OMEGA * DELTA
    ARG = GOMEGA * T
    SINARG = SIN(ARG)
    COSARG = COS(ARG)
    QSI = DZETA / DELTA
    Q0 = (COSARG - QSI * SINARG) * EXPA
    Q1 = OMEGA2 * SINARG / GOMEGA * EXPA
    Q2 = (1# + (QSI * SINARG - COSARG) * EXPA) / T
    R1 = SINARG / GOMEGA * EXPA
    R0 = COSARG * EXPA + DZETA * OMEGA * R1
    R2 = (1 - DZETA * OMEGA * T) * R1 / T - COSARG * EXPA
    R3 = 1 - R1 / T
    XPNP1 = VIT(1)
    YPPNP1 = 0
    YPNP1 = 0
    ACC(1) = Q2 * XPNP1
    FOR I = 2 TO NDATA
	YPN = YPNP1
	YPPN = YPPNP1
	XPN = XPNP1
	XPNP1 = VIT(I)
	YPPNP1 = Q0 * YPPN + Q1 * (XPN - YPN) + Q2 * (XPNP1 - XPN)
	ACC(I) = YPPNP1
	YPNP1 = R0 * YPN + R1 * YPPN + R2 * XPN + R3 * XPNP1
    NEXT I

END SUB 'OSCIL

DEFSNG I-N
SUB PRINTX

END SUB

DEFINT I-N
SUB READATA (VIT(), TS, NDATA)
    '------------------------------------------------
    'Read ndata: number of points of input signal
    '        ts: sampling duration of signal in sec.
    '    VIT(i): Table with ndata speed values
    '------------------------------------------------
    INPUT #1, NDATA
    PRINT #2, STR$(NDATA) + " POINTS READ."
    INPUT #1, TS
    PRINT #2, " SAMPLING DURATION: " + STR$(TS)
    FOR I = 1 TO NDATA
	INPUT #1, TEMP, VIT(I)
    NEXT I

END SUB

SUB SPEC (VIT(), TS, FMIN, FMAX, DZETA, NDATA, NFREQ, S())
    '-------------------------------------------------------------
    ' ACCELERATION SHOCK SPECTRUM OF A SPEED
    '
    ' Inputs:
    '            VIT(I)  speed of basis (from deconvolution or else)
    '            TS      sampling time of speed
    '            FMIN    begin frequency of spectrum (<>0)
    '            FMAX    end frequency of spectrum (<=Nyquist)
    '            DZETA   relative damping factor
    '            NDATA   number of points of VIT(I)
    '            NFREQ   number of frequencies of spectrum (<=1024)
    ' Output:
    '            S(I)    shock spectrum (negative & positive)
    '-------------------------------------------------------------
    NMAX = 2048
    DIM ACC(NMAX)
    DELTAF = (FMAX - FMIN) / (NFREQ - 1)
    ' The algorithm fails if FMIN=0
    IF FMIN < DELTAF THEN FMIN = DELTAF
    FOR II = 1 TO NFREQ
	FR = FMIN + (II - 1) * DELTAF
	OSCIL FR, VIT(), ACC(), TS, DZETA, NDATA
	EXTREM ACC(), YMIN, YMAX, NDATA
	S(2 * II - 1) = YMAX
	S(2 * II) = YMIN
    NEXT II

END SUB