DECLARE SUB imaggetiav (imag!, simag!)
DECLARE SUB imagramp (inom!)
DECLARE SUB blngetblnav (nh%, bln!(), sbln!(), thn!(), sthn!())
DECLARE SUB blnlogcalc (logfile$)
DECLARE SUB harcalcratio (nh%, bln!(), sbln!(), nmain%, rn!(), srn!())
DECLARE SUB hardathar (datfile$, imag!, simag!, nmain%, nh%, bln!(), sbln!(), thn!(), sthn!(), r!(), sr!())
DECLARE SUB harlogcalc (logfile$)
DECLARE SUB harplthar (pltfile$, legend$, nh%, rn!(), srn!())
DECLARE SUB spolcalcsl (rcoil!, bl3!, sbl3!, sl!, ssl!)
DECLARE SUB spolcalcxy (rcoil!, nh%, bln!(), sbln!(), thn!(), sthn!(), x!, sx!, y!, sy!)
DECLARE SUB spoldatsl (datfile$, imag!, simag!, sl!, ssl!, tf!, stf!)
DECLARE SUB spollogcalc (logfile$)
DECLARE SUB spollogsl (logfile$, bl3!, sbl3!, sl!, ssl!)
DECLARE SUB spollogxy (logfile$, thn!(), x!, sx!, y!, sy!)
DECLARE SUB spolpltsl (pltfile$, imag!, sl!, ssl!, tf!, stf!)
DECLARE SUB filenewpage (filename$)

'****************************************************************************
'Module SEXTUPOL.BAS
'
'These subroutines perform calculations related to sextupole magnet
'measurements.
'
'Zachary Wolf
'1/2/95
'****************************************************************************

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

'The required parameters are
'coilradiusmP!
'coilnturnsP%
'coilnrotaveP%
'imagtestP!()
'imagharflagP!()
'nimagtestP%
'logfileP$
'sldatP$
'slpltP$
'hardatP$
'harpltP$
'magnameP$
'runP$
'blnnmeasaveP%
'spolnhardisplayP%
'
'The required parameters are in
'COIL
'IMAG
'FILE
'SPOL
'GETPARAM
'BLN
'
'Sample SPOL parameters
'SPOL
'CONST spolnhardisplayP% = 16      'the number of harmonics to print out < nsampleP%/2

SUB spolcalcsl (rcoil!, bl3!, sbl3!, sl!, ssl!)

'****************************************************************************
'This subroutine calculates the integrated sextupole strength from the
'integrated field strength at the coil radius.
'
'Input:
' rcoil!, the coil radius
' bl3!, the sextupole integrated field strength at the coil radius (Tm)
' sbl3!, rms variation of bl3 measurements
'
'Output:
' sl!, the integrated sextupole strength (T/m)
' ssl!, rms variation of sl measurements
'
'Zachary Wolf
'5/5/95
'****************************************************************************

'Calculation Of SL:
'The sextupole strength S is defined by the vertical
'field on the x-axis (median plane).
'It is the quadratic term in the Taylor expansion.
' By(x) = 1/2 S x^2
'On the x-axis, By(x) = Bth(r=x,th=0)
'From the expression for Bth (see blnlogcalc),
' By(x) = - B3 (x/rref)^2 sin(-3 thspole)
'Take rref = Rcoil, thspole = pi/6, then at x=Rcoil,
' By(Rcoil) = B3 = 1/2 S Rcoil^2
'So,
' S = 2 B3 / Rcoil^2
'The integrated sextupole strength is
' SL = 2 BL3 / Rcoil^2

'Compute SL
sl! = 2 * bl3! / rcoil! ^ 2
ssl! = 2 * sbl3! / rcoil! ^ 2

'Print the result to the screen
PRINT "Integrated Field Strength BL(3) = "; bl3!; " +- "; sbl3!; " Tm"
PRINT "Integrated Sextupole Strength SL = "; sl!; " +- "; ssl!; " T/m"

