DECLARE SUB calchar (ns%, v!(), vm!(), vp!())
DECLARE SUB getcoilfreq (c%, f!, sf!)
DECLARE SUB gettestparam ()
DECLARE SUB getv2 (c%, v2!, sv2!)
DECLARE SUB getvcoil (c%, ns%, v!(), sv!())
DECLARE SUB getvn (c%, nh%, vm!(), svm!(), vp!(), svp!())
DECLARE SUB getvtach (c%, ns%, v!(), sv!())
DECLARE SUB measglvsi ()
DECLARE SUB measharvsi ()
DECLARE SUB measrt ()
DECLARE SUB standardize ()
DECLARE SUB logcalc (logfile$)
DECLARE SUB logfreq (logfile$, freq!)
DECLARE SUB loggl (logfile$, rcoil!, nturns%, nrot%, nmeas%, imag!, freq!, sfreq!, v2!, sv2!, gl!, sgl!, tf!, stf!)
DECLARE SUB logglvsi (logfile$, rcoil!, nturns%, nrot%, nmeas%, ni%, imag!(), freq!(), v2!(), sv2!(), gl!(), sgl!(), tf!(), stf!())
DECLARE SUB loghar (logfile$, imag!, freq!, sfreq!, rcoil!, nturns%, nrot%, nmeas%, gl!, sgl!, nh%, vm!(), svm!(), vp!(), svp!(), bl!(), sbl!(), th!(), sth!(), r!(), sr!(), x!, sx!, y!, sy!)
DECLARE SUB logimag (logfile$, imag!)
DECLARE SUB lognewpage (logfile$)
DECLARE SUB logpltglvsi (logfile$, magname$, run$, ni%, imag!(), gl!(), sgl!(), tf!(), stf!())
DECLARE SUB logpltharvsi (logfile$, magname$, run$, nh%, ni%, inom!(), rvsi!(), srvsi!())
DECLARE SUB logrt (logfile$, nrt%, tim$(), imag!(), vmag!(), rmag!(), twin!(), twout!(), tiron!(), tcoil!(), tamb!())
DECLARE SUB logvsumm (logfile$, ns%, v!(), sv!(), a!, b!, c!, phi!, s!)
DECLARE SUB fileglvsi (ok$)
DECLARE SUB fileharvsi (ok$)
DECLARE SUB filert (ok$)
DECLARE SUB imagmeas (imag!)
DECLARE SUB imagramp (rr!, ifin!)
DECLARE SUB mathfitcos (n%, y!(), a!, b!, c!, phi!, s!)
DECLARE SUB gpibinit ()
DECLARE SUB dac488init ()
DECLARE SUB dac488setupv ()
DECLARE SUB hp3457cf (c%, freq!)
DECLARE SUB hp3457cnvtrig (c%, n%, v!())
DECLARE SUB hp3457init ()
DECLARE SUB hp3457preset ()
DECLARE SUB hp3457readtc (tc!)
DECLARE SUB hp3457readv (v!)
DECLARE SUB hp3457setuptc (c%)
DECLARE SUB hp3457setupv (c%)
DECLARE SUB k7011closecc (card%, chan%)
DECLARE SUB k7011init ()
DECLARE SUB k7011openall ()
DECLARE SUB FFTCalc (xreal!(), yimag!(), numdat%)

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

'****************************************************************************
'This module contains the main HER quadrupole measurement routines.
'
'Zachary Wolf
'9/10/94
'****************************************************************************

'Include the constants used in the module
REM $INCLUDE: 'param.inc'

