/*	**************************************************************	*/

/*
 *	Module BLHAR
 *	This module contains functions to determine the integrated field
 *	strength harmonics in a magnet.  The VTcoil system must be
 *	initialized prior to using this module. 
 *
 *	Zachary Wolf
 *	10/12/98
 */

/*	**************************************************************	*/

/*	INCLUDES  */
#include <analysis.h>
#include <utility.h>
#include <formatio.h>
#include <ansi_c.h>
#include "vtcoil.h"
#include "blhar.h"

/*	**************************************************************	*/

/*	PRIVATE PARAMETERS  */
char log_file[80];
char str_dat_file[80];
char str_plt_file[80];
char har_dat_file[80];
char har_plt_file[80];
char mag_ctr_file[80];
static struct blhar_param_struct blhar_param;

/*	**************************************************************	*/

/*	PRIVATE GLOBAL VARIABLES  */
static char msg[80];

/*	**************************************************************	*/

/*	PRIVATE FUNCTIONS  */
int blhar_check_param(void);
void blhar_error(char* message);
void blhar_calc_magnitude_angle(int num_har, double blhar_re_ave[], double blhar_re_rms[], double blhar_im_ave[], double blhar_im_rms[],
		double blhar_mag_ave[], double blhar_mag_rms[], double blhar_spang_ave[], double blhar_spang_rms[]);
void blhar_calc_ratios(int num_har, double blhar_mag_ave[], double blhar_mag_rms[], double blhar_ratio_ave[], double blhar_ratio_rms[]);
void blhar_log_measurement(char* meas);
void blhar_log_calc_blhar(int Nm, double Rm, double Rref, int num_har, double f[], double vthar_re_ave[], double vthar_im_ave[], double blhar_re_ave[], double blhar_im_ave[]);
void blhar_log_blhar_ave(struct blhar_data_struct blhar_data);
void blhar_log_str_ave(int num_str_har, double coil_const, double vt_re_ave, double vt_re_rms, double vt_im_ave, double vt_im_rms, double sl_ave, double sl_rms, double th_ave, double th_rms);

/*	**************************************************************	*/

/*	PUBLIC FUNCTIONS  */

/*	**************************************************************	*/
/*
 *	blhar_init
 *	This function is used to initialize the integrated field strength 
 *	harmonics measurement system.
 *
 *	Zachary Wolf
 *	10/12/98
 */

