DECLARE SUB gethar ()
DECLARE SUB meashar ()
DECLARE SUB calchar (ns%, v!(), vm!(), vp!())
DECLARE SUB getcoilfreq (c%, f!, sf!)
DECLARE SUB gettestparam ()
DECLARE SUB getvcoil (c%, ns%, v!(), sv!())
DECLARE SUB getvn (c%, nh%, vm!(), svm!(), vp!(), svp!())
DECLARE SUB logcalc (logfile$)
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 logpltharvsi (logfile$, magname$, run$, nh%, ni%, inom!(), rvsi!(), srvsi!())
DECLARE SUB logv (logfile$, ns%, v!(), vm!(), vp!())
DECLARE SUB fileharvsi (ok$)
DECLARE SUB filemkdir (destdir$)
DECLARE SUB imagmeas (imag!)
DECLARE SUB dac488init ()
DECLARE SUB dac488setupv ()
DECLARE SUB gpibinit ()
DECLARE SUB hp3457cf (c%, freq!)
DECLARE SUB hp3457init ()

'****************************************************************************
'Program TEST
'This program is used to make the rotating coil measurements of the
'B-Factory HER quadrupoles.  It samples the voltage from the rotating
'coil and Fourier analyzes the samples to get the strengths of
'each component of the field.
'
'Zachary Wolf
'7/16/94
'****************************************************************************

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

'Prepare the screen
CLS

'Initialize the GPIB
CALL gpibinit

'Initialize the DAC
CALL dac488init
CALL dac488setupv

'Initialize the HP3457
CALL hp3457init

'Initialize the DACHR
'CALL dac488hrinit
'CALL dac488hrsetupv
'STOP

'Open a log and plot file
'CALL filetest

'Make coil voltage measurements
'CALL measvtest

'Measure the harmonics
CALL meashar

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

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

'Message
PRINT "The test is complete."

END

SUB filetest

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

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

END SUB

SUB gethar

'****************************************************************************
'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

'There is only one current for this test
i% = 1

'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%

'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 meashar

'Get test parameters from the operator
testparam:
CALL gettestparam

'Create the \magdata subdirectory to store the data and log files
subdir$ = "c:\magdata\" + magnameP$
CALL filemkdir(subdir$)

'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 data files."
 PRINT "Make sure you are using a new run number."
 GOTO testparam:
END IF

'Measure the harmonics
CALL gethar

END SUB

SUB measvtest

'Parameters
ns% = nsampleP%
c% = coilhpchanP%
cf% = freqhpchanP%

'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

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

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

'Get the measurement coil rotation frequency
CALL hp3457cf(cf%, freq!)
freq! = freq! / ns%                'take out the ns% pulses per revolution
PRINT "cf = "; cf%
PRINT "freq! = "; freq!

'Plot the results
CALL logv(pltfileP$, ns%, v!(), vm!(), vp!())

END SUB