'Log the result
CALL spollogsl(logfileP$, bl3!, sbl3!, sl!, ssl!)

END SUB

SUB spolcalcxy (rcoil!, nh%, bln!(), sbln!(), thn!(), sthn!(), x!, sx!, y!, sy!)

'****************************************************************************
'This subroutine calculates the (x, y) position of the sextupole magnetic
'center relative to the coil axis.  X is horizontal to the right (viewed
'from the encoder end) and Y is up.
'
'Input:
' rcoil!, the coil radius
' nh%, the number of voltage harmonics recorded
' bln!(1 to nh%), the integrated strength of each harmonic at the coil radius (T)
' sbln!(1 to nh%), rms variation of bl
' thn!(1 to nh%), the angle of the first south pole of each harmonic wrt horizontal (deg)
' sthn!(1 to nh%), rms variation of th
'
'Output:
' x!, x coordinate of the magnetic center wrt the coil axis (m)
' sx!, calculated error for x
' y!, y coordinate of the magnetic center wrt the coil axis (m)
' sy!, calculated error for y
'
'Zachary Wolf
'1/2/94
'****************************************************************************

'Simplify the notation
pi! = 3.1415926#

'Compute SL from bl(3)
sl! = 2 * bln!(3) / rcoil! ^ 2
ssl! = 2 * sbln!(3) / rcoil! ^ 2

'Compute the magnetic center relative to the coil
'For a normal sextupole (th(3) = 30 deg):
x! = -(bln!(2) / (sl! * rcoil!)) * SIN(2 * thn!(2) * pi! / 180!)
y! = -(bln!(2) / (sl! * rcoil!)) * COS(2 * thn!(2) * pi! / 180!)
'With the leads reversed, the signs change
IF thn!(3) < 0 THEN
  x! = -x!
  y! = -y!
END IF

'Compute the error on the magnetic center
sx! = ABS(x!) * SQR((ssl! / sl!) ^ 2 + (sbln!(2) / bln!(2)) ^ 2 + (COS(2 * thn!(2) * pi! / 180!) / (SIN(2 * thn!(2) * pi! / 180!) + .000001)) ^ 2 * (2 * sthn!(2) * pi! / 180!) ^ 2)
sy! = ABS(y!) * SQR((ssl! / sl!) ^ 2 + (sbln!(2) / bln!(2)) ^ 2 + (SIN(2 * thn!(2) * pi! / 180!) / (COS(2 * thn!(2) * pi! / 180!) + .000001)) ^ 2 * (2 * sthn!(2) * pi! / 180!) ^ 2)

END SUB

SUB spoldatsl (datfile$, imag!, simag!, sl!, ssl!, tf!, stf!)

'****************************************************************************
'This subroutine writes an SL measurement at Imag to the data file.
'
'Input:
' datfile$, the name of the data file
' imag!, the measured current
' simag!, the standard deviation of imag
' sl!, the integrated sextupole strength at the measured current
' ssl!, the standard deviation of SL
' tf!, SL/I
' stf!, the standard deviation of TF
'
'Zachary Wolf
'8/14/94, 1/4/95
'****************************************************************************

'Keep track of Ncall so a header can be written on the first call
STATIC ncall%
ncall% = ncall% + 1

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

'On the first call, write a header for the hysteresis loop
IF ncall% = 1 THEN
 PRINT #filenum%,
 PRINT #filenum%,
 PRINT #filenum%,
 PRINT #filenum%, "            INTEGRATED SEXTUPOLE STRENGTH VS CURRENT"
 PRINT #filenum%,
 PRINT #filenum%, "Coil: Radius = "; coilradiusmP!; " m,"; " # Turns = "; coilnturnsP%
 PRINT #filenum%, "Average: # Rotations/Measurement = "; coilnrotaveP%; ", # Measurements = "; blnnmeasaveP%
 PRINT #filenum%,
 PRINT #filenum%, "  Imag    sigImag     SL     sigSL      SL/I   sigSL/I "
 PRINT #filenum%, "  (A)       (A)     (T/m)    (T/m)    (T/m/kA) (T/m/kA)"
 PRINT #filenum%, "--------+--------  --------+--------  --------+--------"