void blhar_init(char log_file_in[], char str_dat_file_in[], char har_dat_file_in[], char str_plt_file_in[], char har_plt_file_in[], char mag_ctr_file_in[], struct blhar_param_struct blhar_param_in)
{
	/* Declare variables */
	int err;
	
	/* Save parameters for future use */
	strcpy(log_file, log_file_in);
	strcpy(str_dat_file, str_dat_file_in);
	strcpy(har_dat_file, har_dat_file_in);
	strcpy(str_plt_file, str_plt_file_in);
	strcpy(har_plt_file, har_plt_file_in);
	strcpy(mag_ctr_file, mag_ctr_file_in);
	blhar_param = blhar_param_in;
	
	/* Check all parameter values to find any problems */
	err = blhar_check_param();
	if (err != 0)
	{
		blhar_error("Problem with parameter values");
		return;
	}
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_get_blhar_ave
 *	This function measures the field harmonics in a magnet several times.
 *	It computes the average and rms variation of the measurements.
 *
 *	Output:
 *	blhar_data, structure containing BL harmonics measurement results
 *
 *	Zachary Wolf
 *	10/12/98
 */

void blhar_get_blhar_ave(struct blhar_data_struct* blhar_data)
{
	/* Declare variables */
	int err;
	int i, j;
	int num_har;
	int num_print;
	double vthar_re_ave[BLHAR_MAX_NUM_HAR + 1];
	double vthar_re_rms[BLHAR_MAX_NUM_HAR + 1];
	double vthar_im_ave[BLHAR_MAX_NUM_HAR + 1];
	double vthar_im_rms[BLHAR_MAX_NUM_HAR + 1];
	double blhar_re_ave[BLHAR_MAX_NUM_HAR + 1];
	double blhar_re_rms[BLHAR_MAX_NUM_HAR + 1];
	double blhar_im_ave[BLHAR_MAX_NUM_HAR + 1];
	double blhar_im_rms[BLHAR_MAX_NUM_HAR + 1];
	double blhar_mag_ave[BLHAR_MAX_NUM_HAR + 1];
	double blhar_mag_rms[BLHAR_MAX_NUM_HAR + 1];
	double blhar_spang_ave[BLHAR_MAX_NUM_HAR + 1];
	double blhar_spang_rms[BLHAR_MAX_NUM_HAR + 1];
	double blhar_ratio_ave[BLHAR_MAX_NUM_HAR + 1];
	double blhar_ratio_rms[BLHAR_MAX_NUM_HAR + 1];
	int Nm;
	double Rm;
	double Rref;
	double f[BLHAR_MAX_NUM_HAR + 1];
	
	/* Check all module lever parameter values */
	err = blhar_check_param();
	if (err != 0)
	{
		blhar_error("blhar_get_blhar_ave: Problem with parameter values");
		return;
	}
	
	/* Initialize */
	num_har = blhar_param.num_har;
	
	/*
	 * First measure the harmonics with the main coil.
	 * This will give the bucked harmonics if a bucking coil was requested.
	 */
	
	/* Message */
	printf("\nUsing the main coil to measure harmonics...\n");
	 
	/* Log the measurement */
	blhar_log_measurement("Field harmonics with the main coil:");
	
	/* Start with the vtcoil system in a safe state */
	vtcoil_protect();
	
	/* Route the appropriate signal */
	vtcoil_route_signal(blhar_param.coil_param.main_coil_mux_chan, blhar_param.integ_chan);

	/* Set the integrator gain based on the signal level */
	vtcoil_auto_gain();
	
	/* Measure the integrated voltage harmonics */
	vtcoil_get_vthar_ave(blhar_param.num_meas_ave, blhar_param.num_rev_per_meas,
		blhar_param.num_har, vthar_re_ave, vthar_re_rms, vthar_im_ave, vthar_im_rms);

	/* Calculate the integrated field strengths */
	Nm = blhar_param.coil_param.main_coil_num_turns;
	Rm = blhar_param.coil_param.main_coil_radius;			/* (meters) */
	Rref = blhar_param.reference_radius;		/* (meters) */
	for (i = 1; i <= num_har; i++) f[i] = blhar_param.coil_param.main_coil_sens_factors[i];
	for (i = 1; i <= num_har; i++)
	{
		if (f[i] == 0.)
		{
			blhar_re_ave[i] = 0.;
			blhar_re_rms[i] = 0.;
			blhar_im_ave[i] = 0.;
			blhar_im_rms[i] = 0.;
		}
		else
		{
			blhar_re_ave[i] = i * vthar_re_ave[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_re_rms[i] = i * vthar_re_rms[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_im_ave[i] = i * vthar_im_ave[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_im_rms[i] = i * vthar_im_rms[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
		}
	}
	
	/* Log the parameters and results */
	blhar_log_calc_blhar(Nm, Rm, Rref, num_har, f, vthar_re_ave, vthar_im_ave, blhar_re_ave, blhar_im_ave);
	
	/*
	 * If a bucking coil was requested, remeasure the non-bucked harmonics
	 */
	 
	/* Do this part only if a bucking coil was requested */
	if (strcmp(blhar_param.coil_param.coil_for_harmonics, "main") == 0) goto mainonly;
	
	/* Message */
	printf("\nUsing the bucking coil to measure harmonics...\n");
	 
	/* Log the measurement */
	blhar_log_measurement("Field harmonics with the bucking coil:");
	
	/* Put the vtcoil system in a safe state */
	vtcoil_protect();
	
	/* Route the appropriate signal */
	vtcoil_route_signal(blhar_param.coil_param.bucked_coil_mux_chan, blhar_param.integ_chan);

	/* Set the integrator gain based on the signal level */
	vtcoil_auto_gain();
	
	/* Measure the integrated voltage harmonics */
	vtcoil_get_vthar_ave(blhar_param.num_meas_ave, blhar_param.num_rev_per_meas,
		blhar_param.num_har, vthar_re_ave, vthar_re_rms, vthar_im_ave, vthar_im_rms);

	/* Calculate the integrated field strengths for the non-bucked terms */
	Nm = blhar_param.coil_param.main_coil_num_turns;
	Rm = blhar_param.coil_param.main_coil_radius;			/* (meters) */
	Rref = blhar_param.reference_radius;		/* (meters) */
	for (i = 1; i <= num_har; i++) f[i] = blhar_param.coil_param.bucked_coil_sens_factors[i];
	for (i = 1; i <= num_har; i++)
	{
		if (f[i] != 0.)
		{
			blhar_re_ave[i] = i * vthar_re_ave[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_re_rms[i] = i * vthar_re_rms[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_im_ave[i] = i * vthar_im_ave[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
			blhar_im_rms[i] = i * vthar_im_rms[i] / (Nm * Rm * pow(Rm/Rref, i-1) * f[i]);
		}
	}
	
	/* Log the parameters and results */
	blhar_log_calc_blhar(Nm, Rm, Rref, num_har, f, vthar_re_ave, vthar_im_ave, blhar_re_ave, blhar_im_ave);
	
	/* End of bucked coil measurements */
	mainonly:
	
	/* Calculate the magnitudes and south pole angles for each harmonic */
	blhar_calc_magnitude_angle(num_har, blhar_re_ave, blhar_re_rms, blhar_im_ave, blhar_im_rms,
		blhar_mag_ave, blhar_mag_rms, blhar_spang_ave, blhar_spang_rms);
		
	/* Calculate the strength ratios of each harmonic relative to the main harmonic */
	blhar_calc_ratios(num_har, blhar_mag_ave, blhar_mag_rms, blhar_ratio_ave, blhar_ratio_rms);
	
	/* Fill the data structure */
	blhar_data->reference_radius = blhar_param.reference_radius;
	blhar_data->main_coil_radius = blhar_param.coil_param.main_coil_radius;
	blhar_data->num_main_harmonic = blhar_param.num_main_harmonic;
	blhar_data->num_har = num_har;
	for (i = 1; i <= num_har; i++)
	{
		blhar_data->BLn_norm_ave[i] = -blhar_re_ave[i];
		blhar_data->BLn_norm_rms[i] = blhar_re_rms[i];
		blhar_data->BLn_skew_ave[i] = -blhar_im_ave[i];
		blhar_data->BLn_skew_rms[i] = blhar_im_rms[i];
		blhar_data->BLn_mag_ave[i] = blhar_mag_ave[i];
		blhar_data->BLn_mag_rms[i] = blhar_mag_rms[i];
		blhar_data->BLn_spang_ave[i] = blhar_spang_ave[i];
		blhar_data->BLn_spang_rms[i] = blhar_spang_rms[i];
		blhar_data->BLn_ratio_ave[i] = blhar_ratio_ave[i];
		blhar_data->BLn_ratio_rms[i] = blhar_ratio_rms[i];
	}
	
	/* Log the results */
	blhar_log_blhar_ave(*blhar_data);
	
	/* Write the results to the screen */
	printf("\n                 FIELD HARMONICS                 \n");
	printf(" n       BLn          sigBLn       BLn/BLN sBLn/BLN\n");
	printf("         (Tm)          (Tm)          (%%)      (%%)\n");
	printf("---  -------------+-------------  --------+--------\n");
	if (num_har < 10) num_print = num_har;
	else num_print = 10;
	for (i = 1; i <= num_print; i++)
	{
		printf("%3i  %13.6e %13.6e  %8.4f %8.4f\n", i, blhar_mag_ave[i], blhar_mag_rms[i], blhar_ratio_ave[i] * 100., blhar_ratio_rms[i] * 100.);
	}
	
	/* Leave the vtcoil system in a safe state */
	vtcoil_protect();
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_calc_magnetic_center
 *	This function calculates the position of the magnetic center
 *	from the field harmonics measurements.
 *
 *	Input:
 *	blhar_data, structure containing the BL harmonics
 *
 *	Output:
 *	x_ave, x position of the magnetic center (m)
 *	x_rms, rms variation of x position measurements (m)
 *	y_ave, y position of the magnetic center (m)
 *	y_rms, rms variation of y position measurements (m)
 *
 *	Zachary Wolf
 *	8/19/99
 */

void blhar_calc_magnetic_center(struct blhar_data_struct blhar_data, double* x_ave, double* x_rms, double* y_ave, double* y_rms)
{
	/* Declare variables */
	const double pi = 3.1415927;
	double reference_radius = blhar_data.reference_radius;
	int num_main_harmonic = blhar_data.num_main_harmonic;
	double x_rms_sq, y_rms_sq;
	
	/* Dipole */
	if (num_main_harmonic == 1)
	{
		*x_ave = 0.;
		*x_rms = 0.;
		*y_ave = 0.;
		*y_rms = 0.;
	}
	
	/* Higher multipoles */
	else
	{
		*x_ave = - reference_radius * (1. / (num_main_harmonic - 1)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 2)) * 
			(blhar_data.BLn_norm_ave[num_main_harmonic - 1] * blhar_data.BLn_norm_ave[num_main_harmonic] + blhar_data.BLn_skew_ave[num_main_harmonic - 1] * blhar_data.BLn_skew_ave[num_main_harmonic]);
		*y_ave = - reference_radius * (1. / (num_main_harmonic - 1)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 2)) * 
			(blhar_data.BLn_skew_ave[num_main_harmonic - 1] * blhar_data.BLn_norm_ave[num_main_harmonic] - blhar_data.BLn_norm_ave[num_main_harmonic - 1] * blhar_data.BLn_skew_ave[num_main_harmonic]);
		x_rms_sq = pow(reference_radius, 2) * (1. / pow(num_main_harmonic - 1, 2)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 4)) *
			(pow(blhar_data.BLn_skew_ave[num_main_harmonic - 1], 2) * pow(blhar_data.BLn_skew_rms[num_main_harmonic], 2) +
			 pow(blhar_data.BLn_skew_ave[num_main_harmonic], 2) * pow(blhar_data.BLn_skew_rms[num_main_harmonic - 1], 2) +
			 pow(blhar_data.BLn_norm_ave[num_main_harmonic - 1], 2) * pow(blhar_data.BLn_norm_rms[num_main_harmonic], 2) +
			 pow(blhar_data.BLn_norm_ave[num_main_harmonic], 2) * pow(blhar_data.BLn_norm_rms[num_main_harmonic - 1], 2) )
			 +
			 4. * pow(reference_radius, 2) * (1. / pow(num_main_harmonic - 1, 2)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 6)) *
			 pow((blhar_data.BLn_skew_ave[num_main_harmonic - 1] * blhar_data.BLn_skew_ave[num_main_harmonic] + blhar_data.BLn_norm_ave[num_main_harmonic - 1] * blhar_data.BLn_norm_ave[num_main_harmonic]), 2) *
			 pow(blhar_data.BLn_mag_rms[num_main_harmonic], 2);
		y_rms_sq = pow(reference_radius, 2) * (1. / pow(num_main_harmonic - 1, 2)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 4)) *
			(pow(blhar_data.BLn_norm_ave[num_main_harmonic - 1], 2) * pow(blhar_data.BLn_skew_rms[num_main_harmonic], 2) +
			 pow(blhar_data.BLn_skew_ave[num_main_harmonic], 2) * pow(blhar_data.BLn_norm_rms[num_main_harmonic - 1], 2) +
			 pow(blhar_data.BLn_skew_ave[num_main_harmonic - 1], 2) * pow(blhar_data.BLn_norm_rms[num_main_harmonic], 2) +
			 pow(blhar_data.BLn_norm_ave[num_main_harmonic], 2) * pow(blhar_data.BLn_skew_rms[num_main_harmonic - 1], 2) )
			 +
			 4. * pow(reference_radius, 2) * (1. / pow(num_main_harmonic - 1, 2)) * (1. / pow(blhar_data.BLn_mag_ave[num_main_harmonic], 6)) *
			 pow((blhar_data.BLn_norm_ave[num_main_harmonic - 1] * blhar_data.BLn_skew_ave[num_main_harmonic] - blhar_data.BLn_skew_ave[num_main_harmonic - 1] * blhar_data.BLn_norm_ave[num_main_harmonic]), 2) *
			 pow(blhar_data.BLn_mag_rms[num_main_harmonic], 2);
		*x_rms = sqrt(x_rms_sq);
		*y_rms = sqrt(y_rms_sq);
	}
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_dat_blhar_ave
 *	This function writes the results of a BL harmonics measurement
 *	to the data file.
 *
 *	Input:
 *	magnet_current_ave, magnet current (A)
 *	magnet_current_rms, rms variation of Imag measurements (A)
 *	blhar_data, structure containing the BL harmonics
 *	x_ave, x positions of the magnetic center (m)
 *	x_rms, rms variation of the x position measurements (m)
 *	y_ave, y positions of the magnetic center (m)
 *	y_rms, rms variation of the y position measurements (m)
 *
 *	Zachary Wolf
 *	11/28/98
 */

void blhar_dat_blhar_ave(double magnet_current_ave, double magnet_current_rms, struct blhar_data_struct blhar_data, double x_ave, double x_rms, double y_ave, double y_rms)
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	static int num_call = 0;
	
	/* Open the dat file */
	file_ptr = fopen(har_dat_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_dat_blhar_ave: Unable to open harmonics dat file\n");
		return;
	}
	
	/* Write a header on the front page */
	num_call++;
	if (num_call == 1)
	{
		fprintf(file_ptr, "\n\n\n");
		fprintf(file_ptr, "                          Field Harmonics Data\n");
		fprintf(file_ptr, "                        (on the following pages)\n");
	}
	
	/* Write a form feed character to the file so each measurement starts on a new page */
	fprintf(file_ptr, "\f");
	
	/* Write the measured values to the dat file */
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "                          Field Harmonics\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Magnet Current, Imag = %f +- %f A\n", magnet_current_ave, magnet_current_rms);
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Harmonics given at R = %f m\n", blhar_data.reference_radius);
	fprintf(file_ptr, "Coil Radius, Rcoil = %f m\n", blhar_data.main_coil_radius);
	fprintf(file_ptr, "Main Harmonic N = %i\n", blhar_data.num_main_harmonic);
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, " n        BLn         sigBLn     THspole   sigTH    BLn/BLN   sBLn/BLN\n");
	fprintf(file_ptr, "         (Tm)          (Tm)       (deg)    (deg)      (%%)       (%%)   \n");
	fprintf(file_ptr, "---  ------------+------------  --------+--------  ---------+---------\n");
	for (i = 1; i <= blhar_data.num_har; i++)
	{
		fprintf(file_ptr, "%3i  %12.4e %12.4e  %8.2f %8.2f  %9.5f %9.5f\n",
			i,
			blhar_data.BLn_mag_ave[i], blhar_data.BLn_mag_rms[i],
			blhar_data.BLn_spang_ave[i], blhar_data.BLn_spang_rms[i],
			blhar_data.BLn_ratio_ave[i] * 100., blhar_data.BLn_ratio_rms[i] * 100.);
	}
	
	/* Print the location of the magnetic center */
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Magnetic Center:\n");
	fprintf(file_ptr, "Xcenter = %13.8f +- %13.8f m\n", x_ave, x_rms);
	fprintf(file_ptr, "Ycenter = %13.8f +- %13.8f m\n", y_ave, y_rms);
	
	/* Make a picture of the magnet */

	/* Describe where the picture is viewed from */
	fprintf(file_ptr, "\n%s\n", blhar_param.magnet_view);

	/* Dipoles */

	/* Skew dipole at 180 deg */
	if (blhar_data.num_main_harmonic == 1 && (blhar_data.BLn_spang_ave[1] > 135. || blhar_data.BLn_spang_ave[1] < -135.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        S     N       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S     N --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S     N       \n");
	}

	/* Normal dipole at 90 deg */
	else if (blhar_data.num_main_harmonic == 1 && (blhar_data.BLn_spang_ave[1] <= 135. && blhar_data.BLn_spang_ave[1] >= 45.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        S  S  S       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N  N  N       \n");
	}

	/* Skew dipole at 0 deg */
	else if (blhar_data.num_main_harmonic == 1 && (blhar_data.BLn_spang_ave[1] < 45. && blhar_data.BLn_spang_ave[1] > -45.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        N     S       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N     S --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N     S       \n");
	}

	/*Normal dipole at -90 deg */
	else if (blhar_data.num_main_harmonic == 1 && (blhar_data.BLn_spang_ave[1] <= -45. && blhar_data.BLn_spang_ave[1] >= -135.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        N  N  N       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S  S  S       \n");
	}

	/* Quadrupoles */

	/* Skew quadrupole at 90 deg */
	else if (blhar_data.num_main_harmonic == 2 && (blhar_data.BLn_spang_ave[2] > 60. || blhar_data.BLn_spang_ave[2] < -60.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           S          \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N     N --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           S          \n");
	}

	/* Normal quadrupole at 45 deg */
	else if (blhar_data.num_main_harmonic == 2 && (blhar_data.BLn_spang_ave[2] <= 60. && blhar_data.BLn_spang_ave[2] >= 20.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        N     S       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S     N       \n");
	}

	/* Skew quadrupole at 0 deg */
	else if (blhar_data.num_main_harmonic == 2 && (blhar_data.BLn_spang_ave[2] < 20. && blhar_data.BLn_spang_ave[2] > -20.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           N          \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S     S --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           N          \n");
	}

	/* Normal quadrupole at -45 deg */
	else if (blhar_data.num_main_harmonic == 2 && (blhar_data.BLn_spang_ave[2] >= -60. && blhar_data.BLn_spang_ave[2] <= -20.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "        S     N       \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N     S       \n");
	}

	/* Sextupoles */

	/* Skew sextupole at 60 deg */
	else if (blhar_data.num_main_harmonic == 3 && (blhar_data.BLn_spang_ave[3] > 50. || blhar_data.BLn_spang_ave[3] < -50.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "         N   S        \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        S     N --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "         N   S        \n");
	}

	/* Normal sextupole at 30 deg */
	else if (blhar_data.num_main_harmonic == 3 && (blhar_data.BLn_spang_ave[3] <= 50. && blhar_data.BLn_spang_ave[3] >= 10.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           N          \n");
		fprintf(file_ptr, "        S     S       \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "        N     N       \n");
		fprintf(file_ptr, "           S          \n");
	}

	/* Skew sextupole at 0 deg */
	else if (blhar_data.num_main_harmonic == 3 && (blhar_data.BLn_spang_ave[3] < 10. && blhar_data.BLn_spang_ave[3] > -10.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "         S   N        \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "        N     S --> x \n");
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "         S   N        \n");
	}

	/* Normal sextupole at -60 deg */
	else if (blhar_data.num_main_harmonic == 3 && (blhar_data.BLn_spang_ave[3] <= -10. && blhar_data.BLn_spang_ave[3] >= -50.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           S          \n");
		fprintf(file_ptr, "        N     N       \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "        S     S       \n");
		fprintf(file_ptr, "           N          \n");
	}

	/* Octupoles */

	/* Skew octupole at 45 deg */
	else if (blhar_data.num_main_harmonic == 4 && (blhar_data.BLn_spang_ave[4] > 35. || blhar_data.BLn_spang_ave[4] < -35.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           N          \n");
		fprintf(file_ptr, "        S     S       \n");
		fprintf(file_ptr, "       N       N--> x \n");
		fprintf(file_ptr, "        S     S       \n");
		fprintf(file_ptr, "           N          \n");
	}

	/* Normal octupole at 45/2 deg */
	else if (blhar_data.num_main_harmonic == 4 && (blhar_data.BLn_spang_ave[4] <= 35. && blhar_data.BLn_spang_ave[4] >= 10.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "         S   N        \n");
		fprintf(file_ptr, "       N       S      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "       S       N      \n");
		fprintf(file_ptr, "         N    S       \n");
	}

	/* Skew octupole at 0 deg */
	else if (blhar_data.num_main_harmonic == 4 && (blhar_data.BLn_spang_ave[4] < 10. && blhar_data.BLn_spang_ave[4] > -10.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "           S          \n");
		fprintf(file_ptr, "        N     N       \n");
		fprintf(file_ptr, "       S       S--> x \n");
		fprintf(file_ptr, "        N     N       \n");
		fprintf(file_ptr, "           S          \n");
	}

	/* Normal octupole at -45/2 deg */
	else if (blhar_data.num_main_harmonic == 4 && (blhar_data.BLn_spang_ave[4] <= -10. && blhar_data.BLn_spang_ave[4] >= -35.))
	{
		fprintf(file_ptr, "                      \n");
		fprintf(file_ptr, "           y          \n");
		fprintf(file_ptr, "           ^          \n");
		fprintf(file_ptr, "           |          \n");
		fprintf(file_ptr, "         N   S        \n");
		fprintf(file_ptr, "       S       N      \n");
		fprintf(file_ptr, "                --> x \n");
		fprintf(file_ptr, "       N       S      \n");
		fprintf(file_ptr, "         S    N       \n");
	}

	/* Unknown pole configuration */
	else
	{
		/* Don't print anything */
	}
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_dat_mag_ctr
 *	This function writes the results of a magnetic center determination
 *	to the magnetic center data file.
 *
 *	Input:
 *	i_mag, magnet current (A)
 *	x_ave, x positions of the magnetic center (m)
 *	x_rms, rms variation of the x position measurements (m)
 *	y_ave, y positions of the magnetic center (m)
 *	y_rms, rms variation of the y position measurements (m)
 *
 *	Zachary Wolf
 *	11/10/00
 */

void blhar_dat_mag_ctr(double i_mag, double x_ave, double x_rms, double y_ave, double y_rms)
{
	/* Declare variables */
	FILE* file_ptr;
	static int num_call = 0;
	
	/* Open the dat file */
	file_ptr = fopen(mag_ctr_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_dat_mag_ctr: Unable to open magnetic center dat file\n");
		return;
	}
	
	/* Write a header on the first call */
	num_call++;
	if (num_call == 1)
	{
		fprintf(file_ptr, "\n\n\n");
		fprintf(file_ptr, "                          Magnetic Center Data\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "  Time      Imag        X       sigX         Y       sigY   \n");
		fprintf(file_ptr, "            (A)          (microns)            (microns)     \n");
		fprintf(file_ptr, "--------  --------  ---------+---------  ---------+---------\n");
	}
	
	/* Print the location of the magnetic center */
	fprintf(file_ptr, "%s  %8.3f  %9.2f %9.2f  %9.2f %9.2f\n", TimeStr(), i_mag, x_ave * 1.e6, x_rms * 1.e6, y_ave * 1.e6, y_rms * 1.e6);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_plt_blhar_ave
 *	This function writes the results of a BL harmonics measurement
 *	to the plotting file.
 *
 *	Input:
 *	legend, legend for the new harmonics being added to the plot file 
 *	blhar_data, structure containing the BL harmonics
 *
 *	Zachary Wolf
 *	12/3/98
 */

void blhar_plt_blhar_ave(char* legend, struct blhar_data_struct blhar_data)
{
	/* Local parameters */
	#define MAX_NUM_CHAR_IN 1000

	/* Declare variables */
	FILE* file_ptr;
	FILE* temp_ptr;
	int i;
	static int num_call = 0;
	char har_1 = '0';
	char buffer[MAX_NUM_CHAR_IN];
	char* buf_ptr;
	int buf_len;
	int err;
	
	/* Increment num_call */
	num_call++;
	
	/* Open the plt file */
	file_ptr = fopen(har_plt_file, "r");
	if (file_ptr == NULL)
	{
		printf("blhar_plt_blhar_ave: Unable to open harmonics plt file\n");
		return;
	}
	
	/* Open the temp file */
	temp_ptr = fopen("hartemp.dat", "w");
	if (temp_ptr == NULL)
	{
		printf("blhar_plt_blhar_ave: Unable to open temp file\n");
		return;
	}
	
	/* On the first call, write data directly to the plt file */
	if (num_call == 1)
	{
		fprintf(temp_ptr, ";n r sr r sr... for different currents\n");		/* Column contents */
		fprintf(temp_ptr, "/sa l '%s' 2\n", legend);						/* Legend for first data columns */
		for (i = 1; i <= blhar_data.num_har; i++)							/* Harmonics in columns */
		{
			fprintf(temp_ptr, "%3i %9.6f %9.6f\n", i, blhar_data.BLn_ratio_ave[i], blhar_data.BLn_ratio_rms[i]);
		}
	}
	
	/* After the first call, read the data in the plt file, append the new data, */
	/* and save the results in the temporary file*/
	if (num_call > 1)
	{
		while (har_1 != '1')		/* Copy the plot file to the temp file until the harmonics are reached */
		{
			buf_ptr = fgets(buffer, MAX_NUM_CHAR_IN, file_ptr);
			if (buf_ptr == NULL)
			{
				printf("blhar_plt_blhar_ave: Problem reading harmonics plot file.\n");
				return;
			}
			har_1 = buffer[2];
			if (har_1 != '1') fputs(buffer, temp_ptr);
		}
		fprintf(temp_ptr, "/sa l '%s' %i\n", legend, 2 * num_call);			/* Legend for data columns */
		buf_len = strlen(buffer);
		buffer[buf_len - 1] = '\0';											/* Remove newline character */
		fprintf(temp_ptr, "%s %9.6f %9.6f\n", buffer, blhar_data.BLn_ratio_ave[1], blhar_data.BLn_ratio_rms[1]);	/* Harmonic 1 */
		for (i = 2; i <= blhar_data.num_har; i++)							/* Harmonics 2 to N in columns */
		{
			buf_ptr = fgets(buffer, MAX_NUM_CHAR_IN, file_ptr);
			if (buf_ptr == NULL)
			{
				printf("blhar_plt_blhar_ave: Problem reading harmonics plot file.\n");
				return;
			}
			buf_len = strlen(buffer);
			buffer[buf_len - 1] = '\0';											/* Remove newline character */
			fprintf(temp_ptr, "%s %9.6f %9.6f\n", buffer, blhar_data.BLn_ratio_ave[i], blhar_data.BLn_ratio_rms[i]);
		}
	}
	
	/* Close the plt file */
	fclose(file_ptr);
	
	/* Close the temp file */
	fclose(temp_ptr);
	
	/* Make the temp file the plt file */
	err = CopyFile("hartemp.dat", har_plt_file);
	if (err != 0) blhar_error("blhar_plt_blhar_ave: Problem copying file");
	
	/* Delete the temp file */
	err = DeleteFile("hartemp.dat");
	if (err != 0) blhar_error("blhar_plt_blhar_ave: Problem deleting file");
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_get_str_ave
 *	This function measures the strength of a magnet.
 *	It returns the integrated field strength of a dipole, the
 *	integrated gradient of a quadrupole, ...
 *
 *	Output:
 *	num_str_harm, number of the field harmonic the strength is being returned for
 *	sl_ave, integrated strength (Tm^2/m^n), average
 *	sl_rms, integrated strength (Tm^2/m^n), rms variation
 *	th_ave, angle of the first south pole (deg), average
 *	th_rms, angle of the first south pole (deg), rms variation
 *
 *	Zachary Wolf
 *	7/8/99
 */

void blhar_get_str_ave(int* num_str_harm, double* sl_ave, double* sl_rms, double* th_ave, double* th_rms)
{
	/* Declare variables */
	int err;
	int num_har;
	int num_str_har;
	double vthar_re_ave[BLHAR_MAX_NUM_HAR + 1];
	double vthar_re_rms[BLHAR_MAX_NUM_HAR + 1];
	double vthar_im_ave[BLHAR_MAX_NUM_HAR + 1];
	double vthar_im_rms[BLHAR_MAX_NUM_HAR + 1];
	double vt_re_ave, vt_re_rms;
	double vt_im_ave, vt_im_rms;
	double vt_mag_ave, vt_mag_rms;
	
	/* Check all module lever parameter values */
	err = blhar_check_param();
	if (err != 0)
	{
		blhar_error("blhar_get_dipole_str_ave: Problem with parameter values");
		return;
	}
	
	/* The type of strength coil is given in the coil parameters */
	if (strcmp(blhar_param.coil_param.str_coil_type, "dipole") == 0) num_str_har = 1;
	else if (strcmp(blhar_param.coil_param.str_coil_type, "quad") == 0) num_str_har = 2;
	else if (strcmp(blhar_param.coil_param.str_coil_type, "sext") == 0) num_str_har = 3;
	else if (strcmp(blhar_param.coil_param.str_coil_type, "oct") == 0) num_str_har = 4;
	else
	{
		blhar_error("blhar_get_str_ave: unknown strength coil type");
		return;
	}
	
	/* Return the number of the strength harmonic */
	*num_str_harm = num_str_har;
	
	/* Determine how many harmonics to measure with the strength coil */
	num_har = num_str_har + 2;
	
	/* Message */
	printf("\nMeasuring the magnet strength...\n");
	 
	/* Log the measurement */
	blhar_log_measurement("Magnet strength measurement:");
	
	/* Start with the vtcoil system in a safe state */
	vtcoil_protect();
	
	/* Route the appropriate signal */
	vtcoil_route_signal(blhar_param.coil_param.str_coil_mux_chan, blhar_param.integ_chan);

	/* Set the integrator gain based on the signal level */
	vtcoil_auto_gain();
	
	/* Measure the integrated voltage harmonics */
	vtcoil_get_vthar_ave(blhar_param.num_meas_ave, blhar_param.num_rev_per_meas,
		num_har, vthar_re_ave, vthar_re_rms, vthar_im_ave, vthar_im_rms);

	/* Select the appropriate integrated voltage harmonic */
	vt_re_ave = vthar_re_ave[num_str_har];
	vt_re_rms = vthar_re_rms[num_str_har];
	vt_im_ave = vthar_im_ave[num_str_har];
	vt_im_rms = vthar_im_rms[num_str_har];
	
	/* Calculate the magnitude of the integrated voltage */
	vt_mag_ave = sqrt( pow(vt_re_ave, 2) + pow(vt_im_ave, 2) );
	if (vt_mag_ave == 0.)
		vt_mag_rms = 0.;
	else
		vt_mag_rms = (1. / vt_mag_ave) * sqrt( pow(vt_re_ave, 2) * pow(vt_re_rms, 2) + pow(vt_im_ave, 2) * pow(vt_im_rms, 2) );

	/* Calculate the magnet strength */
	*sl_ave = vt_mag_ave * blhar_param.coil_param.str_coil_constant;
	*sl_rms = vt_mag_rms * blhar_param.coil_param.str_coil_constant;
	
	/* Calculate the south pole angle */
	*th_ave = - (1. / num_str_har) * (atan2(vt_im_ave, vt_re_ave) + 3.1415927 / 2.);
	if (*th_ave > 3.1415927 / num_str_har) *th_ave = *th_ave - 2. * 3.1415927 / num_str_har;
	if (*th_ave < -3.1415927 / num_str_har) *th_ave = *th_ave + 2. * 3.1415927 / num_str_har;
	if (vt_mag_ave == 0.)
		*th_rms = 0.;
	else
		*th_rms = (1. / num_str_har) * (1. / pow(vt_mag_ave, 2)) * sqrt( pow(vt_im_ave, 2) * pow(vt_re_rms, 2) + pow(vt_re_ave, 2) * pow(vt_im_rms, 2) );
	
	/* Convert the south pole angle to degrees */
	*th_ave = *th_ave * 180. / 3.1415927;
	*th_rms = *th_rms * 180. / 3.1415927;
	
	/* Log the results */
	blhar_log_str_ave(num_str_har, blhar_param.coil_param.str_coil_constant, vt_re_ave, vt_re_rms, vt_im_ave, vt_im_rms, *sl_ave, *sl_rms, *th_ave, *th_rms);
	
	/* Write the results to the screen */
	if (num_str_har == 1) printf("\nIntegrated Field Strength BL = %f +- %f Tm\n", *sl_ave, *sl_rms);
	if (num_str_har == 2) printf("\nIntegrated Gradient GL = %f +- %f T\n", *sl_ave, *sl_rms);
	if (num_str_har == 3) printf("\nIntegrated Sextupole Strength SL = %f +- %f T/m\n", *sl_ave, *sl_rms);
	if (num_str_har == 4) printf("\nIntegrated Octupole Strength OL = %f +- %f T/m^2\n", *sl_ave, *sl_rms);
	printf("South Pole Angle TH = %f +- %f deg\n", *th_ave, *th_rms);
	
	/* Leave the vtcoil system in a safe state */
	vtcoil_protect();
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_dat_str_ave
 *	This function writes the results of a magnet strength measurement
 *	to the data file.
 *
 *	Input:
 *	num_str_har, number of the harmonic the strength is measured for
 *	magnet_current_ave, magnet current (A)
 *	magnet_current_rms, rms variation of Imag measurements (A)
 *	sl_ave, integrated strength (Tm^2/m^n), average
 *	sl_rms, integrated strength (Tm^2/m^n), rms variation
 *	th_ave, angle of the first south pole (deg), average
 *	th_rms, angle of the first south pole (deg), rms variation
 *
 *	Zachary Wolf
 *	6/29/99
 */

void blhar_dat_str_ave(int num_str_har, double magnet_current_ave, double magnet_current_rms, double sl_ave, double sl_rms, double th_ave, double th_rms)
{
	/* Declare variables */
	FILE* file_ptr;
	static int call_num;
	double tf_ave, tf_rms;
	
	/* Check input parameters */
	if (num_str_har < 1 || num_str_har > 4)
	{
		printf("blhar_dat_str_ave: unknown num_str_har");
		return;
	}
	
	/* Open the dat file */
	file_ptr = fopen(str_dat_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_dat_str_ave: Unable to open magnet strength dat file\n");
		return;
	}
	
	/* Increment the number of calls */
	call_num++;
	
	/* Write a header on the first call */
	if (call_num == 1 && num_str_har == 1)
	{
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "                Integrated Field Strength\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "  Imag    sigImag      BL       sigBL        BL/I     sigBL/I   \n");
		fprintf(file_ptr, "  (A)       (A)       (Tm)       (Tm)       (Tm/A)    (Tm/A)    \n");
		fprintf(file_ptr, "--------+--------  ----------+----------  ----------+---------- \n");
	}
	if (call_num == 1 && num_str_har == 2)
	{
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "                   Integrated Gradient\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "  Imag    sigImag      GL       sigGL        GL/I     sigGL/I   \n");
		fprintf(file_ptr, "  (A)       (A)        (T)       (T)         (T/A)     (T/A)    \n");
		fprintf(file_ptr, "--------+--------  ----------+----------  ----------+---------- \n");
	}
	if (call_num == 1 && num_str_har == 3)
	{
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "             Integrated Sextupole Strength\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "  Imag    sigImag      SL       sigSL        SL/I     sigSL/I   \n");
		fprintf(file_ptr, "  (A)       (A)       (T/m)     (T/m)       (T/m/A)   (T/m/A)   \n");
		fprintf(file_ptr, "--------+--------  ----------+----------  ----------+---------- \n");
	}
	if (call_num == 1 && num_str_har == 4)
	{
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "             Integrated Octupole Strength\n");
		fprintf(file_ptr, "\n");
		fprintf(file_ptr, "  Imag    sigImag      OL       sigOL        OL/I     sigOL/I   \n");
		fprintf(file_ptr, "  (A)       (A)      (T/m^2)   (T/m^2)          (T/m^2/A)       \n");
		fprintf(file_ptr, "--------+--------  ----------+----------  ----------+---------- \n");
	}
	
	/* Compute the transfer function */
	if (fabs(magnet_current_ave) > .01)
	{
		tf_ave = sl_ave / magnet_current_ave;
		tf_rms = sl_rms / magnet_current_ave;
	}
	else
	{
		tf_ave = 0.;
		tf_rms = 0.;
	}
	
	/* Write the data */
	fprintf(file_ptr, "%8.3f %8.3f  %10.5f %10.5f  %10.7f %10.7f\n",
		magnet_current_ave, magnet_current_rms,
		sl_ave, sl_rms, tf_ave, tf_rms);
	
	/* Close the dat file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_plt_str_ave
 *	This function writes the results of a magnet strength measurement
 *	to the plot file.
 *
 *	Input:
 *	num_str_har, number of the harmonic the strength is measured for
 *	magnet_current_ave, magnet current (A)
 *	magnet_current_rms, rms variation of Imag measurements (A)
 *	sl_ave, integrated strength (Tm^2/m^n), average
 *	sl_rms, integrated strength (Tm^2/m^n), rms variation
 *	th_ave, angle of the first south pole (deg), average
 *	th_rms, angle of the first south pole (deg), rms variation
 *
 *	Zachary Wolf
 *	6/29/99
 */

void blhar_plt_str_ave(int num_str_har, double magnet_current_ave, double magnet_current_rms, double sl_ave, double sl_rms, double th_ave, double th_rms)
{
	/* Declare variables */
	FILE* file_ptr;
	static int call_num;
	double tf_ave, tf_rms;
	
	/* Check input parameters */
	if (num_str_har < 1 || num_str_har > 4)
	{
		printf("blhar_plt_str_ave: unknown num_str_har");
		return;
	}
	
	/* Open the plt file */
	file_ptr = fopen(str_plt_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_plt_str_ave: Unable to open magnet strength plt file\n");
		return;
	}
	
	/* Increment the number of calls */
	call_num++;
	
	/* Write a header on the first call */
	if (call_num == 1)
	{
		fprintf(file_ptr, ";Imag, sImag, SL, sSL, TF, sTF\n");
	}
	
	/* Compute the transfer function */
	if (fabs(magnet_current_ave) > .01)
	{
		tf_ave = sl_ave / magnet_current_ave;
		tf_rms = sl_rms / magnet_current_ave;
	}
	else
	{
		tf_ave = 0.;
		tf_rms = 0.;
	}
	
	/* Write the data */
	fprintf(file_ptr, "%8.3f %8.3f  %8.5f %8.5f  %10.7f %10.7f\n",
		magnet_current_ave, magnet_current_rms,
		sl_ave, sl_rms, tf_ave, tf_rms);
	
	/* Close the dat file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/

/*	INTERNAL FUNCTIONS  */

/*	**************************************************************	*/
/*
 *	blhar_error
 *	This function handles errors for the IMAG module.
 *
 *	Input:
 *	message, string to display in standard I/O
 *
 *	Zachary Wolf
 *	10/11/98
 */
 
void blhar_error(char* message)
{
	/* Declare variables */
	char buf[80];
	
	/* Notify the operator of the error */
	printf("\nERROR: %s\n", message);
	Beep();
	Delay(.5);
	Beep();
	
	/* Terminate the program if the operator desires */
	printf("Press ENTER to continue.\nPress any key then ENTER to terminate program.\n");
	fgets(buf, 80, stdin);
	if (buf[0] == '\n') return;
	else exit(0);
}

/*	**************************************************************	*/
/*
 *	blhar_check_param
 *	This function checks that all required parameters have been assigned
 *	appropriate values.
 *
 *	Output:
 *	err, 0 = no error, non-zero = error
 *
 *	Zachary Wolf
 *	8/11/98
 */

int blhar_check_param()
{
	/* Declare variables */
	int i;
	
	/* log_file */
	if (CompareStrings(log_file, 0, "", 0, 0) == 0)
	{
		blhar_error("blhar_check_param: Log file not defined.");
		return -1;
	}
	
	/* str_dat_file */
	if (CompareStrings(str_dat_file, 0, "", 0, 0) == 0)
	{
		blhar_error("blhar_check_param: Magnet strength dat file not defined.");
		return -1;
	}
	
	/* har_dat_file */
	if (CompareStrings(har_dat_file, 0, "", 0, 0) == 0)
	{
		blhar_error("blhar_check_param: Harmonics dat file not defined.");
		return -1;
	}

	/* str_plt_file */
	if (CompareStrings(str_plt_file, 0, "", 0, 0) == 0)
	{
		blhar_error("blhar_check_param: Magnet strength plt file not defined.");
		return -1;
	}

	/* har_plt_file */
	if (CompareStrings(har_plt_file, 0, "", 0, 0) == 0)
	{
		blhar_error("blhar_check_param: Harmonics plt file not defined.");
		return -1;
	}

	/* num_har */
	if (blhar_param.num_har < 1 || blhar_param.num_har > 40)
	{
		blhar_error("blhar_check_param: num_har is invalid");
		return -1;
	}

	/* num_main_harmonic */
	if (blhar_param.num_main_harmonic < 1 || blhar_param.num_main_harmonic > 12)
	{
		blhar_error("blhar_check_param: num_main_harmonic is invalid");
		return -1;
	}

	/* num_meas_ave */
	if (blhar_param.num_meas_ave < 1 || blhar_param.num_meas_ave > 10)
	{
		blhar_error("blhar_check_param: num_meas_ave is invalid");
		return -1;
	}

	/* integ_chan */
	if (blhar_param.integ_chan != 'A' && blhar_param.integ_chan != 'B')
	{
		blhar_error("blhar_check_param: integ_chan is invalid");
		return -1;
	}

	/* num_rev_per_meas */
	if (blhar_param.num_rev_per_meas < 1 || blhar_param.num_rev_per_meas > 32)
	{
		blhar_error("blhar_check_param: num_rev_per_meas is invalid");
		return -1;
	}

	/* reference_radius */
	if (blhar_param.reference_radius <= 0. || blhar_param.reference_radius > 5.)
	{
		blhar_error("blhar_check_param: reference_radius is invalid");
		return -1;
	}
	
	/* str_coil_type */
	if (strcmp(blhar_param.coil_param.str_coil_type, "dipole") != 0 &&
		strcmp(blhar_param.coil_param.str_coil_type, "quad") != 0 &&
		strcmp(blhar_param.coil_param.str_coil_type, "sext") != 0 &&
		strcmp(blhar_param.coil_param.str_coil_type, "oct") != 0)
	{
		blhar_error("blhar_check_param: str_coil_type is invalid");
		return -1;
	}

	/* str_coil_mux_chan */
	if (blhar_param.coil_param.str_coil_mux_chan < 1 || blhar_param.coil_param.str_coil_mux_chan > 40)
	{
		blhar_error("blhar_check_param: str_coil_mux_chan is invalid");
		return -1;
	}

	/* str_coil_constant */
	if (blhar_param.coil_param.str_coil_constant < 0. || blhar_param.coil_param.str_coil_constant > 1.e6)
	{
		blhar_error("blhar_check_param: str_coil_constant is invalid");
		return -1;
	}

	/* main_coil_radius */
	if (blhar_param.coil_param.main_coil_radius <= 0. || blhar_param.coil_param.main_coil_radius > 5.)
	{
		blhar_error("blhar_check_param: main_coil_radius is invalid");
		return -1;
	}

	/* main_coil_num_turns */
	if (blhar_param.coil_param.main_coil_num_turns < 1 || blhar_param.coil_param.main_coil_num_turns > 500)
	{
		blhar_error("blhar_check_param: main_coil_num_turns is invalid");
		return -1;
	}

	/* main_coil_mux_chan */
	if (blhar_param.coil_param.main_coil_mux_chan < 1 || blhar_param.coil_param.main_coil_mux_chan > 40)
	{
		blhar_error("blhar_check_param: main_coil_mux_chan is invalid");
		return -1;
	}

	/* bucked_coil_mux_chan */
	if (blhar_param.coil_param.bucked_coil_mux_chan < 1 || blhar_param.coil_param.bucked_coil_mux_chan > 40)
	{
		blhar_error("blhar_check_param: bucked_coil_mux_chan is invalid");
		return -1;
	}
	
	/* coil_for_harmonics */
	if (strcmp(blhar_param.coil_param.coil_for_harmonics, "main") != 0 &&
		strcmp(blhar_param.coil_param.coil_for_harmonics, "bucked") != 0)
	{
		blhar_error("blhar_check_param: coil_for_harmonics is invalid");
		return -1;
	}

	/* num_sens_factors */
	if (blhar_param.coil_param.num_sens_factors < blhar_param.num_har || blhar_param.coil_param.num_sens_factors > BLHAR_MAX_NUM_HAR)
	{
		blhar_error("blhar_check_param: num_sens_factors is invalid");
		return -1;
	}
	
	/* main_coil_sens_factors */
	for (i = 1; i <= blhar_param.coil_param.num_sens_factors; i++)
	{
		if (fabs(blhar_param.coil_param.main_coil_sens_factors[i]) > 10.)
		{
			blhar_error("blhar_check_param: main_coil_sens_factor is invalid");
			return -1;
		}
	}
	
	/* bucked_coil_sens_factors */
	for (i = 1; i <= blhar_param.coil_param.num_sens_factors; i++)
	{
		if (fabs(blhar_param.coil_param.bucked_coil_sens_factors[i]) > 10.)
		{
			blhar_error("blhar_check_param: bucked_coil_sens_factor is invalid");
			return -1;
		}
	}
	
	/* If we made it this far, all parameters have values */
	return 0;
}

/*	**************************************************************	*/
/*
 *	blhar_calc_magnitude_angle
 *	This function calculates the field harmonics in a magnet.
 *	It uses the real and imaginary field expansion coefficients to
 *	calculate the harmonic field amplitude and the angle of the first
 *	south pole.
 *
 *	Input:
 *	num_har, number of harmonics to record
 *	blhar_re_ave[0 to num_har], bl harmonics (Tm), real part, [0] is unused
 *	blhar_re_rms[0 to num_har], rms variation of bl harmonics (Tm), real part, [0] is unused
 *	blhar_im_ave[0 to num_har], bl harmonics (Tm), imaginary part, [0] is unused
 *	blhar_im_rms[0 to num_har], rms variation of bl harmonics (Tm), imaginary part, [0] is unused
 *
 *	Output:
 *	blhar_mag_ave[0 to num_har], bl harmonics magnitudes (Tm), [0] is unused
 *	blhar_mag_rms[0 to num_har], rms variation of bl harmonics magnitudes (Tm), [0] is unused
 *	blhar_spang_ave[0 to num_har], south pole angles (deg), [0] is unused
 *	blhar_spang_rms[0 to num_har], rms variation of south pole angles (deg), [0] is unused
 *
 *	Zachary Wolf
 *	11/11/98
 */

void blhar_calc_magnitude_angle(int num_har, double blhar_re_ave[], double blhar_re_rms[], double blhar_im_ave[], double blhar_im_rms[],
		double blhar_mag_ave[], double blhar_mag_rms[], double blhar_spang_ave[], double blhar_spang_rms[])
{
	/* Declare variables */
	int i;
	
	/* Check input parameters */
	if (num_har > BLHAR_MAX_NUM_HAR)
	{
		blhar_error("blhar_calc_magnitude_angle: problem with num_har, exiting");
		return;
	}
	
	/* Calculate the field harmonic magnitudes */
	for (i = 1; i <= num_har; i++)
	{
		blhar_mag_ave[i] = sqrt(pow(blhar_re_ave[i], 2) + pow(blhar_im_ave[i], 2));
		if (blhar_mag_ave[i] == 0.)
			blhar_mag_rms[i] = 0.;
		else
			blhar_mag_rms[i] = (1. / blhar_mag_ave[i]) * sqrt(pow(blhar_re_ave[i],2) * pow(blhar_re_rms[i], 2) + pow(blhar_im_ave[i], 2) * pow(blhar_im_rms[i], 2));
	}
	
	/* Calculate the field harmonic south pole angles, convert to degrees */
	for (i = 1; i <= num_har; i++)
	{
		blhar_spang_ave[i] = - (1. / i) * (atan2(blhar_im_ave[i], blhar_re_ave[i]) + 3.1415927 / 2.);
		if (blhar_spang_ave[i] > 3.1415927 / i) blhar_spang_ave[i] = blhar_spang_ave[i] - 2. * 3.1415927 / i;
		if (blhar_spang_ave[i] < -3.1415927 / i) blhar_spang_ave[i] = blhar_spang_ave[i] + 2. * 3.1415927 / i;
		if (blhar_mag_ave[i] == 0.)
			blhar_spang_rms[i] = 0.;
		else
			blhar_spang_rms[i] = (1. / i) * (1. / pow(blhar_mag_ave[i], 2)) * sqrt( pow(blhar_im_ave[i], 2) * pow(blhar_re_rms[i], 2) + pow(blhar_re_ave[i], 2) * pow(blhar_im_rms[i], 2) );
		blhar_spang_ave[i] = blhar_spang_ave[i] * 180. / 3.1415927;
		blhar_spang_rms[i] = blhar_spang_rms[i] * 180. / 3.1415927;
	}
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_calc_ratios
 *	This function calculates the ratios of the field harmonic strengths
 *	to the strength of the main harmonic.
 *
 *	Input:
 *	num_har, number of harmonics to record
 *	blhar_mag_ave[0 to num_har], bl harmonics magnitudes (Tm), [0] is unused
 *	blhar_mag_rms[0 to num_har], rms variation of bl harmonics magnitudes (Tm), [0] is unused
 *
 *	Output:
 *	blhar_ratio_ave[0 to num_har], BL(n) / BL(num_main_harmonic), [0] is unused
 *	blhar_ratio_rms[0 to num_har], rms variation of BL(n) / BL(num_main_harmonic), [0] is unused
 *
 *	Zachary Wolf
 *	11/16/98
 */

void blhar_calc_ratios(int num_har, double blhar_mag_ave[], double blhar_mag_rms[], double blhar_ratio_ave[], double blhar_ratio_rms[])
{
	/* Declare variables */
	int i;
	int num_main_harmonic = blhar_param.num_main_harmonic;
	
	/* Check input parameters */
	if (num_har > BLHAR_MAX_NUM_HAR)
	{
		blhar_error("blhar_calc_ratios: problem with num_har, exiting");
		return;
	}
	
	/* Calculate the field harmonic ratios to the main harmonic */
	for (i = 1; i <= num_har; i++)
	{
		blhar_ratio_ave[i] = blhar_mag_ave[i] / blhar_mag_ave[num_main_harmonic];
		blhar_ratio_rms[i] = blhar_ratio_ave[i] * sqrt( pow((blhar_mag_rms[i] / (blhar_mag_ave[i] + .00000001)), 2) + pow((blhar_mag_rms[num_main_harmonic] / (blhar_mag_ave[num_main_harmonic] + .00000001)), 2) );
	}
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_log_measurement
 *	This function logs a measurement to the log file.
 *
 *	Input:
 *	meas, a description of the measurement for the log file
 *
 *	Zachary Wolf
 *	7/8/99
 */
 
void blhar_log_measurement(char* meas)
{
	/* Declare variables */
	FILE* file_ptr;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_log_measurement: Unable to open log file\n");
		return;
	}
	
	/* Write to the log file */
	fprintf(file_ptr, "%s  Magnetic Field Measurement\n", TimeStr());
	fprintf(file_ptr, "          %s\n", meas);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_log_calc_blhar
 *	This function logs the result of a BL harmonics calculation.
 *
 *	Input:
 *	Nm, number of turns on the main winding
 *	Rm, radius of the outer turns on the main winding (m)
 *	Rref, reference radius (m)
 *	num_har, number of harmonics to record
 *	f[0 to num_har], sensitivity factors, real, [0] is unused
 *	vthar_re_ave[0 to num_har], vt harmonics (Vs), real part, [0] is unused
 *	vthar_im_ave[0 to num_har], vt harmonics (Vs), imaginary part, [0] is unused
 *	blhar_re_ave[0 to num_har], BL harmonics (Tm), real part, [0] is unused
 *	blhar_im_ave[0 to num_har], BL harmonics (Tm), imaginary part, [0] is unused
 *
 *	Zachary Wolf
 *	7/8/99
 */

void blhar_log_calc_blhar(int Nm, double Rm, double Rref, int num_har, double f[], double vthar_re_ave[], double vthar_im_ave[], double blhar_re_ave[], double blhar_im_ave[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_log_calc_blhar: Unable to open log file\n");
		return;
	}
	
	/* Write the results to the log file */
	fprintf(file_ptr, "%s  Field Harmonics Calculation\n", TimeStr());
	fprintf(file_ptr, "          Main Coil Radius, Rm = %f m\n", Rm);
	fprintf(file_ptr, "          Main Coil Number Of Turns, Nm = %i\n", Nm);
	fprintf(file_ptr, "          Reference Radius, Rref = %f m\n", Rref);
	fprintf(file_ptr, "           n   vthar_re[i]   vthar_im[i]     Bucking     blhar_re[i]   blhar_im[i] \n");
	fprintf(file_ptr, "                   (Vs)          (Vs)        Factor          (Tm)          (Tm)    \n");
	fprintf(file_ptr, "          ---  ------------  ------------  ------------  ------------  ------------\n");
	for (i = 1; i <= num_har; i++)
	{
		fprintf(file_ptr, "          %3i  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f\n", i, vthar_re_ave[i], vthar_im_ave[i], f[i], blhar_re_ave[i], blhar_im_ave[i]);
	}
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_log_blhar_ave
 *	This function logs the result of a BL harmonics measurement.
 *	It writes the average and rms variation of BL harmonics measurements
 *	to the log file.
 *
 *	Input:
 *	blhar_data, results of a BL harmonics measurement
 *
 *	Zachary Wolf
 *	11/10/98
 */

void blhar_log_blhar_ave(struct blhar_data_struct blhar_data)
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_log_blhar_ave: Unable to open log file\n");
		return;
	}
	
	/* Write the measured values to the log file */
	fprintf(file_ptr, "%s  Field Harmonics, Average and RMS of Measurements\n", TimeStr());
	fprintf(file_ptr, "           n   blhar_nm[i]   sblhar_nm[i]  blhar_sk[i]   sblhar_sk[i]  blhar_mag[i]  sblhar_mg[i]     THsp[i]      sTHsp[i]       BLn/BLN      sBLn/BLN  \n");
	fprintf(file_ptr, "                   (Tm)          (Tm)          (Tm)          (Tm)          (Tm)          (Tm)          (deg)         (deg)          (%%)           (%%)    \n");
	fprintf(file_ptr, "          ---  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------  ------------\n");
	for (i = 1; i <= blhar_data.num_har; i++)
	{
		fprintf(file_ptr, "          %3i  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f  %12.7f\n",
			i,
			blhar_data.BLn_norm_ave[i], blhar_data.BLn_norm_rms[i],
			blhar_data.BLn_skew_ave[i], blhar_data.BLn_skew_rms[i],
			blhar_data.BLn_mag_ave[i], blhar_data.BLn_mag_rms[i],
			blhar_data.BLn_spang_ave[i], blhar_data.BLn_spang_rms[i],
			blhar_data.BLn_ratio_ave[i] * 100., blhar_data.BLn_ratio_rms[i] * 100.);
	}
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blhar_log__str_ave
 *	This function logs the result of a magnet strength measurement.
 *
 *	Input:
 *	num_str_har, number of the field harmonic the strength is being measured for
 *	coil_const, coil constant
 *	vt_re_ave, integrated voltage (Vs), real part, average
 *	vt_re_rms, integrated voltage (Vs), real part, rms variation
 *	vt_im_ave, integrated voltage (Vs), imaginary part, average
 *	vt_im_rms, integrated voltage (Vs), imaginary part, rms variation
 *	sl_ave, integrated strength (Tm^2/m^n), average
 *	sl_rms, integrated strength (Tm^2/m^n), rms variation
 *	th_ave, angle of the first south pole (deg), average
 *	th_rms, angle of the first south pole (deg), rms variation
 *
 *	Zachary Wolf
 *	6/29/99
 */

void blhar_log_str_ave(int num_str_har, double coil_const, double vt_re_ave, double vt_re_rms, double vt_im_ave, double vt_im_rms, double sl_ave, double sl_rms, double th_ave, double th_rms)
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blhar_log_str_ave: Unable to open log file\n");
		return;
	}
	
	/* Write the measured values to the log file */
	fprintf(file_ptr, "%s  Magnet Strength Measurement\n", TimeStr());
	fprintf(file_ptr, "          Harmonic Number = %i\n", num_str_har);
	fprintf(file_ptr, "          Coil Constant = %12.7f Tm^2/m^n/Vs\n", coil_const);
	fprintf(file_ptr, "          VT_re = %12.7f +- %12.7f Vs\n", vt_re_ave, vt_re_rms);
	fprintf(file_ptr, "          VT_im = %12.7f +- %12.7f Vs\n", vt_im_ave, vt_im_rms);
	fprintf(file_ptr, "          SL    = %12.7f +- %12.7f Tm^2/m^n\n", sl_ave, sl_rms);
	fprintf(file_ptr, "          THsp  = %12.7f +- %12.7f deg\n", th_ave, th_rms);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}