DECLARE SUB hp3457extsig (sig%)
DECLARE SUB fileopenout (destdir$, basename$, ext$, filename$, ok$)
DECLARE SUB coilcalcvn (ns%, v!(), vm!(), vp!())
DECLARE SUB coilfitv2 (n%, y!(), a!, b!, c!, phi!, s!)
DECLARE SUB coilgetf (f!)
DECLARE SUB coilgetv (ns%, v!(), sv!())
DECLARE SUB coillogf (logfile$, freq!)
DECLARE SUB coillogfitv2 (logfile$, a!, b!, c!, phi!, s!)
DECLARE SUB coillogplotv (logfile$, ns%, v!(), sv!(), vm!(), vp!(), vres!())
DECLARE SUB coillogv (logfile$, rawptr%)
DECLARE SUB coillogvsamp (logfile$, rawptr%, ns%, v!())
DECLARE SUB hp3457cf (c%, freq!)
DECLARE SUB hp3457cnvtrig (c%, n%, v!())
DECLARE SUB FFTCalc (xreal!(), yimag!(), numdat%)

DECLARE FUNCTION FFTMagnitude! (xr(), yi(), n%, i%)
DECLARE FUNCTION FFTPhase! (xr(), yi(), n%, i%)

'****************************************************************************
'Module COIL.BAS
'
'Zachary Wolf
'1/4/95
'****************************************************************************

'Open the parameter file
REM $INCLUDE: 'param.inc'

SUB coilcalcvn (ns%, v!(), vm!(), vp!())

'****************************************************************************
'This subroutine calculates the magnitude and phase of each voltage harmonic.
'
'Input:
' ns%, the number of encoder pulses per revolution
' v!(0 to ns%-1), the voltage samples for one coil revolution
'
'Output:
' vm!(0 to ns%/2), the magnitude of the voltage at each harmonic
' vp!(0 to ns%/2), the phase of the voltage at each harmonic in degrees -180 to 180
'
'Zachary Wolf
'7/9/94
'****************************************************************************

'Initialize
DIM xr!(0 TO ns% - 1)
DIM yi!(0 TO ns% - 1)
FOR i% = 0 TO ns% - 1
 xr!(i%) = v!(i%)
 yi!(i%) = 0!
NEXT i%

'Get the voltage at each harmonic using an FFT
CALL FFTCalc(xr!(), yi!(), ns%)

'Calculate the magnitude and phase of each voltage harmonic
FOR i% = 0 TO ns% / 2
 vm!(i%) = FFTMagnitude(xr!(), yi!(), ns%, i%)
 vp!(i%) = FFTPhase(xr!(), yi!(), ns%, i%)
NEXT i%

'Convert from radians to degrees
'Convert from (0,360) to (-180,180)
FOR i% = 0 TO ns% / 2
 vp!(i%) = vp!(i%) * 180! / 3.141592654#
 IF vp!(i%) > 180! THEN vp!(i%) = -(360! - vp!(i%))
NEXT i%

END SUB

SUB coilfitv2 (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%)

'Log the results
CALL coillogfitv2(logfileP$, a!, b!, c!, phi!, s!)

END SUB

SUB coilgetf (f!)

'****************************************************************************
'This subroutine is used to get the rotation frequency of a measurement coil.
'
'Output:
' f!, the average of the frequency measurements
'
'Zachary Wolf
'10/6/94, 11/15/94
'****************************************************************************

'Simplify the notation
ns% = nsampleP%            '# encoder pulses per revolution
cf% = freqhpchanP%         'HP3457 channel for the coil frequency

'Perform the frequency measurement
CALL hp3457cf(cf%, f!)

'Take out ns% pulses per revolution
f! = f! / ns%

'Log the coil rotation frequency
CALL coillogf(logfileP$, f!)

END SUB

SUB coilgetv (ns%, v!(), sv!())

'****************************************************************************
'This subroutine gets the triggered voltages from the coil.
'
'Input:
' ns%, the number of voltage samples per revolution
'
'Output:
' v!(0 to ns%-1), the average voltage at each sample point
' sv!(0 to ns%-1), the standard deviation of v at each sample point
'
'Zachary Wolf
'8/14/94, 11/15/94
'****************************************************************************

'Simplify the notation
c% = coilhpchanP%
nr% = nrotaveP%

'Initialize data arrays
DIM vnr!(0 TO nr% * ns% - 1)        'voltage samples
DIM v2d!(0 TO ns% - 1, 1 TO nr%)    'put samples in 2d array

'Make the measurement over nr% rotations
CALL hp3457cnvtrig(c%, nr% * ns%, vnr!())