SUB calchar (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 corrspeed (ns%, v!(), vc!())

'****************************************************************************
'This subroutine corrects the coil voltage samples to the values if the
'speed were constant at the average value.
'
'Input:
' ns%, the number of samples per revolution
' v!(0 to ns%-1), the measured voltage samples
'
'Output:
' vc!(0 to ns%-1), the corrected voltage samples
'
'Zachary Wolf
'9/7/94
'****************************************************************************

'Initialize
DIM vt!(0 TO ns% - 1)             'tachometer voltage
DIM svt!(0 TO ns% - 1)            'sigma vt
DIM corr!(0 TO ns% - 1)           'factor to correct for speed variations

'Get the voltage samples from the tachometer
CALL getvtach(tachhpchanP%, ns%, vt!(), svt!())

'Find the average tachometer voltage
avevt! = 0!
FOR i% = 0 TO ns% - 1
 avevt! = avevt! + vt!(i%)
NEXT i%
avevt! = avevt! / ns%

'Find the correction factor to take out speed variations
FOR i% = 0 TO ns% - 1
 corr!(i%) = avevt! / vt!(i%)
NEXT i%

'Apply the correction
FOR i% = 0 TO ns% - 1
 vc!(i%) = v!(i%) * corr!(i%)
NEXT i%

END SUB

SUB getcoilfreq (c%, f!, sf!)

'****************************************************************************
'This subroutine is used to get the average frequency and standard
'deviation of the rotation frequency of a measurement coil.
'
'Input:
' c%, HP3457 channel receiving ns% encoder pulses per revolution
'
'Output:
' f!, the average of the frequency measurements
' sf!, the standard deviation of the frequency measurements
'
'Zachary Wolf
'10/6/94
'****************************************************************************

'Simplify the notation
nf% = nfreqaveP%           '# measurements for the average
ns% = nsampleP%            '# encoder pulses per revolution

'Initialization
DIM fmeas!(1 TO nf%)

'Take frequency readings for the average
'Each measurement is the average of 10 readings
FOR i% = 1 TO nf%
 CALL hp3457cf(c%, fmeas!(i%))
 CALL logfreq(logfileP$, fmeas!(i%) / ns%)
NEXT i%

'Find the average frequency
f! = 0
FOR i% = 1 TO nf%
 f! = f! + fmeas!(i%)
NEXT i%
f! = f! / nf%

'Find the standard deviation
sf! = 0
FOR i% = 1 TO nf%
 sf! = sf! + (fmeas!(i%) - f!) ^ 2              'S
NEXT i%
sf! = sf! / nf%                                 'MS
sf! = SQR(sf!)                                  'RMS

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

END SUB

SUB gettestparam

'****************************************************************************
'This subroutine prompts the operator for the magnet name and other
'parameters of the test.
'
'Outputs put in the /testparam/ common block:
' magnameP$, name of magnet
' operatorP$, name of the operator doing the test
' teststandP$, the teststand being used
' coilP$, the ID of the measurement coil
' runP$, run number
' commentP$, comment about the test
'
'Zachary Wolf
'7/8/94
'****************************************************************************

begin:
'Get the name of the magnet
'PRINT "What is the name of the quadrupole being tested?"
'PRINT "(8 characters or less): ["; magnameP$; "]";
'INPUT magname$
'magname$ = UCASE$(magname$)
'IF magname$ <> "" THEN magnameP$ = magname$
PRINT "Beampipe Test"
magnameP$ = "beampipe"

'Get the operator's name(s)
PRINT "Enter the operator(s) name(s): ["; operatorP$; "]";
INPUT operator$
IF operator$ <> "" THEN operatorP$ = operator$

'Get the name of the test stand
'stand:
'PRINT "Enter which test stand is being used (Q1, Q2): ["; teststandP$; "]";
'INPUT teststand$
'teststand$ = UCASE$(teststand$)
'IF teststand$ <> "" THEN teststandP$ = teststand$
'IF teststandP$ <> "Q1" AND teststandP$ <> "Q2" THEN GOTO stand

'Get the name of the measurement coil
'meascoil:
'PRINT "Enter the name of the measurement coil (QC1, or QC2): ["; coilP$; "]";
'INPUT coil$
'coil$ = UCASE$(coil$)
'IF coil$ <> "" THEN coilP$ = coil$
'IF coilP$ <> "QC1" AND coilP$ <> "QC2" THEN GOTO meascoil

'Get the run number
PRINT "What is the run number (1, 2, 3, ...): ["; runP$; "]";
INPUT run$
IF run$ <> "" THEN runP$ = run$

'Get a comment
PRINT "Enter any comment about the run: ["; commentP$; "]";
LINE INPUT comment$
IF comment$ <> "" THEN commentP$ = comment$

'See if everything was entered correctly
PRINT
INPUT "Do you wish to make any changes to the values entered (Y or N): ", yn$
IF yn$ = "Y" OR yn$ = "y" THEN GOTO begin

END SUB

SUB getv2 (c%, v2!, sv2!)

'****************************************************************************
'This subroutine gets the voltage at the second rotation harmonic from
'the coil.
'
'Input:
' c%, DVM channel to sample
'
'Output:
' v2!, the voltage at the second harmonic
' sv2!, the standard deviation of v2!
'
'Zachary Wolf
'8/14/94
'****************************************************************************

'Simplify the notation
ns% = nsampleP%
nm% = nmeasaveP%

'Initialize data arrays
DIM v!(0 TO ns% - 1)                  'voltage samples
DIM sv!(0 TO ns% - 1)                 'standard deviation of v!(n%)
DIM vc!(0 TO ns% - 1)                 'corrected 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 v2nm!(1 TO nm%)                   'V2 for each measurement

'Perform the measurement nm% times
FOR i% = 1 TO nm%

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

'Apply a correction for speed variations
'CALL corrspeed(ns%, v!(), vc!())

'Calculate the harmonics in the corrected voltage sample
CALL calchar(ns%, v!(), vm!(), vp!())

'Save the value of V2
v2nm!(i%) = vm!(2)

'End loop over measurements
NEXT i%

'Find the average value of V2
v2! = 0!
FOR i% = 1 TO nm%
 v2! = v2! + v2nm!(i%)
NEXT i%
v2! = v2! / nm%

'Find the standard deviation of V2
sv2! = 0!
FOR i% = 1 TO nm%
 sv2! = sv2! + (v2nm!(i%) - v2!) ^ 2
NEXT i%
sv2! = SQR(sv2! / nm%)

END SUB

SUB getvcoil (c%, ns%, v!(), sv!())

'****************************************************************************
'This subroutine gets the triggered voltages from the coil.
'
'Input:
' c%, DVM channel to sample
' 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
'****************************************************************************

'Simplify the notation
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%

'Convention
'The convention is that the coil voltage is positive when the coil turns
'in front of a South pole.  It does not do this.  Rather than reverse
'the leads, I fix this in software:
FOR i% = 0 TO ns% - 1
 v!(i%) = -v!(i%)
NEXT i%

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

'Log the results
CALL logvsumm(logfileP$, ns%, v!(), sv!(), a!, b!, c!, phi!, s!)

END SUB

SUB getvn (c%, nh%, vm!(), svm!(), vp!(), svp!())

'****************************************************************************
'This subroutine gets the voltage harmonics from the coil.
'
'Input:
' c%, DVM channel to sample
' nh%, the number of harmonics to record
'
'Output:
' vm!(1 to nh%), average voltage magnitude at each harmonic
' svm!(1 to nh%), standard deviation of the voltage magnitude measurements
' vp!(1 to nh%), average voltage phase in degrees at each harmonic
' svp!(1 to nh%), standard deviation of the voltage phase measurements
'
'Zachary Wolf
'7/13/94, 8/14/94
'****************************************************************************

'Simplify the notation
ns% = nsampleP%
nm% = nmeasaveP%

'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 vc!(0 TO ns% - 1)                 'corrected voltage samples
DIM vmfft!(0 TO ns% / 2)              'voltage magnitude at each harmonic
DIM vpfft!(0 TO ns% / 2)              'voltage phase at each harmonic
DIM mstore!(0 TO ns% / 2, 1 TO nm%)   'store magnitudes for average
DIM pstore!(0 TO ns% / 2, 1 TO nm%)   'store phases for average

'Perform the measurement nm% times
FOR j% = 1 TO nm%

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

'Apply a correction for speed variations
'CALL corrspeed(ns%, v!(), vc!())

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

'Save the values for averaging
FOR i% = 0 TO ns% / 2
 mstore!(i%, j%) = vmfft!(i%)
 pstore!(i%, j%) = vpfft!(i%)
NEXT i%

'End loop over measurements
NEXT j%

'Compute the average magnitudes and phases
FOR i% = 1 TO nh%
 vm!(i%) = 0!
 vp!(i%) = 0!
FOR j% = 1 TO nm%
 vm!(i%) = vm!(i%) + mstore!(i%, j%)
 vp!(i%) = vp!(i%) + pstore!(i%, j%)
NEXT j%
 vm!(i%) = vm!(i%) / nm%
 vp!(i%) = vp!(i%) / nm%
NEXT i%

'Compute the standard deviations of the magnitudes and phases
FOR i% = 1 TO nh%
 svm!(i%) = 0!
 svp!(i%) = 0!
FOR j% = 1 TO nm%
 svm!(i%) = svm!(i%) + (mstore!(i%, j%) - vm!(i%)) ^ 2
 svp!(i%) = svp!(i%) + (pstore!(i%, j%) - vp!(i%)) ^ 2
NEXT j%
 svm!(i%) = SQR(svm!(i%) / nm%)
 svp!(i%) = SQR(svp!(i%) / nm%)
NEXT i%

END SUB

SUB getvtach (c%, ns%, v!(), sv!())

'****************************************************************************
'This subroutine gets the triggered voltages from the tachometer.
'These readings allow a correction for speed variations.
'
'Input:
' c%, DVM channel to sample
' 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
'9/7/94
'****************************************************************************

'Simplify the notation
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%

END SUB

SUB main

'****************************************************************************
'This is the main HER quadrupole routine.  It calls all the subroutines
'which do the testing.
'
'Zachary Wolf
'9/10/94
'****************************************************************************

'Prepare the screen
CLS

'Initialize the GPIB
CALL gpibinit

'Initialize the DAC
CALL dac488init
CALL dac488setupv

'Initialize the HP3457
CALL hp3457init

'Initialize the multiplexer
CALL k7011init

'Print a message about the program
PRINT
PRINT "Program BEAMPIPE"
PRINT "This program performs Uli's beampipe test."
PRINT

'Get test parameters from the operator
testparam:
CALL gettestparam

'Initialize the log file for the RT measurements
rt:
CALL filert(ok$)
IF ok$ <> "y" THEN
 PRINT
 PRINT "Problem initializing the data files."
 PRINT "Make sure you are using a new run number."
 GOTO testparam:
END IF

'Perform the RT measurements
CALL measrt
GOTO done

'Initialize the log, data, and plot files for the GL vs I measurement
gl:
PRINT
CALL fileglvsi(ok$)
IF ok$ <> "y" THEN
 PRINT
 PRINT "Problem initializing the harmonics data files."
 PRINT "This error should never happen."
 GOTO done
END IF

'Standardize the magnet
CALL standardize

'Measure GL vs I
CALL measglvsi

'Initialize the log, data, and plot files for the harmonics measurement
har:
PRINT
CALL fileharvsi(ok$)
IF ok$ <> "y" THEN
 PRINT
 PRINT "Problem initializing the harmonics data files."
 PRINT "This error should never happen."
 GOTO done
END IF

'Measure the harmonics
CALL measharvsi

'Message
done:
PRINT
PRINT "Ramping to 0 Amps."

'Ramp to 0 amps
CALL imagramp(imagramprateP!, 0!)

'Message
PRINT "The test is complete."

END SUB

SUB measglvsi

'****************************************************************************
'This subroutine measures the integrated gradient as a function of current
'for the B-Factory HER quads.
'Constants come from param.inc.
'
'Zachary Wolf
'7/8/94, 8/14/94
'****************************************************************************

'Message
PRINT
PRINT "Beginning the GL vs I measurement cycle..."

'Simplify the notation
pi! = 3.1415926#
ni% = imagnglvsiP%
ns% = nsampleP%
nturns% = nturnscoilP%
c% = coilhpchanP%
cf% = freqhpchanP%
rcoil! = rcoilP!

'Initialize data arrays
DIM imag!(1 TO ni%)               'measured magnet current
DIM freq!(1 TO ni%)               'coil rotation frequency
DIM sfreq!(1 TO ni%)              'standard deviation of coil rotation frequency
DIM v2!(1 TO ni%)                 'average voltage magnitude V2
DIM sv2!(1 TO ni%)                'standard deviation of V2 values
DIM gl!(1 TO ni%)                 'average integrated gradient at each current
DIM sgl!(1 TO ni%)                'standard deviation of gl
DIM tf!(1 TO ni%)                 'average transfer function
DIM stf!(1 TO ni%)                'standard deviation of tf

'Loop over currents
FOR i% = 1 TO ni%

'Message
PRINT
PRINT "Ramping to "; imagglvsiP!(i%); " Amps..."

'Ramp to the desired current
CALL imagramp(imagramprateP!, imagglvsiP!(i%))

'Make an initial measurement of the magnet current
CALL imagmeas(iinit!)

'Log the magnet current
CALL logimag(logfileP$, iinit!)

'Message
PRINT "Measuring the coil voltage..."

'Get the voltage at the second harmonic
CALL getv2(c%, v2!(i%), sv2!(i%))

'Get the coil rotation frequency
CALL getcoilfreq(cf%, freq!(i%), sfreq!(i%))

'Repeat the magnet current measurement and use this value in calculations
CALL imagmeas(imag!(i%))

'Some Math
'Vn = Nturns * velocity * Bn * L
'   = Nturns * Rcoil * ang_freq * (BLn)
'BLn = Vn / (Nturns * Rcoil * ang_freq)
'GL = BL2 / Rcoil = V2 / (Nturns * Rcoil^2 * ang_freq)

'Compute the integrated gradient
gl!(i%) = v2!(i%) / (nturns% * rcoil! ^ 2 * 2! * pi! * freq!(i%))
siggl! = (sv2!(i%) / v2!(i%)) ^ 2 + (sfreq!(i%) / freq!(i%)) ^ 2
siggl! = SQR(siggl!)
siggl! = gl!(i%) * siggl!
sgl!(i%) = siggl!

'Compute the transfer function
tf!(i%) = gl!(i%) / imag!(i%)
stf!(i%) = sgl!(i%) / imag!(i%)

'Print the results to the screen
PRINT
PRINT "Integrated Gradient Measurement:"
PRINT "Imag = "; imag!(i%); " A"
'PRINT "Coil Radius = "; rcoil!; " m"
'PRINT "Coil Turns = "; nturns%
PRINT "Coil Rotation Frequency = "; freq!(i%); " +- "; sfreq!(i%); " 1/s"
PRINT "Coil Voltage V2 = "; v2!(i%); " +- "; sv2!(i%); " V"
PRINT "Integrated Gradient = "; gl!(i%); " +- "; sgl!(i%); " T"
PRINT "Transfer Function = "; tf!(i%); " +- "; stf!(i%); " T/A"

'Write these results to the log file
CALL loggl(logfileP$, rcoil!, nturns%, nrotaveP%, nmeasaveP%, imag!(i%), freq!(i%), sfreq!(i%), v2!(i%), sv2!(i%), gl!(i%), sgl!(i%), tf!(i%), stf!(i%))

'End loop over currents
NEXT i%

'Write a summary of the results to the dat file and to the log file
CALL logglvsi(datfileP$, rcoil!, nturns%, nrotaveP%, nmeasaveP%, ni%, imag!(), freq!(), v2!(), sv2!(), gl!(), sgl!(), tf!(), stf!())
CALL logglvsi(logfileP$, rcoil!, nturns%, nrotaveP%, nmeasaveP%, ni%, imag!(), freq!(), v2!(), sv2!(), gl!(), sgl!(), tf!(), stf!())

'Write the results to the plotting file
CALL logpltglvsi(pltfileP$, magnameP$, runP$, ni%, imag!(), gl!(), sgl!(), tf!(), stf!())

'Write a summary of the calculations to the log file
CALL lognewpage(logfileP$)
CALL logcalc(logfileP$)

END SUB

SUB measharvsi

'****************************************************************************
'This subroutine controls the harmonics measurements of the B-Factory
'HER quads.
'Constants come from param.inc.
'
'Zachary Wolf
'7/8/94, 8/14/94
'****************************************************************************

'Message
PRINT
PRINT "Beginning the harmonics measurement cycle..."

'Simplify the notation
pi! = 3.1415926#
ni% = imagnharvsiP%
nh% = nhardisplayP%
ns% = nsampleP%
nturns% = nturnscoilP%
c% = coilhpchanP%
cf% = freqhpchanP%
rcoil! = rcoilP!

'Initialize data arrays
DIM vm!(1 TO nh%)                 'average voltage magnitudes
DIM svm!(1 TO nh%)                'standard deviation of voltage magnitudes
DIM vp!(1 TO nh%)                 'average voltage phases in degrees
DIM svp!(1 TO nh%)                'standard deviation of voltage phases
DIM bl!(1 TO nh%)                 'field strength at Rcoil of the n'th harmonic
DIM sbl!(1 TO nh%)                'standard deviation of Bn
DIM th!(1 TO nh%)                 'field angle
DIM sth!(1 TO nh%)                'standard deviation of THn
DIM r!(1 TO nh%)                  'BLn/BL2
DIM sr!(1 TO nh%)                 'standard deviation of r
DIM rvsi!(1 TO nh%, 1 TO ni%)     'BLn/BL2 vs I
DIM srvsi!(1 TO nh%, 1 TO ni%)    'standard deviation of BLn/BL2 vs I

'Loop over currents
FOR i% = 1 TO ni%

'Message
PRINT
PRINT "Ramping to "; imagharvsiP!(i%); " Amps..."

'Ramp to the desired current
CALL imagramp(imagramprateP!, imagharvsiP!(i%))

'Make an initial measurement of the magnet current
CALL imagmeas(iinit!)

'Log the magnet current
CALL logimag(logfileP$, iinit!)

'Message
PRINT "Measuring the coil voltage..."

'Get the voltage harmonics
CALL getvn(c%, nh%, vm!(), svm!(), vp!(), svp!())

'Get the coil rotation frequency
CALL getcoilfreq(cf%, freq!, sfreq!)

'Repeat the magnet current measurement and use this value in calculations
CALL imagmeas(imag!)

'Some Math
'Vn = Nturns * velocity * Bn * L
'   = Nturns * Rcoil * ang_freq * (BLn)
'BLn = Vn / (Nturns * Rcoil * ang_freq)
'GL = BL2 / Rcoil = V2 / (Nturns * Rcoil^2 * ang_freq)
'Rn = BLn / BL2 = Vn / V2

'Compute the strength of each field harmonic
FOR j% = 1 TO nh%
 bl!(j%) = vm!(j%) / (nturns% * rcoil! * 2! * pi! * freq!)
 sigbl! = (svm!(j%) / vm!(j%)) ^ 2 + (sfreq! / freq!) ^ 2
 sigbl! = SQR(sigbl!)
 sigbl! = bl!(j%) * sigbl!
 sbl!(j%) = sigbl!
NEXT j%

'Compute the angle of the first south pole wrt horizontal for each harmonic
FOR j% = 1 TO nh%
 th!(j%) = -vp!(j%) / j%
 sth!(j%) = svp!(j%) / j%
NEXT j%

'Compute the field strength ratios at the coil radius
FOR j% = 1 TO nh%
 r!(j%) = vm!(j%) / vm!(2)
 sr!(j%) = r!(j%) * SQR((svm!(j%) / vm!(j%)) ^ 2 + (svm!(2) / vm!(2)) ^ 2)
NEXT j%
sr!(2) = 0     'ratio is 1, sigma=0

'Compute the integrated gradient
gl! = bl!(2) / rcoil!
sgl! = sbl!(2) / rcoil!

'Compute the transfer function
tf! = gl! / imag!
stf! = sgl! / imag!

'Compute the magnetic center relative to the coil
'For a normal quad:
x! = -(1! / gl!) * bl!(1) * SIN(th!(1) * pi! / 180!)
y! = -(1! / gl!) * bl!(1) * COS(th!(1) * pi! / 180!)
'With the leads reversed, the signs change
IF th!(2) < 0 THEN
  x! = -x!
  y! = -y!
END IF

'Compute the error on the magnetic center
sx! = (sgl! / gl!) ^ 2 + (sbl!(1) / bl!(1)) ^ 2 + (COS(th!(1) * pi! / 180!) / (SIN(th!(1) * pi! / 180!) + .000001)) ^ 2 * (sth!(1) * pi! / 180!) ^ 2
sx! = SQR(sx!)
sx! = ABS(x!) * sx!
sy! = (sgl! / gl!) ^ 2 + (sbl!(1) / bl!(1)) ^ 2 + (SIN(th!(1) * pi! / 180!) / (COS(th!(1) * pi! / 180!) + .000001)) ^ 2 * (sth!(1) * pi! / 180!) ^ 2
sy! = SQR(sy!)
sy! = ABS(y!) * sy!

'Print the results to the screen
PRINT
PRINT "Harmonics:"
PRINT "Imag = "; imag!; " A"
'PRINT "Coil Radius = "; rcoil!; " m"
'PRINT "Number of Turns = "; nturns%
PRINT "Coil Rotation Frequency = "; freq!; " +- "; sfreq!; " 1/s"
PRINT "Coil Voltage V2 = "; vm!(2); " +- "; svm!(2); " V"
PRINT "Integrated Gradient = "; gl!; " +- "; sgl!; " T"
PRINT "Transfer Function = "; tf!; " +- "; stf!; " T/A"
PRINT
PRINT " N    |V|   sig|V|  V Phase  sigVPh   BLn   THspole BLn/BL2   sig  "
PRINT "      (V)     (V)    (deg)   (deg)    (Tm)   (deg)    (%)   BLn/BL2"
PRINT "--- -------+------- -------+------- ------- ------- -------+-------"
IF nh% > 10 THEN nprint% = 10 ELSE nprint% = nh%
FOR j% = 1 TO nprint%
 PRINT USING "###"; j%;
 PRINT USING " ##.####"; vm!(j%);
 PRINT USING " ##.####"; svm!(j%);
 PRINT USING " ####.##"; vp!(j%);
 PRINT USING " ####.##"; svp!(j%);
 PRINT USING " ###.###"; bl!(j%);
 PRINT USING " ####.##"; th!(j%);
 PRINT USING " ###.###"; r!(j%) * 100!;
 PRINT USING " ###.###"; sr!(j%) * 100!
NEXT j%

'Write the results to the dat file
IF i% > 1 THEN CALL lognewpage(datfileP$)
CALL loghar(datfileP$, imag!, freq!, sfreq!, rcoil!, nturns%, nrotaveP%, nmeasaveP%, gl!, sgl!, nh%, vm!(), svm!(), vp!(), svp!(), bl!(), sbl!(), th!(), sth!(), r!(), sr!(), x!, sx!, y!, sy!)

'Write the results to the log file
CALL loghar(logfileP$, imag!, freq!, sfreq!, rcoil!, nturns%, nrotaveP%, nmeasaveP%, gl!, sgl!, nh%, vm!(), svm!(), vp!(), svp!(), bl!(), sbl!(), th!(), sth!(), r!(), sr!(), x!, sx!, y!, sy!)

'Save the results for the plots
FOR j% = 1 TO nh%
 rvsi!(j%, i%) = r!(j%)
 srvsi!(j%, i%) = sr!(j%)
NEXT j%

'End loop over currents
NEXT i%

'Plot the results
CALL logpltharvsi(pltfileP$, magnameP$, runP$, nh%, ni%, imagharvsiP!(), rvsi!(), srvsi!())

'Write a summary of the calculations to the log file
CALL lognewpage(logfileP$)
CALL logcalc(logfileP$)

END SUB

SUB measrt

'****************************************************************************
'This subroutine controls the resistance and temperature measurements.
'
'Zachary Wolf
'9/6/94
'****************************************************************************

'Message
PRINT
PRINT "Beginning the magnet resistance and temperature measurement cycle..."

'Get the desired current
getinom:
PRINT
INPUT "What current do you wish to ramp to? (0<I<2000 A) ", inom!
IF inom! < 0 OR inom! > 2000 THEN
 PRINT : PRINT "The current must be between 0<I<2000 A."
 GOTO getinom:
END IF

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

'Ramp the magnet
PRINT "Ramping to "; inom!; " Amps..."
CALL imagramp(imagramprateP!, inom!)

'Calculate the number of measurements which will be made
nrt% = tsec! / tbtwnmeassecP%

'Initialize the data arrays
DIM tim$(1 TO nrt%)
DIM imag!(1 TO nrt%)
DIM vmag!(1 TO nrt%)
DIM rmag!(1 TO nrt%)
DIM twin!(1 TO nrt%)
DIM twout!(1 TO nrt%)
DIM tiron!(1 TO nrt%)
DIM tcoil!(1 TO nrt%)
DIM tamb!(1 TO nrt%)

'Record the temperatures and resistance at regular intervals up to tsec
FOR i% = 1 TO nrt%

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

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

'TEMPERATURE MEASUREMENTS

'Set up the HP3457 to read the temperatures
CALL hp3457setuptc(mxhpchanP%)

'Measure the water in temperature
CALL k7011closecc(mxcardP%, twinmxchanP%)     'close the multiplexer channel
CALL hp3457readtc(twin!(i%))                  'measure the temperature
CALL k7011openall                             'open all multiplexer channels

'Measure the water out temperature
CALL k7011closecc(mxcardP%, twoutmxchanP%)
CALL hp3457readtc(twout!(i%))
CALL k7011openall

'Measure the temperature of the magnet iron
CALL k7011closecc(mxcardP%, tironmxchanP%)
CALL hp3457readtc(tiron!(i%))
CALL k7011openall

'Measure the temperature of a magnet coil
CALL k7011closecc(mxcardP%, tcoilmxchanP%)
CALL hp3457readtc(tcoil!(i%))
CALL k7011openall

'Measure the ambient temperature
CALL k7011closecc(mxcardP%, tambmxchanP%)
CALL hp3457readtc(tamb!(i%))
CALL k7011openall

'Put the DVM back in a high impedance state
CALL hp3457preset

'MAGNET RESISTANCE MEASUREMENTS

'Measure and store the magnet current
CALL imagmeas(imag!(i%))

'Set up the DVM to measure the magnet voltage
CALL hp3457setupv(mxhpchanP%)

'Measure the magnet voltage
CALL k7011closecc(mxcardP%, vmagmxchanP%)   'close the multiplexer channel
CALL hp3457readv(vmag!(i%))                 'measure the magnet voltage
CALL k7011openall                           'open all multiplexer channels

'Compute the magnet resistance
rmag!(i%) = vmag!(i%) / imag!(i%)

'Write the results to the screen
PRINT
PRINT "Resistance and Temperature Measurements"
PRINT " Imag    Vmag    Rmag     T1      T2      T3      T4    Ambient"
PRINT "  (A)     (V)   (mohm)  (deg C) (deg C) (deg C) (deg C) (deg C)"
PRINT "------- ------- ------- ------- ------- ------- ------- -------"
PRINT USING "####.##"; imag!(i%);
PRINT USING " ###.###"; vmag!(i%);
PRINT USING " ###.###"; rmag!(i%) * 1000!;
PRINT USING " ####.##"; twin!(i%);
PRINT USING " ####.##"; twout!(i%);
PRINT USING " ####.##"; tiron!(i%);
PRINT USING " ####.##"; tcoil!(i%);
PRINT USING " ####.##"; tamb!(i%)

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

'End of measurement loop
NEXT i%

'Write the results to the log file
CALL logrt(logfileP$, nrt%, tim$(), imag!(), vmag!(), rmag!(), twin!(), twout!(), tiron!(), tcoil!(), tamb!())

'Ramp to 0 amps
CALL imagramp(imagramprateP!, 0!)

END SUB

SUB standardize

'****************************************************************************
'This procedure controls the standardization of the magnet.
'Parameters come from param.inc.
'
'Zachary Wolf
'4/23/94
'****************************************************************************

'Let the operator know what is happening
PRINT
PRINT "Standardizing the magnet..."

'Simplify the notation
ni% = imagnstandP%

'Loop over currents
FOR i% = 1 TO ni%

'Message
PRINT
PRINT "Ramping to "; imagstandP!(i%); " Amps."

'Ramp to the desired current
CALL imagramp(imagramprateP!, imagstandP!(i%))

'Measure the magnet current
CALL imagmeas(imag!)

'Log the magnet current
CALL logimag(logfileP$, imag!)

'End loop over currents
NEXT i%

END SUB