DECLARE SUB coilfitv1 (n%, y!(), a!, b!, c!, phi!, s!)
DECLARE SUB coillogfitv1 (logfile$, a!, b!, c!, phi!, s!)
DECLARE SUB coillogplotf (logfile$, nf%, tim$(), f!())
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$, logptr%)
DECLARE SUB coillogvsamp (logfile$, logptr%, 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 coilfitv1 (n%, y!(), a!, b!, c!, phi!, s!)

'****************************************************************************
'This subroutine fits n data points to a sine wave having 1 cycle in the
'n points.
'The fit is to y = a!*sin(1*2pi*i/n) + b!*cos(1*2pi*i/n)
'                = c!*cos(1*2pi*i/n + phi!)
'The rms of the residuals is also returned.
'y = c!*cos(1*2pi*i/n - nphi!)
'nphi! = -phi!
'y = c!*cos(1*2pi*i/n)*cos(nphi!) + c!*sin(1*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(1 * 2 * pi! * i% / n%) * SIN(1 * 2 * pi! * i% / n%)
 Scc! = Scc! + COS(1 * 2 * pi! * i% / n%) * COS(1 * 2 * pi! * i% / n%)
 Ssc! = Ssc! + SIN(1 * 2 * pi! * i% / n%) * COS(1 * 2 * pi! * i% / n%)
 Sys! = Sys! + y!(i%) * SIN(1 * 2 * pi! * i% / n%)
 Syc! = Syc! + y!(i%) * COS(1 * 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(1 * 2 * pi! * i% / n%) - b! * COS(1 * 2 * pi! * i% / n%)) ^ 2
NEXT i%
s! = SQR(s! / n%)

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

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$, logptr%)

'Save the values in the raw data file
CALL coillogvsamp(rawfileP$, logptr%, 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 coillogfitv1 (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(1*2pi*i/n + ";
PRINT #filenum%, USING "####.##"; phi!;
PRINT #filenum%, " deg)"

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

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

'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 coillogplotf (logfile$, nf%, tim$(), f!())

'****************************************************************************
'This subroutine writes coil frequency samples to the log file.
'
'Inputs:
' logfile$, the name of the log file
' nf%, the number of frequency samples
' tim$(1 to nf%), the time of each sample
' f!(1 to nf%), the coil frequency
'
'Zachary Wolf
'1/19/95
'****************************************************************************

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

'Write the samples to the log file
FOR k% = 1 TO nf%
 PRINT #filenum%, USING "###"; k%;
 PRINT #filenum%, USING " ###.#####"; f!(k%)
NEXT k%

'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 dipol term
'
'Zachary Wolf
'7/16/94, 1/5/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$, logptr%)

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

'Relate measurement to line in raw data file
STATIC index%
index% = index% + 1
logptr% = index%

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

'Print the results
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Coil Voltage Measurement Made, #"; logptr%; " in raw data file"

'Close the log file
CLOSE filenum%

END SUB

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

'****************************************************************************
'This subroutine writes voltage samples to the log file.
'
'Inputs:
' logfile$, the name of the log file
' logptr%, correlates voltage samples in raw file with log file
' ns%, the number of voltage samples
' v!(0 to ns%-1), the voltage samples
'
'Zachary Wolf
'1/10/95
'****************************************************************************

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

'Write the samples to the log file
PRINT #filenum%, USING "###"; logptr%;
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

'Perform the frequency and Vcoil measurements in succession

'Get the coil rotation frequency
CALL coilgetf(f1!)

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

'Get the coil rotation frequency again to form an average
CALL coilgetf(f2!)

'Average f
f! = (f1! + f2!) / 2!

'Fit a sin wave to v! for the log file
CALL coilfitv1(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 coilplotf

'Message
PRINT
PRINT "Beginning the coil frequency measurement cycle..."

'Get the measurement time
gettime:
INPUT "How long do you wish to remain at the desired current? (min) ", tmin!
tsec! = tmin! * 60

'Calculate the number of measurements which will be made
tbtwn% = 15             'make a measurement every 15 sec
nf% = 1 + tsec! / tbtwn%

'Initialize the data arrays
DIM tim$(1 TO nf%)
DIM f!(1 TO nf%)

'Record the coil frequency at regular intervals up to tsec
FOR i% = 1 TO nf%

'Message
PRINT
PRINT "Beginning a measurement cycle..."
PRINT "Time: "; TIME$

'Record the time of the measurement
tim$(i%) = TIME$
tm! = TIMER

'Measure the coil rotation frequency
CALL coilgetf(f!(i%))

'Write the results to the screen
PRINT
PRINT TIME$; "  Coil Rotation Frequency = "; f!(i%)

'Wait between measurements
WHILE TIMER - tm! < tbtwn%
WEND

'End of measurement loop
NEXT i%

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

'Write the results to the plt file
CALL coillogplotf(pltfilecoilP$, nf%, tim$(), f!())

END SUB

SUB coilplotv

'Set the parameters for the harmonics measurement
coilhpchanP% = coil1hpchanP%
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