Attribute VB_Name = "Module2"
Dim X(100) As Double, Y(100) As Double
Dim a As Double, b As Double, c As Double, d As Double
Dim dx As Double, xmn As Double, xmx As Double
Dim xx As Double, ymn As Double, ymx As Double
Dim iorder As Integer, iz As Integer
Dim m As Integer, n As Integer
Dim s1 As String * 12, s2 As String * 12
Dim buf As String, filename As String

'draw a plus sign at physical location a,b
Sub CrossXY(ByVal a#, ByVal b#)
  dx = (xmx - xmn) / 50
  dy = (ymx - ymn) / 50
  MoveXY a - dx, b
  LineXY 0, a + dx, b
  MoveXY a, b - dy
  LineXY 0, a, b + dy
End Sub

'*******************************************************
'* Least Squares of order 1 or 2 Demonstration Program *
'* --------------------------------------------------- *
'*   Reference: BASIC Scientific Subroutines, Vol. II  *
'*   By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
'*                                                     *
'*      Visual Basic 4,0 release by J-P Moreau, Paris  *
'*      (with graphic representation of results).      *
'*                  (www.jpmoreau.fr)                  *
'* --------------------------------------------------- *
'* FIRST SAMPLE RUN:                                   *
'*                                                     *
'* This program calculates a linear or parabolic least *
'* squares fit to a given data set.                    *
'*                                                     *
'* INSTRUCTIONS                                        *
'* ------------                                        *
'*                                                     *
'* The number of data coordinates provided must be     *
'* greater than three. Otherwise, a divide by zero     *
'* error may result.                                   *
'*                                                     *
'* Order of fit (1 or 2): 1                            *
'* Number of data points: 4                            *
'*                                                     *
'* There are two input options:                        *
'* 1. input coordinate pairs                           *
'* 2. first input the independant variable values,     *
'*    then input dependant ones.                       *
'*    Your choice (1 or 2): 1                          *
'*                                                     *
'* ? 1 1                                               *
'* ? 2 2                                               *
'* ? 3 3                                               *
'* ? 5 5.01                                            *
'*                                                     *
'* Fitted equation is:                                 *
'*   Y = -0.004571 + 1.002571 X                        *
'*                                                     *
'* Standard deviation of fit: 0.002928                 *
'*                                                     *
'* SECOND SAMPLE RUN:                                  *
'*                                                     *
'* Order of fit (1 or 2): 2                            *
'* Number of data points: 4                            *
'*                                                     *
'* There are two input options:                        *
'* 1. input coordinate pairs                           *
'* 2. first input the independant variable values,     *
'*    then input dependant ones.                       *
'*    Your choice (1 or 2): 1                          *
'*                                                     *
'* ? 1 1                                               *
'* ? 2 4                                               *
'* ? 3 9                                               *
'* ? 5 24.95                                           *
'*                                                     *
'* Fitted equation is:                                 *
'*   Y = -0.017727 + 0.022045 X + 0.994318 X^2         *
'*                                                     *
'* Standard deviation of fit: 0.004767                 *
'*                                                     *
'******************************************************}
'Note: use with graph2d.bas and lstsqr.frm


'*********************************************************
'*          Linear least squares subroutine              *
'* ----------------------------------------------------- *
'* In:  integer n := number of points                    *
'*      n values x[i), y(i) shared with main             *
'* Out: coefficients a,b of fit (a+b*x) shared with main *
'*      standard deviation d shared with main            *
'*********************************************************
Sub lstsqr1(ByVal n%)
Dim m As Integer
  a1 = 0: a2 = 0: b0 = 0: b1 = 0
  For m = 1 To n
    a1 = a1 + X(m)
    a2 = a2 + X(m) * X(m)
    b0 = b0 + Y(m)
    b1 = b1 + Y(m) * X(m)
  Next m
  a1 = a1 / n: a2 = a2 / n: b0 = b0 / n: b1 = b1 / n
  d = a1 * a1 - a2
  a = a1 * b1 - a2 * b0: a = a / d
  b = a1 * b0 - b1: b = b / d
  'Evaluation of standard deviation d (unbiased estim%ate)
  d = 0
  For m = 1 To n
    d1 = Y(m) - a - b * X(m)
    d = d + d1 * d1
  Next m
  d = Sqr(d / (n - 2))
End Sub

'*****************************************************************
'*           Parabolic least squares subroutine                  *
'* ------------------------------------------------------------- *
'* In:  integer n := number of points                            *
'*      n values x[i), y(i) shared with main                     *
'* Out: coefficients a,b,c of fit (a+b*x+c*x^2) shared with main *
'*      standard deviation d shared with main                    *
'*****************************************************************
Sub lstsqr2(ByVal n%)
Dim m As Integer
  a0 = 1: a1 = 0: a2 = 0: a3 = 0: a4 = 0
  b0 = 0: b1 = 0: b2 = 0
  For m = 1 To n
    a1 = a1 + X(m)
    a2 = a2 + X(m) * X(m)
    a3 = a3 + X(m) * X(m) * X(m)
    a4 = a4 + X(m) * X(m) * X(m) * X(m)
    b0 = b0 + Y(m)
    b1 = b1 + Y(m) * X(m)
    b2 = b2 + Y(m) * X(m) * X(m)
  Next m
  a1 = a1 / n: a2 = a2 / n: a3 = a3 / n: a4 = a4 / n
  b0 = b0 / n: b1 = b1 / n: b2 = b2 / n
  d = a0 * (a2 * a4 - a3 * a3) - a1 * (a1 * a4 - a2 * a3) + a2 * (a1 * a3 - a2 * a2)
  a = b0 * (a2 * a4 - a3 * a3) + b1 * (a2 * a3 - a1 * a4) + b2 * (a1 * a3 - a2 * a2)
  a = a / d
  b = b0 * (a2 * a3 - a1 * a4) + b1 * (a0 * a4 - a2 * a2) + b2 * (a1 * a2 - a0 * a3)
  b = b / d
  c = b0 * (a1 * a3 - a2 * a2) + b1 * (a2 * a1 - a0 * a3) + b2 * (a0 * a2 - a1 * a1)
  c = c / d
  'Evaluation of standard deviation d
  d = 0
  For m = 1 To n
    d1 = Y(m) - a - b * X(m) - c * X(m) * X(m)
    d = d + d1 * d1
  Next m
  d = Sqr(d / (n - 3))
End Sub

'***********************************************************
'* Driver subroutine for least square of order one or two, *
'* with data read from input text file.                    *
'***********************************************************
Sub Lstsqr()
Dim m As Integer, n As Integer, iorder As Integer, iz As Integer
  Form1.AutoRedraw = False
  Form1.Cls
  'Examples of valid input files:
  'lstsqr.dat           lstsqr1.dat
  '----------           -----------
  '1                    2
  '4                    4
  '1                    1
  '1, 1                 1, 1
  '2, 2                 2, 4
  '3, 3                 3, 9
  '5, 5.01              5, 24.95
  '-----------------------------------
  'select input file name (*.dat) from common
  'dialog box and open input file.
  Form1.OpenDialog.CancelError = False
  Form1.OpenDialog.Flags = cdlOFNCreatePrompr Or cdlOFNShowHelp
  Form1.OpenDialog.Filter = "Data files ( *.DAT ) | *.dat"
  Form1.OpenDialog.ShowOpen
  filename = Form1.OpenDialog.filename
  
  If Len(filename) = 0 Then
    MsgBox "No data file is open!", 48, "ERROR"
    Exit Sub
  Else
    Open filename For Input As #1
  End If

  'read from input file iorder (1 or 2)
  Input #1, iorder
  'read from input file number of points
  Input #1, n
  ReDim V(n)        'vector V is declared in graph2d.bas
  'read from input file input option (1 or 2)
  ' There are two input options:
  ' iz=1. input coordinate pairs (Example;1,2.5)
  ' iz=2. first input the independant variable values, then input dependant ones.
  Input #1, iz

  If iz = 2 Then
    'read data from input file (option 2)
    For m = 1 To n
      Input #1, X(m)
    Next m
    For m = 1 To n
      Input #1, Y(m)
    Next m
  Else
    'read data from input file (option 1)
    For m = 1 To n
      Input #1, X(m), Y(m)
    Next m
  End If
  
  Close #1   'close input file

  'Call linear or parabolic least squares subroutine
  If iorder = 1 Then
    lstsqr1 n
  Else
    lstsqr2 n
  End If

  'draw graph of results
  
  'find range in Ox and Oy
  For m = 0 To n - 1
    V(m) = X(m + 1)
  Next
  Minmax n
  xmn = vmini: xmx = vmaxi
  For m = 0 To n - 1
    V(m) = Y(m + 1)
  Next
  Minmax n
  ymn = vmini: ymx = vmaxi
    
  'init graphic window in physical coordinates
  '0=screen  10=predefined window #10 in pixels (see graph2d.bas)
  InitWindow 0, 10, xmn, xmx, ymn, ymx
  
  'Represent data points by plus signs
  For m = 1 To n
    CrossXY X(m), Y(m)
  Next m
  
  If iorder = 1 Then
    'Draw least squares line
    MoveXY xmn, b * xmn + a
    LineXY 0, xmx, b * xmx + a
  Else  'iorder=2
    'Draw least squares parabola a+bx+cx² with 50 points
    MoveXY xmn, c * xmn * xmn + b * xmn + a
    dx = (xmx - xmn) / 49: xx = xmn
    For m = 1 To 50
      xx = xx + dx
      LineXY 0, xx, c * xx * xx + b * xx + a
    Next m
  End If
  
  'Prepare graph caption
  If Abs(a) < 1 Then
    a = Int(a * 100000) / 100000
  End If
  If Abs(b) < 1 Then
    b = Int(b * 100000) / 100000
  End If
  s1 = Str$(a): s2 = Str$(b)
  If b > 0 Then
    buf = " Y = " & s1 & " +" & s2 & " X"
  Else
    buf = " Y = " & s1 & s2 & " X"
  End If
  
  If iorder = 2 Then
    If Abs(c) < 1 Then
      c = Int(c * 100000) / 100000
    End If
    s2 = Str$(c)
    If c > 0 Then
      buf = buf & " +" & s2 & " X^2"
    Else
      buf = buf & s2 & " X^2"
    End If
  End If
  
  'print graph caption and names of axes
  Legends 0, 10, buf, "X", "Y"
  
  'prepare standard deviation display
  If Abs(d) < 1 Then
    d = Int(d * 1000000) / 1000000
  End If
  buf = "Standard deviation =" & Str$(d)
  Form1.Font.Name = "Arial"
  Form1.Font.Size = 10
  Form1.Font.Bold = True
  Display 0, 3250, 550, buf
  
  'change label of button "Go"
  Form1.Command2.Caption = "Continue"

End Sub