'Put the samples in a 2d array for ease of handling
FOR i% = 0 TO ns% - 1
FOR j% = 1 TO nr%
 v2d!(i%, j%) = vnr!(i% + (j% - 1) * ns%)
NEXT j%
NEXT i%

'Compute the average at each sample point
FOR i% = 0 TO ns% - 1
 v!(i%) = 0!
FOR j% = 1 TO nr%
 v!(i%) = v!(i%) + v2d!(i%, j%)
NEXT j%
 v!(i%) = v!(i%) / nr%
NEXT i%

'Compute the standard deviation of the voltage at each sample point
FOR i% = 0 TO ns% - 1
 sv!(i%) = 0!
FOR j% = 1 TO nr%
 sv!(i%) = sv!(i%) + (v2d!(i%, j%) - v!(i%)) ^ 2
NEXT j%
 sv!(i%) = SQR(sv!(i%) / nr%)
NEXT i%

'Log that a measurement was made
CALL coillogv(logfileP$, rawptr%)

'Save the samples in the raw data file
CALL coillogvsamp(rawfileP$, rawptr%, ns%, v!())

END SUB

SUB coillogf (logfile$, freq!)

'****************************************************************************
'This subroutine writes the coil rotation frequency to the log file.
'
'Input:
' logfile$, the name of the log file
' freq!, the measured coil rotation frequency
'
'Zachary Wolf
'10/3/94
'****************************************************************************

'Open the log file
filenum% = FREEFILE
OPEN logfile$ FOR APPEND AS filenum%

'Print the results
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Coil Rotation Frequency, Fcoil = "; freq!; " Hz"

'Close the log file
CLOSE filenum%

END SUB

SUB coillogfitv2 (logfile$, a!, b!, c!, phi!, s!)

'****************************************************************************
'This subroutine logs the parameters of the fit to the voltage samples.
'
'Input:
' logfile$, the name of the log file
' 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
'1/5/95
'****************************************************************************

'Open the log file
filenum% = FREEFILE
OPEN logfile$ FOR APPEND AS filenum%

'Write the fit parameters
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Coil Voltage: "
PRINT #filenum%, "Vfit(i) = ";
PRINT #filenum%, USING "##.#####"; c!;
PRINT #filenum%, " * cos(2*2pi*i/n + ";
PRINT #filenum%, USING "####.##"; phi!;
PRINT #filenum%, " deg)"

PRINT #filenum%, "        = ";
PRINT #filenum%, USING "##.#####"; a!;
PRINT #filenum%, " * sin(2*2pi*i/n) + ";
PRINT #filenum%, USING "##.#####"; b!;
PRINT #filenum%, " * cos(2*2pi*i/n)"

PRINT #filenum%, "RMS of the residuals = ";
PRINT #filenum%, USING "##.#####"; s!

'Close the log file
CLOSE filenum%

END SUB

SUB coillogplotv (logfile$, ns%, v!(), sv!(), vm!(), vp!(), vres!())

'****************************************************************************
'This subroutine writes voltage samples to the log file.
'
'Inputs:
' logfile$, the name of the log file
' ns%, the number of voltage samples
' v!(0 to ns%-1), the voltage samples
' sv!(0 to ns%-1), standard deviation on each voltage sample
' vm!(0 to ns%/2), magnitudes of the harmonics
' vp!(0 to ns%/2), phases of the harmonics in degrees
' vres!(0 to ns%/2), difference between the samples and the quad term
'
'Zachary Wolf
'1/7/95
'****************************************************************************

'Open the log file
filenum% = FREEFILE
OPEN logfile$ FOR APPEND AS filenum%

'Write the samples to the log file
FOR k% = 0 TO ns% - 1
 PRINT #filenum%, USING "###"; k%;
 PRINT #filenum%, USING " ###.####"; v!(k%);
 PRINT #filenum%, USING " ###.####"; sv!(k%);
 IF k% <= ns% / 2 THEN
  PRINT #filenum%, USING " ##.#####"; vm!(k%);
  PRINT #filenum%, USING " ####.###"; vp!(k%);
 ELSE
  PRINT #filenum%, USING " ##.#####"; 0!;
  PRINT #filenum%, USING " ##.#####"; 0!;
 END IF
 PRINT #filenum%, USING " #.######"; vres!(k%);
 PRINT #filenum%,
NEXT k%

'Close the log file
CLOSE filenum%

END SUB

SUB coillogv (logfile$, rawptr%)

'****************************************************************************
'This subroutine writes to the log file that a coil voltage measurement
'was made.
'
'Input:
' logfile$, the name of the log file
'
'Output:
' rawptr%, relates the line in the raw data file to the log file
'
'Zachary Wolf
'1/5/95
'****************************************************************************

'Relate log file to raw data file
STATIC index%
index% = index% + 1
rawptr% = index%