END IF

'Write the result
PRINT #filenum%, USING "####.###"; imag!;
PRINT #filenum%, USING " ####.###"; simag!;
PRINT #filenum%, USING "  ##.#####"; sl!;
PRINT #filenum%, USING " ##.#####"; ssl!;
PRINT #filenum%, USING "  ###.####"; tf! * 1000!;
PRINT #filenum%, USING " ###.####"; stf! * 1000!

'Close the log file
CLOSE filenum%

END SUB

SUB spollogcalc (logfile$)

'****************************************************************************
'This subroutine summarizes the calculations used in the program.
'
'Inputs:
' logfile$, the name of the log file
'
'Zachary Wolf
'9/22/94, 12/3/94, 5/4/95
'****************************************************************************

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

'SL calculations
PRINT #filenum%,
PRINT #filenum%, "Calculation Of SL:"
PRINT #filenum%, "The sextupole strength S is defined by the vertical"
PRINT #filenum%, "field on the x-axis (median plane)."
PRINT #filenum%, "It is the quadratic term in the Taylor expansion."
PRINT #filenum%, " By(x) = 1/2 S x^2"
PRINT #filenum%, "On the x-axis, By(x) = Bth(r=x,th=0)"
PRINT #filenum%, "From the expression for Bth above,"
PRINT #filenum%, " By(x) = - B3 (x/rref)^2 sin(-3 thspole)"
PRINT #filenum%, "Take rref = Rcoil, thspole = pi/6, then at x=Rcoil,"
PRINT #filenum%, " By(Rcoil) = B3 = 1/2 S Rcoil^2"
PRINT #filenum%, "So,"
PRINT #filenum%, " S = 2 B3 / Rcoil^2"
PRINT #filenum%, "The integrated sextupole strength is"
PRINT #filenum%, " SL = 2 BL3 / Rcoil^2"

'Magnetic center calculations
PRINT #filenum%,
PRINT #filenum%, "Calculation Of The Sextupole Center:"
PRINT #filenum%, "In the sextupole's frame,"
PRINT #filenum%, " Bx' = S * x' * y'"
PRINT #filenum%, " By' = 1/2 * S * (x'^2 - y'^2)."
PRINT #filenum%, "In the coil's frame (unprimed frame) the magnetic center is at (x0, y0)."
PRINT #filenum%, "In the coil's frame,"
PRINT #filenum%, " Bx = S * (x - x0) * (y - y0)"
PRINT #filenum%, " By = 1/2 * S * [(x - x0)^2 - (y - y0)^2]."
PRINT #filenum%, "The magnetic center can be found in terms of the measured"
PRINT #filenum%, "quadrupole field.  Compare Br and Btheta on the x-axis to"
PRINT #filenum%, " Bxquad = -S * (x0 * y + y0 * x)"
PRINT #filenum%, " Byquad = -S * (x0 * x - y0 * y)."
PRINT #filenum%, "Evaluate at x = Rcoil, y = 0 and compare to"
PRINT #filenum%, " Br(theta = 0) = Bx, Btheta(theta = 0) = By."
PRINT #filenum%, "In terms of the measured integrated strengths, this gives"
PRINT #filenum%, " Xcenter = x0 = - (1/(SL * Rcoil)) * BL2 * sin(2 * THspole2)"
PRINT #filenum%, " Ycenter = y0 = - (1/(SL * Rcoil)) * BL2 * cos(2 * THspole2)"

'Close the log file
CLOSE filenum%

END SUB

SUB spollogsl (logfile$, bl3!, sbl3!, sl!, ssl!)

