SUB mathfitcos (n%, y!(), a!, b!, c!, phi!, s!)

'****************************************************************************
'This subroutine fits n data points to a sine wave having 2 cycles in the
'n points.
'The fit is to y = a!*sin(2*2pi*i/n) + b!*cos(2*2pi*i/n)
'                = c!*cos(2*2pi*i/n + phi!)
'The rms of the residuals is also returned.
'y = c!*cos(2*2pi*i/n - nphi!)
'nphi! = -phi!
'y = c!*cos(2*2pi*i/n)*cos(nphi!) + c!*sin(2*2pi*i/n)*sin(nphi!)
'a! = c!*sin(nphi!), b! = c!*cos(nphi!)
'c!^2 = a!^2 + b!^2, c! = sqr(a!^2 + b!^2)
'a!/b! = tan(nphi!), nphi! = atan(a!/b!), consider different quadrants
'
'Input:
' n%, the number of data points
' y!(0 to n%-1), the data points
'
'Output:
' a!, the coeffecient of the sin term from the fit
' b!,            "           cos        "
' c!, the fitted amplitude
' phi!, the fitted phase in degrees
' s!, the rms of the residuals
'
'Zachary Wolf
'8/19/94
'****************************************************************************

'Parameters
pi! = 3.1415926#

'Compute various sums
Sss! = 0!
Scc! = 0!
Ssc! = 0!
Sys! = 0!
Syc! = 0!
FOR i% = 0 TO n% - 1
 Sss! = Sss! + SIN(2 * 2 * pi! * i% / n%) * SIN(2 * 2 * pi! * i% / n%)
 Scc! = Scc! + COS(2 * 2 * pi! * i% / n%) * COS(2 * 2 * pi! * i% / n%)
 Ssc! = Ssc! + SIN(2 * 2 * pi! * i% / n%) * COS(2 * 2 * pi! * i% / n%)
 Sys! = Sys! + y!(i%) * SIN(2 * 2 * pi! * i% / n%)
 Syc! = Syc! + y!(i%) * COS(2 * 2 * pi! * i% / n%)
NEXT i%

'Compute the determinant
d! = Sss! * Scc! - Ssc! * Ssc!

'Compute the sin and cos coefficients
a! = (Sys! * Scc! - Ssc! * Syc!) / d!
b! = (Sss! * Syc! - Sys! * Ssc!) / d!

'Compute the amplitude coefficient
c! = SQR(a! ^ 2 + b! ^ 2)

'Compute the phase coefficient, nphi! = -phi!
IF a! > 0 AND b! > 0 THEN nphi! = ATN(a! / b!)
IF a! > 0 AND b! < 0 THEN nphi! = pi! - ATN(-a! / b!)
IF a! < 0 AND b! > 0 THEN nphi! = -ATN(-a! / b!)
IF a! < 0 AND b! < 0 THEN nphi! = -pi! + ATN(a! / b!)
IF a! = 0 AND b! >= 0 THEN nphi! = 0!
IF a! = 0 AND b! < 0 THEN nphi! = pi!
IF b! = 0 AND a! > 0 THEN nphi! = pi! / 2
IF b! = 0 AND a! < 0 THEN nphi! = -pi! / 2
phi! = -nphi!
phi! = phi! * 180! / pi!        'convert to degrees

'Compute the rms of the residuals
s! = 0!
FOR i% = 0 TO n% - 1
 s! = s! + (y!(i%) - a! * SIN(2 * 2 * pi! * i% / n%) - b! * COS(2 * 2 * pi! * i% / n%)) ^ 2
NEXT i%
s! = SQR(s! / n%)

END SUB

SUB mathlinfit (n%, x!(), y!(), m!, b!)

'****************************************************************************
'This subroutine does a least squares linear fit to the (x,y) data points.
'It returns the slope and intercept of the line.
'
'Input:
' n%, the number of (x,y) points
' x!(1 to n%), the x values
' y!(1 to n%), the y values
'
'Output:
' m!, the fitted line's slope
' b!, the fitted line's intercept
'
'Zachary Wolf
'5/12/94
'****************************************************************************

'Compute various sums
Sxy! = 0!
Sxx! = 0!
Sx! = 0!
Sy! = 0!
FOR i% = 1 TO n%
 Sxy! = Sxy! + x!(i%) * y!(i%)
 Sxx! = Sxx! + x!(i%) * x!(i%)
 Sx! = Sx! + x!(i%)
 Sy! = Sy! + y!(i%)
NEXT i%

'Compute the determinant
d! = n% * Sxx! - Sx! * Sx!

'Compute the slope and intercept
m! = (Sxy! * n% - Sx! * Sy!) / d!
b! = (Sxx! * Sy! - Sxy! * Sx!) / d!

END SUB