'Open the log file
filenum% = FREEFILE
OPEN logfile$ FOR APPEND AS filenum%

'Print the results
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Coil Voltage Measurement # "; rawptr%; " Made, Stored in RAWFILE"

'Close the log file
CLOSE filenum%

END SUB

SUB coillogvsamp (logfile$, rawptr%, ns%, v!())

'****************************************************************************
'This subroutine writes voltage samples to the log file.
'
'Inputs:
' logfile$, the name of the log file
' rawptr%, relates raw data file to log file
' ns%, the number of voltage samples
' v!(0 to ns%-1), the voltage samples
'
'Zachary Wolf
'1/7/95
'****************************************************************************

'Open the log file
filenum% = FREEFILE
OPEN logfile$ FOR APPEND AS filenum%

'Write the samples to the log file
PRINT #filenum%, USING "###"; rawptr%;
FOR i% = 0 TO ns% - 1
 PRINT #filenum%, USING " ###.####"; v!(i%);
NEXT i%
PRINT #filenum%,

'Close the log file
CLOSE filenum%

END SUB

SUB coilmeas (nh%, f!, vm!(), vp!())

'****************************************************************************
'This subroutine performs one coil measurement.
'
'Input:
' nh%, the number of harmonics to record
'
'Output:
' f!, the coil rotation frequency
' vm!(1 to nh%), voltage magnitude at each harmonic
' vp!(1 to nh%), voltage phase in degrees at each harmonic
'
'Zachary Wolf
'11/15/94, 1/1/95
'****************************************************************************

'Simplify the notation
ns% = nsampleP%

'Initialize data arrays
DIM v!(0 TO ns% - 1)                  'voltage samples
DIM sv!(0 TO ns% - 1)                 'standard deviation of v at each sample point
DIM vmfft!(0 TO ns% / 2)              'voltage magnitude at each harmonic
DIM vpfft!(0 TO ns% / 2)              'voltage phase at each harmonic

'Make sure there are trigger pulses for f1
CALL hp3457extsig(0)

'Perform the frequency and Vcoil measurements in succession

'Get the coil rotation frequency before the voltage measurement
CALL coilgetf(f1!)
SLEEP 1
CALL coilgetf(f2!)
'PRINT "Frequency = "; f1!; f2!

'Sample the voltage from the coil
CALL coilgetv(ns%, v!(), sv!())

'Get the coil rotation frequency again to form an average
CALL coilgetf(f3!)
SLEEP 1
CALL coilgetf(f4!)

'Average f
f! = (f1! + f2! + f3! + f4!) / 4!
'PRINT "Frequency = "; f!

'Fit a sin wave to v! for the log file
CALL coilfitv2(ns%, v!(), a!, b!, c!, phi!, s!)

'Calculate the harmonics in the voltage sample
CALL coilcalcvn(ns%, v!(), vmfft!(), vpfft!())

'Return the requested values
FOR i% = 1 TO nh%
 vm!(i%) = vmfft!(i%)
 vp!(i%) = vpfft!(i%)
NEXT i%

END SUB

SUB coilplotv

'Initialize a log file
logfileP$ = "vcoillog.dat"
filenum% = FREEFILE
OPEN logfileP$ FOR OUTPUT AS filenum%
CLOSE filenum%

'Initialize a raw data file
rawfileP$ = "vcoilraw.dat"
filenum% = FREEFILE
OPEN rawfileP$ FOR OUTPUT AS filenum%
CLOSE filenum%

'Simplify the notation
ns% = nsampleP%

'Initialize data arrays
DIM v!(0 TO ns% - 1)                  'voltage samples
DIM sv!(0 TO ns% - 1)                 'voltage samples
DIM vm!(0 TO ns% / 2)                 'voltage magnitude at each harmonic
DIM vp!(0 TO ns% / 2)                 'voltage phase at each harmonic
DIM vres!(0 TO ns% - 1)               'residuals after the fundamental is removed

'Sample the coil voltage
CALL coilgetv(ns%, v!(), sv!())

'Calculate the harmonics
CALL coilcalcvn(ns%, v!(), vm!(), vp!())

'Calculate the residuals
FOR i% = 0 TO ns% - 1
  vres!(i%) = v!(i%) - vm!(2) * COS(2 * 2 * 3.141592 * i% / ns% + vp!(2) * 3.141592 / 180!)
NEXT i%

'Initialize a plot file
pltfilecoilP$ = "vcoilplt.dat"
filenum% = FREEFILE
OPEN pltfilecoilP$ FOR OUTPUT AS filenum%
CLOSE filenum%

'Plot the results
CALL coillogplotv(pltfilecoilP$, ns%, v!(), sv!(), vm!(), vp!(), vres!())

END SUB