'****************************************************************************
'This subroutine write a summary of an average SL measurement to the log file.
'
'Input:
' logfile$, the name of the log file
' bl3!, the integrated field strength of the third harmonic at Rcoil
' sbl3!, the rms variation of bl3
' sl!, the average of the integrated sextupole strength measurements
' ssl!, the rms variation of sl
'
'Zachary Wolf
'8/15/94, 1/3/95, 5/7/95
'****************************************************************************

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

'Print the results
PRINT #filenum%,
PRINT #filenum%, TIME$; "  Integrated Sextupole Strength Calculation:"
PRINT #filenum%, "Coil: Radius = "; coilradiusmP!; " m,"; " # Turns = "; coilnturnsP%
PRINT #filenum%, "Average, # Rotations/Measurement = "; coilnrotaveP%; ", # Measurements = "; blnnmeasaveP%
PRINT #filenum%, "Integrated Field Strength At Rcoil, BL3 = "; bl3!; " +- "; sbl3!; " Tm"
PRINT #filenum%, "Integrated Sextupole Strength, SL = "; sl!; " +- "; ssl!; " T/m"

'Close the log file
CLOSE filenum%

END SUB

SUB spollogxy (logfile$, thn!(), x!, sx!, y!, sy!)

'****************************************************************************
'This subroutine writes the (x, y) location of the magnetic center
'relative to the coil to the log file.
'
'Input:
' logfile$, file to print results
' thn!(1 to nh%), gives the orientation of the sextupole field
' x!, x position of the magnetic center wrt the coil axis
' sx!, error on x
' y!, y position of the magnetic center wrt the coil axis
' sy!, error on y
'
'Zachary Wolf
'1/4/95
'****************************************************************************

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

'Print the magnetic center to the log file
PRINT #filenum%,
PRINT #filenum%,
PRINT #filenum%, "The magnetic center is at"
PRINT #filenum%, "X = ";
PRINT #filenum%, USING "#####.#"; x! * 1000000!;
PRINT #filenum%, " +- ";
PRINT #filenum%, USING "#####.#"; sx! * 1000000!;
PRINT #filenum%, " microns"
PRINT #filenum%, "Y = ";
PRINT #filenum%, USING "#####.#"; y! * 1000000!;
PRINT #filenum%, " +- ";
PRINT #filenum%, USING "#####.#"; sy! * 1000000!;
PRINT #filenum%, " microns"
PRINT #filenum%,
PRINT #filenum%, "The magnetic center is measured relative to the coil axis"
PRINT #filenum%, "effectively at the longitudinal center of the magnet."
PRINT #filenum%, "Pitch and yaw are not seen by the coil."

'Print the following diagram
PRINT #filenum%,
PRINT #filenum%, "View from the non-lead end of the magnet:"
PRINT #filenum%,
IF thn!(3) > 0 THEN
  PRINT #filenum%, "           y          "
  PRINT #filenum%, "           ^          "
  PRINT #filenum%, "           |          "
  PRINT #filenum%, "           N          "
  PRINT #filenum%, "        S     S       "
  PRINT #filenum%, "                --> x "
  PRINT #filenum%, "        N     N       "
  PRINT #filenum%, "           S          "
ELSE
  PRINT #filenum%, "           y          "
  PRINT #filenum%, "           ^          "
  PRINT #filenum%, "           |          "
  PRINT #filenum%, "           S          "
  PRINT #filenum%, "        N     N       "
  PRINT #filenum%, "                --> x "
  PRINT #filenum%, "        S     S       "
  PRINT #filenum%, "           N          "
END IF

'Close the log file
CLOSE filenum%

END SUB

SUB spolmain

'****************************************************************************
'This subroutine controls the magnetic measurements of a sextupole at
'various currents.
'
'Zachary Wolf
'1/4/95, 5/6/95
'****************************************************************************

'Message
PRINT
PRINT "Beginning the sextupole field measurement cycle..."

'Simplify the notation
ni% = nimagtestP%
nh% = spolnhardisplayP%

