Attribute VB_Name = "Module2"

Public freq As Double      'frequency of oscillator
Public dzeta As Double     'reduced damping of oscillator
Public dt As Double        'Increment axix Ox
Public ti As Double, tf As Double 'starting, ending time

Public A() As Double       'Complex vector
Public ndata As Integer    'number of points
Public numfich As Integer  'number of open file(s), 1 or 2.

Public filename1 As String * 80  'name of input file
Public filename2 As String * 80  'name of output file


'*******************************************************************
'* This subroutine calculates the sismic mass acceleration of a one*
'* degree of freedom oscillator the basis of which is submitted to *
'* an input acceleration x". The input signal x" is sampled with a *
'* fixed period dt, stored in vector V() composed of ndata accele- *
'* ration values.                                                  *
'* Inputs:                                                         *
'*       fr#  : eigenfrequency of oscillator in hertz,             *
'*       dt#  : sampling increment in time axis,                   *
'*       dzeta: reduced damping (usually 0.05),                    *
'*       n%   : number of points of signal V(),                    *
'*       V    : n% ordinates of x" (vector defined in graph2d.bas).*
'* Output:                                                         *
'*       A    : n% ordinates of mass response are put in complex   *
'*              vector A(i,0) defined in Declarations section.     *
'*              (here the imaginary part is not used).             *
'* --------------------------------------------------------------- *
'*  Visual Basic 4.0, 04/15/1996 after a pascal release for        *
'*  Windows 3.1, 02/04/1990.                      By J-P Moreau    *
'*                                              (www.jpmoreau.fr)  *
'*******************************************************************
Sub CalculateResponse(ByVal fr#, ByVal dt#, ByVal dzeta As Double, ByVal n%)
Const pi = 3.1415926535
  omega = 2 * pi * fr#
  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(0)
  yn = 0#
  A(0, 0) = yn
  For i% = 1 To n% - 1
    ynm2 = ynm1: ynm1 = yn: xnm1 = xn
    xn = V(i%)
    yn = p0 * xn + p1 * xnm1 - q1 * ynm1 - q2 * ynm2
    A(i%, 0) = yn
  Next
 
End Sub

'*****************************************************
'* This subroutine:                                  *
'* 1) reads a curve F(t) stored on disk (*.SGN),     *
'* 2) draws to screen F(t) in window #5,             *
'* 3) calculates mass response of oscillator with    *
'*    frequency, freq and reduced damping, dzeta,    *
'*    calling subroutine CalculateResponse,          *
'* 4) draws to screen response(t) in window #3,      *
'* 5) stores to disk mass response at format SGN.    *
'***************************************************** 
Sub Response1dof(flp As Integer)
Dim f As String
' flp=0: screen - flp<>0: printer
  ReadInput "*.SGN", False

' read frequency, dzeta from dialog box
  If Not flp Then
    Form3.Show 1                      ' 1 = modal dialog
    Form1.Cls
  End If
  
  ReDim A(ndata, 2)
    
  Form1.AutoRedraw = False
  'Form1.Cls
  
  CurveXY 5, ndata, ti, tf, flp, True
  Legends flp, 5, "Basis acceleration", nomx$, nomy$

  CalculateResponse freq, dt, dzeta, ndata

  For i% = 0 To ndata - 1
    V(i%) = A(i%, 0)
  Next
  CurveXY 3, ndata, ti, tf, flp, True
  Legends flp, 3, "Mass response", nomx$, nomy$
    
  If Not flp Then
    f = "1dof acc. Response ( a*.SGN ) | a*.sgn"
    Form1.OpenDialog.Filter = f
    Form1.OpenDialog.ShowSave
    nom1$ = Form1.OpenDialog.filename
    nom1$ = CFileTitle(nom1$)
    If Mid$(nom1$, Len(nom1$) - 3, 1) <> "." Then nom1$ = nom1$ + ".SGN"
    titre$ = "Response of " + titre$
    DiskWriteData nom1$, ndata, ti, tf, titre$, nomx$, nomy$
  End If
  
  Form1.Command1.Caption = "Continue"
  
End Sub

'*****************************************************
'* Read a curve F(t) stored on disk at format *.SGN. *
'* Nubers may not have a decimal comma, but a deci-  *
'* mal point, else erreor.                           *
'* ------------------------------------------------- *
'* 1st line: ndata, ti, tf                           *
'* (number of points, starting, final abscissas)     *
'* 2nd line: curve caption                           *
'* 3rd line: name axis Ox                            *
'* 4th line: name axis Oy                            *
'* ndata following lines: ordinate values (double)   *
'* (exponant notation accepted).                     *
'*****************************************************
Sub DiskReadData(nom$)
Dim nfich

  nfich = FreeFile
  Open nom$ For Input As nfich
  Input #nfich, ndata, ti, tf
  Input #nfich, titre
  Input #nfich, nomx
  Input #nfich, nomy
  ReDim V(ndata)
  For i% = 0 To ndata - 1
    Input #nfich, V(i%)
  Next
  Close nfich

End Sub

'***************************************************
'* Store a curve F(t) to disk at format *.SGN.     *
'* Numbers with a decimal comma are converted to a *
'* decimal point by function CNumber$.             *
'***************************************************
Sub DiskWriteData(nom$, ByVal n%, ByVal xi#, ByVal xf#, titre$, nomx$, nomy$)
Dim nfich

  nfich = FreeFile
  Open nom$ For Output As nfich
  Print #nfich, n%, CNumber$(xi#), CNumber$(xf#)
  Print #nfich, titre$
  Print #nfich, nomx$
  Print #nfich, nomy$
  For i% = 0 To n% - 1
    Print #nfich, CNumber$(V(i%))
  Next
  Close nfich

End Sub

'******************************************************
'* convert a number with decimal comma into a number  *
'* with a decimal dot before writing to a sequential  *
'* file.                                              *
'******************************************************
Function CNumber$(X As Double)
  xx$ = Str$(X)
  For i% = 1 To Len(xx$)
    If Mid$(xx$, i%, 1) = "," Then Mid$(xx$, i%, 1) = "."
  Next
  CNumber$ = xx$
End Function


' read a curve  F(t) from disk
Sub ReadInput(s As String, flag As Integer)
On Error GoTo cancel
  If numfich = 0 Then OpenFile s
  DiskReadData (filename1)
  dt = (tf - ti) / (ndata - 1)
  Form1.CurrentX = 1: Form1.CurrentY = 1
  If flag <> 0 Then
    Form1.Print
    Form1.Print
    Form1.Print "   Number of points read: "; ndata
    Form1.Print "   Curve caption        : "; titre
    Form1.Print "   Name axis Ox         : "; nomx
    Form1.Print "   Name axis Oy         : "; nomy
    Form1.Print "   Increment axis Ox    : "; Str$(dt)
    Pause
  End If
cancel:  'no action
End Sub

Sub OpenFile(s As String)
Dim f As String
  If s = "*.SGN" Then
    f = "Curves F(t) ( *.SGN ) | *.sgn"
  ElseIf s = "sp*.SGN" Then
    f = "Shock Spectrum ( sp*.SGN ) | sp*.sgn"
  ElseIf s = "a*.SGN" Then
    f = "Acc. Response 1dof ( a*.SGN ) | a*.sgn"
  ElseIf s = "v*.SGN" Then
    f = "Basis Speed ( v*.SGN ) | v*.sgn"
  Else
    f = "Curves F(t) ( *.TXT ) | *.TXT"
  End If
  Form1.OpenDialog.Filter = f
  If numfich = 0 Then   ' aucun fichier ouvert
    Form1.OpenDialog.CancelError = False
    Form1.OpenDialog.Flags = cdlOFNCreatePrompt Or cdlOFNShowHelp
    Form1.OpenDialog.ShowOpen
    filename1 = Form1.OpenDialog.filename
    numfich = numfich + 1
  ElseIf numfich = 1 Then  ' fichier de sortie
    Form1.OpenDialog.CancelError = True
    Form1.OpenDialog.Flags = cdlOFNCreatePrompt Or cdlOFNShowHelp
    Form1.OpenDialog.ShowOpen
    filename2 = Form1.OpenDialog.filename
    numfich = numfich + 1
  Else
    MsgBox "Both files are open!", 48, "WARNING"
    Beep
  End If
  'print name(s) of open file(s) downwards
  Form1.CurrentX = 50: Form1.CurrentY = MaxY + 300
  Form1.Print CFileTitle$(filename1)
  Form1.CurrentX = MaxX / 2 + 50: Form1.CurrentY = MaxY + 500
  Form1.Print CFileTitle$(filename2)
End Sub

'eliminate path from file name
Function CFileTitle$(filename As String)
  For i = Len(filename) To 1 Step -1
    If Mid$(filename, i, 1) = "\" Then
      n = i
      CFileTitle$ = Right$(filename, Len(filename) - n)
      Exit Function
    End If
  Next
End Function

'end of file response.bas