DECLARE SUB READATA (VIT() AS SINGLE, TS AS SINGLE, NDATA AS INTEGER)
DECLARE SUB OSCIL (FREQU AS SINGLE, VIT() AS SINGLE, ACC() AS SINGLE, TS AS SINGLE, DZETA AS SINGLE, NDATA AS INTEGER)
'******************************************************************
'* Calculate the acceleration of the sismic mass of an elementary *
'* oscillator (1 degree of freedom), the basis of which is submit-*
'* ted to a given speed, VIT(t). The input signal VIT is digitali-*
'* zed with a constant time step, TS. The table VIT(I) contains   *
'* NDATA speed values.                                            *
'*                                                                *
'*                      Quick Basic Version By J-P Moreau, Paris  *
'*                              (in simple precision).            *
'* -------------------------------------------------------------- *
'* SAMPLE RUN:                                                    *
'*                                                                *
'* Read data from file speed.dat                                  *
'* Write results to file response.txt.                            *
'******************************************************************
DEFINT I-N
NMAX = 2048  'Maximum number of points of input signal
DIM ACC(NMAX), VIT(NMAX) AS SINGLE
DIM NDATA AS INTEGER
DIM DZETA, F0, FMAX, T, TBEGIN, TEND, TS AS SINGLE

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

        PRINT #2,
        PRINT #2, " Mass Response to speed(T) at Basis"
        PRINT #2,

        F0 = 250!: DZETA = .05

' Frequency of oscillator
        PRINT #2, " FREQUENCY OF oscillator: " + STR$(F0) + " HERTZ"
        PRINT #2,

' rod damping
        PRINT #2, " DAMPING OF oscillator: " + STR$(DZETA)
        PRINT #2,

' READ MEASUREMENTS
        READATA VIT(), TS, NDATA

        CLOSE #1

        TBEGIN = 0: TEND = TBEGIN + (TS * (NDATA - 1))

        FMAX = .5 / TS
        PRINT #2, " NYQUIST FREQUENCY:" + STR$(FMAX)

' Compute mass response
        OSCIL F0, VIT(), ACC(), TS, DZETA, NDATA

        PRINT #2,
        PRINT #2, "    TIME (S)    MASS ACC. (M/S2) "
        PRINT #2, "  ------------  -----------------"

        T = TBEGIN

        FOR I = 1 TO NDATA
            PRINT #2, " " + STR$(T);
            PRINT #2, TAB(18); ACC(I)
            T = T + TS
        NEXT I

        CLOSE #2

        CLS
        PRINT
        PRINT " Results in file reponse.txt..."
        PRINT

    END 'of main program




SUB OSCIL (FREQU AS SINGLE, VIT() AS SINGLE, ACC() AS SINGLE, TS AS SINGLE, DZETA AS SINGLE, NDATA AS INTEGER)
    '*************************************************************
    '* 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, DELTA, EXPA, GOMEGA, OMEGA, OMEGA2, PI, T AS SINGLE
    DIM COSARG, Q0, Q1, Q2, QSI, R0, R1, R2, R3, SINARG AS SINGLE
    DIM XPN, XPNP1, YPN, YPNP1, YPPN, YPPNP1 AS DOUBLE

    PI = 4! * ATN(1!)
    T = TS
    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

SUB READATA (VIT() AS SINGLE, TS AS SINGLE, NDATA AS INTEGER)
    '------------------------------------------------
    'Read ndata: number of points of input signal
    '        ts: sampling duration of signal in sec.
    '    ACC(i): Table with ndata acceleration values
    '------------------------------------------------
    DIM TEMP AS SINGLE
    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