'Initialize data arrays
DIM bln!(1 TO nh%)                'field strength at Rcoil of the n'th harmonic
DIM sbln!(1 TO nh%)               'standard deviation of Bn
DIM thn!(1 TO nh%)                'field angle
DIM sthn!(1 TO nh%)               'standard deviation of THn
DIM rn!(1 TO nh%)                 'BLn/BL3
DIM srn!(1 TO nh%)                'standard deviation of rn

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

'Ramp to the desired current
CALL imagramp(imagtestP!(i%))

'Measure the current
CALL imaggetiav(imag!, simag!)

'Get the field strength at the coil radius from each harmonic
CALL blngetblnav(nh%, bln!(), sbln!(), thn!(), sthn!())

'Compute the integrated sextupole strength
CALL spolcalcsl(coilradiusmP!, bln!(3), sbln!(3), sl!, ssl!)

'Compute the transfer function
tf! = sl! / imag!
stf! = tf! * SQR((ssl! / sl!) ^ 2 + (simag! / imag!) ^ 2)

'Write the strength measurement to the data file
CALL spoldatsl(sldatP$, imag!, simag!, sl!, ssl!, tf!, stf!)

'Write the strength measurement to the plot file
CALL spolpltsl(slpltP$, imag!, sl!, ssl!, tf!, stf!)

'Do the following only if the harmonics are desired
IF imagharflagP!(i%) > 0 THEN

 'Compute the strength ratios of the harmonics
 CALL harcalcratio(nh%, bln!(), sbln!(), nmain%, rn!(), srn!())

 'Compute the location of the magnetic center
 CALL spolcalcxy(coilradiusmP!, nh%, bln!(), sbln!(), thn!(), sthn!(), x!, sx!, y!, sy!)

 'Write the harmonics to the data file
 CALL hardathar(hardatP$, imag!, simag!, nmain%, nh%, bln!(), sbln!(), thn!(), sthn!(), rn!(), srn!())

 'Write the magnetic center to the data file
 CALL spollogxy(hardatP$, thn!(), x!, sx!, y!, sy!)

 'Write the harmonics to the plot file
 legend$ = "I = " + STR$(imagtestP!(i%)) + " A"
 CALL harplthar(harpltP$, legend$, nh%, rn!(), srn!())

'End harmonics check
END IF

'End loop over currents
NEXT i%

'Write a summary of the calculations to the SL data file
CALL filenewpage(sldatP$)
CALL blnlogcalc(sldatP$)
CALL harlogcalc(sldatP$)
CALL spollogcalc(sldatP$)

'Write a summary of the calculations to the harmonics data file
CALL filenewpage(hardatP$)
CALL blnlogcalc(hardatP$)
CALL harlogcalc(hardatP$)
CALL spollogcalc(hardatP$)

END SUB

SUB spolpltsl (pltfile$, imag!, sl!, ssl!, tf!, stf!)

'****************************************************************************
'This subroutine write a summary of the sl vs I measurements for plotting.
'
'Input:
' pltfile$, the name of the plot file
' imag!, the measured current
' sl!, the integrated sextupole strength at the measured current
' ssl!, the standard deviation of sl
' tf!, the transfer function at the measured current
' stf!, the standard deviation of tf
'
'Zachary Wolf
'7/16/94, 5/6/95
'****************************************************************************

'Keep track of Ncall for the header
STATIC ncall%
ncall% = ncall% + 1

'Open the plot file
filenum% = FREEFILE
OPEN pltfile$ FOR APPEND AS filenum%

'Write the title for plotting on the first call
IF ncall% = 1 THEN
 PRINT #filenum%, "/et g '" + magnameP$ + ", Run " + runP$ + "'"
 PRINT #filenum%, ";i sl ssl tf stf"
END IF

'Write the values
PRINT #filenum%, USING "####.#"; imag!;
PRINT #filenum%, USING " ##.#####"; sl!;
PRINT #filenum%, USING " ##.#####"; ssl!;
PRINT #filenum%, USING " ###.#####"; tf!;
PRINT #filenum%, USING " ###.#####"; stf!

'Close the log file
CLOSE filenum%

END SUB