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

/*
 *	Module BLSCAN
 *	This module contains functions to determine the integrated field
 *	strength of a magnet as a function of position.  The VTscan system
 *	must be initialized prior to using this module. 
 *
 *	Zachary Wolf
 *	3/16/00
 */

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

/*	INCLUDES  */
#include <userint.h>
#include <analysis.h>
#include <utility.h>
#include <formatio.h>
#include <ansi_c.h>
#include "vtscan.h"
#include "vtscanui.h"
#include "blscan.h"
#include "blscanui.h"

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

/*	PRIVATE PARAMETERS  */
char log_file[BLSCAN_MAX_NUM_CHAR + 1];
char dat_file[BLSCAN_MAX_NUM_CHAR + 1];
char plt_file[BLSCAN_MAX_NUM_CHAR + 1];
static struct blscan_param_struct blscan_param;

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

/*	PRIVATE GLOBAL VARIABLES  */
static char msg[BLSCAN_MAX_NUM_CHAR + 1];

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

/*	PRIVATE FUNCTIONS  */
void blscan_error(char* message);
int blscan_check_param(void);
double blscan_fact(int n);
void blscan_log_measurement(void);
void blscan_log_get_vt(void);
void blscan_log_calc_vt_coeff(int order, double coeff[]);
void blscan_log_M_double_coil(double M[8][8], double invM[8][8]);
void blscan_log_M_single_coil(double M[9][9], double invM[9][9]);
void blscan_log_calc_bln(int num_har, double bln[]);
void blscan_log_get_bln_ave(int order, double vtn_ave[], double vtn_rms[], int num_har, double bln_ave[], double bln_rms[]);

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

/*	PUBLIC FUNCTIONS  */

/*	**************************************************************	*/
/*
 *	blscan_init
 *	This function is used to initialize the BLwire system. 
 *
 *	Zachary Wolf
 *	11/24/99
 */

void blscan_init(char log_file_in[], char dat_file_in[], char plt_file_in[], struct blscan_param_struct blscan_param_in)
{
	/* Declare variables */
	int err;
	
	/* Save parameters for future use */
	strcpy(log_file, log_file_in);
	strcpy(dat_file, dat_file_in);
	strcpy(plt_file, plt_file_in);
	blscan_param = blscan_param_in;
	
	/* Check all parameter values to find any problems */
	err = blscan_check_param();
	if (err != 0)
	{
		blscan_error("Problem with parameter values");
		return;
	}
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_get_vt
 *	This function measures the integrated voltage from a coil as it
 *	is moved across the magnet.  The coil is moved forward, then backward
 *	to take out the integrator drift.
 *
 *	Input:
 *	x_start, starting position of the wire motion (m)
 *	dx, distance the wire moves (m)
 *	num_samp, number of VT samples to take as the wire is moved
 *
 *	Output:
 *	x[0 to n], x[0] is the starting position, x[1 to n] are the positions of the samples (m)
 *	vt[0 to n], vt[0] = 0, vt[1 to n] = n measured integrated voltage samples (Vs)
 *
 *	Zachary Wolf
 *	11/24/99
 */

void blscan_get_vt(double x_start, double dx, int num_samp, double x[], double vt[])
{
	/* Declare variables */
	int err;
	double vt_for[BLSCAN_MAX_NUM_SAMP + 1];
	double vt_rev[BLSCAN_MAX_NUM_SAMP + 1];
	double cvt_rev[BLSCAN_MAX_NUM_SAMP + 1];
	double x_for[BLSCAN_MAX_NUM_SAMP + 1];
	double x_rev[BLSCAN_MAX_NUM_SAMP + 1];
	double cx_rev[BLSCAN_MAX_NUM_SAMP + 1];
	int i;
	double c0;
	
	/* Check input parameters */
	if (fabs(x_start) > 1.)
	{
		blscan_error("blscan_get_vt: x_start has improper value");
		return;
	}
	if (fabs(dx) < 1.e-6 || fabs(dx) > 1.)
	{
		blscan_error("blscan_get_vt: dx has improper value");
		return;
	}
	if (num_samp < 1 || num_samp > BLSCAN_MAX_NUM_SAMP)
	{
		blscan_error("blscan_get_vt: num_samp has improper value");
		return;
	}
	
	/* Check all module lever parameter values */
	err = blscan_check_param();
	if (err != 0)
	{
		blscan_error("blscan_get_vt: Problem with parameter values");
		return;
	}
	
	/* Perform the forward scan */
	vtscan_get_vt(x_start, dx, num_samp, x_for, vt_for);
	
	/* Clear the graph of the samples */
	vtscanui_clear_vt_vs_x();
	
	/* Plot the samples */
	vtscanui_scale_vt_vs_x(num_samp + 1, x_for);
	vtscanui_plot_vt_vs_x(num_samp, x_for, vt_for, VAL_GREEN);
	
	/* Save the samples */
	blscan_save_vt_samp("vt_for.dat", num_samp, x_for, vt_for);
	
	/* Perform the reverse scan, start at the last sample */
	vtscan_get_vt(x_for[num_samp], -dx, num_samp, x_rev, vt_rev);
	
	/* Plot the samples */
	vtscanui_plot_vt_vs_x(num_samp, x_rev, vt_rev, VAL_BLUE);
	
	/* Save the samples */
	blscan_save_vt_samp("vt_rev.dat", num_samp, x_rev, vt_rev);

	/* Reverse the order of the reverse scan samples */
	for (i = 0; i <= num_samp; i++)
	{
		cx_rev[i] = x_rev[num_samp - i];
		cvt_rev[i] = vt_rev[num_samp - i];
	}
	
	/* Add an offset so the reverse scan samples start at 0 */
	c0 = cvt_rev[0];
	for (i = 0; i <= num_samp; i++)
	{
		cvt_rev[i] = cvt_rev[i] - c0;
	}
	
	/* Plot the samples */
	vtscanui_plot_vt_vs_x(num_samp, cx_rev, cvt_rev, VAL_CYAN);
	
	/* Average the forward scan and the corrected reverse scan to take out drift */
	for (i = 0; i <= num_samp; i++)
	{
		vt[i] = (vt_for[i] + cvt_rev[i]) / 2.;
		x[i] = x_for[i];
	}
	
	/* Plot the samples */
	vtscanui_plot_vt_vs_x(num_samp, x, vt, VAL_RED);
	
	/* Save the samples */
	blscan_save_vt_samp("vt_ave.dat", num_samp, x, vt);
							  
	/* Log the results */
	blscan_log_get_vt();
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_calc_vt_coeff
 *	This function performs a polynomial fit to the data from a
 *	VT vs X measurement.
 *
 *	Input:
 *	num_samp, number of VT samples to take as the wire is moved
 *	x[0 to n], x[0] is the starting position, x[1 to n] are the positions of the samples (m)
 *	vt[0 to n], vt[0] = 0, vt[1 to n] = n measured integrated voltage samples (Vs)
 *	order, order of the polynomial fit
 *
 *	Output:
 *	x_norm[0 to n], x positions divided by the normalization radius R0 (dimensionless)
 *	vt_coeff[0 to order], the VTi in VT(x) = VT0 + VT1*(x/R0) + VT2*(x/R0)^2 + ... + VTorder*(x/R0)^order (Vs)
 *	vt_fit[0 to n], fitted integrated voltage values at the normalized x positions (Vs)
 *
 *	Zachary Wolf
 *	3/17/00
 */

void blscan_calc_vt_coeff(int num_samp, double x[], double vt[], int order, double x_norm[], double vt_coeff[], double vt_fit[])
{
	/* Declare variables */
	int err;
	double mse;
	int i;
	
	/* Check input parameters */
	if (num_samp < 1 || num_samp > BLSCAN_MAX_NUM_SAMP)
	{
		blscan_error("blscan_calc_coeff: num_samp has improper value");
		return;
	}
	if (order < 0 || order > num_samp/2)
	{
		blscan_error("blscan_calc_coeff: order has improper value");
		return;
	}
	
	/* Calculate the normalized x positions */
	for (i = 0; i <= num_samp; i++) x_norm[i] = x[i] / blscan_param.norm_radius;
	
	/* Perform a polynomial fit to the VT(x) measurements */
	err = PolyFit(x_norm, vt, num_samp + 1, order, vt_fit, vt_coeff, &mse);
	if (err != 0)
	{
		blscan_error("blscan_calc_vt_coeff: error in polynomial fit");
		return;
	}
					  
	/* Clear the graphs of the samples */
	blscanui_clear_vt_vs_x();
	blscanui_clear_vt_coeff();
	
	/* Scale the horizontal axes of the graphs */
	blscanui_scale_vt_vs_x(num_samp + 1, x_norm);
	blscanui_scale_vt_coeff(order);
							 
	/* Plot the original and fitted points */
	blscanui_plot_vt_vs_x(num_samp, x_norm, vt, VAL_RED);
	blscanui_plot_vt_vs_x(num_samp, x_norm, vt_fit, VAL_BLUE);
	
	/* Plot the fit coefficients */
	blscanui_scale_vt_coeff(order);
	blscanui_plot_vt_coeff(order, vt_coeff);

	/* Log the calculation */		  
	blscan_log_calc_vt_coeff(order, vt_coeff);
	
	/* Done */
	return;			
}
			
/*	**************************************************************	*/
/*
 *	blscan_calc_bln_double_coil
 *	This function calculates the field expansion coefficients from
 *	the coefficients in the fit to the VT vs X data.  It assumes the
 *	measurement was made with the double coil.  The algorithm requires
 *	that at least an eighth order polynomial was fit to the VT vs X data.
 *
 *	Input:
 *	order, order of the polynomial fit
 *	vt_coeff[0 to order], the VTi in VT(x) = VT0 + VT1*(x/R0) + VT2*(x/R0)^2 + ... + VTorder*(x/R0)^order (Vs)
 *	num_har, number of field harmonics to record
 *
 *	Output:
 *	bln[0 to num_har], the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *
 *	Zachary Wolf
 *	8/18/00
 */

void blscan_calc_bln_double_coil(int order, double vt_coeff[], int num_har, double bln[])
{
	/* Declare variables */
	int err;
	double VT[8][1];
	double BL[8][1];
	double M[8][8];
	double invM[8][8];
	int j, n, ij, in;
	
	/* Check input parameters */
	if (order < 8 || order < num_har - 2)
	{
		blscan_error("blscan_calc_bln_double_coil: order has improper value");
		return;
	}
	if (num_har < 3 || num_har > BLSCAN_MAX_NUM_HAR)
	{
		blscan_error("blscan_calc_bln_double_coil: num_har has improper value");
		return;
	}
	
	/* Form the matrix relating BL to VT */
	for (ij = 0; ij <= 7; ij++)
	{
		for (in = 0; in <= 7; in++)
		{
			j = ij + 1;
			n = in + 3;
			if (in < ij) M[ij][in] = 0.;
			else M[ij][in] = blscan_param.num_turns_coil * blscan_param.norm_radius * (1./n)
							 * (1./pow(blscan_param.norm_radius,n-j))
							 * (blscan_fact(n) / (blscan_fact(j) * blscan_fact(n-j)))
							 * (pow(blscan_param.outer_coil_pos, n-j) - pow(blscan_param.inner_coil_pos, n-j)
							    - pow(-blscan_param.inner_coil_pos, n-j) + pow(-blscan_param.outer_coil_pos, n-j));
		}
	}
	
	/* Prepare the VT polynomial fit coefficient array */
	for (ij = 0; ij < 8; ij++) VT[ij][0] = vt_coeff[ij + 1];
	
	/* Find the inverse of M */
	err = InvMatrix(M, 8, invM);
	if (err != 0)
	{
		blscan_error("blscan_calc_bln_double_coil: error in matrix inversion");
		return;
	}
	
	/* Multiply the inverse of M by VT to get BL */
	err = MatrixMul(invM, VT, 8, 8, 1, BL);
	if (err != 0)
	{
		blscan_error("blscan_calc_bln_double_coil: error in matrix multiplication");
		return;
	}
	
	/* Get the BL coefficients from BL */
	bln[0] = 0.;
	bln[1] = 0.;
	bln[2] = 0.;
	for (n = 3; n <= num_har; n++) bln[n] = BL[n-3][0];

	/* Log the calculation */		  
	blscan_log_M_double_coil(M, invM);
	blscan_log_calc_bln(num_har, bln);
	
	/* Done */
	return;			
}
			
/*	**************************************************************	*/
/*
 *	blscan_calc_bln_single_coil
 *	This function calculates the field expansion coefficients from
 *	the coefficients in the fit to the VT vs X data.  It assumes the
 *	measurement was made with the outer single coil.  The algorithm requires
 *	that at least a ninth order polynomial was fit to the VT vs X data.
 *
 *	Input:
 *	order, order of the polynomial fit
 *	vt_coeff[0 to order], the VTi in VT(x) = VT0 + VT1*(x/R0) + VT2*(x/R0)^2 + ... + VTorder*(x/R0)^order (Vs)
 *	num_har, number of field harmonics to record
 *
 *	Output:
 *	bln[0 to num_har], the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *
 *	Zachary Wolf
 *	8/21/00
 */

void blscan_calc_bln_single_coil(int order, double vt_coeff[], int num_har, double bln[])
{
	/* Declare variables */
	int err;
	double VT[9][1];
	double BL[9][1];
	double M[9][9];
	double invM[9][9];
	int j, n, ij, in;
	
	/* Check input parameters */
	if (order < 9 || order < num_har - 1)
	{
		blscan_error("blscan_calc_bln_single_coil: order has improper value");
		return;
	}
	if (num_har < 3 || num_har > BLSCAN_MAX_NUM_HAR)
	{
		blscan_error("blscan_calc_bln_double_coil: num_har has improper value");
		return;
	}
	
	/* Form the matrix relating BL to VT */
	for (ij = 0; ij <= 8; ij++)
	{
		for (in = 0; in <= 8; in++)
		{
			j = ij + 1;
			n = in + 2;
			if (in < ij) M[ij][in] = 0.;
			else M[ij][in] = blscan_param.num_turns_coil * blscan_param.norm_radius * (1./n)
							 * (1./pow(blscan_param.norm_radius, n-j))
							 * (blscan_fact(n) / (blscan_fact(j) * blscan_fact(n-j)))
							 * (pow((blscan_param.outer_coil_pos - blscan_param.inner_coil_pos) / 2., n-j) - pow(-(blscan_param.outer_coil_pos - blscan_param.inner_coil_pos) / 2., n-j));
		}
	}
	
	/* Prepare the VT polynomial fit coefficient array */
	for (ij = 0; ij < 9; ij++) VT[ij][0] = vt_coeff[ij + 1];
	
	/* Find the inverse of M */
	err = InvMatrix(M, 9, invM);
	if (err != 0)
	{
		blscan_error("blscan_calc_bln_double_coil: error in matrix inversion");
		return;
	}
	
	/* Multiply the inverse of M by VT to get BL */
	err = MatrixMul(invM, VT, 9, 9, 1, BL);
	if (err != 0)
	{
		blscan_error("blscan_calc_bln_double_coil: error in matrix multiplication");
		return;
	}
	
	/* Get the BL coefficients from BL */
	bln[0] = 0.;
	bln[1] = 0.;
	for (n = 2; n <= num_har; n++) bln[n] = BL[n-2][0];

	/* Log the calculation */		  
	blscan_log_M_single_coil(M, invM);
	blscan_log_calc_bln(num_har, bln);
	
	/* Done */
	return;			
}
			
/*	**************************************************************	*/
/*
 *	blscan_calc_bln_wire
 *	This function calculates the field expansion coefficients from
 *	the coefficients in the fit to the VT vs X data.  It assumes the
 *	measurement was made with a wire or wire bundle.  The return
 *	path of the circuit is fixed.
 *
 *	Input:
 *	order, order of the polynomial fit
 *	vt_coeff[0 to order], the VTi in VT(x) = VT0 + VT1*(x/R0) + VT2*(x/R0)^2 + ... + VTorder*(x/R0)^order (Vs)
 *	num_har, number of field harmonics to record
 *
 *	Output:
 *	bln[0 to num_har], the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *
 *	Zachary Wolf
 *	9/19/00
 */

void blscan_calc_bln_wire(int order, double vt_coeff[], int num_har, double bln[])
{
	/* Declare variables */
	int err;
	int n;
	
	/* Check input parameters */
	if (order < 1)
	{
		blscan_error("blscan_calc_bln_wire: order has improper value");
		return;
	}
	if (num_har < 1 || num_har > BLSCAN_MAX_NUM_HAR)
	{
		blscan_error("blscan_calc_bln_wire: num_har has improper value");
		return;
	}
	
	/* Calculate the BL coefficients from the VT coefficients */
	for (n = 1; n <= order; n++) bln[n] = (n / (blscan_param.num_turns_wire * blscan_param.norm_radius)) * vt_coeff[n];
	for (n = order + 1; n <= num_har; n++) bln[n] = 0.;

	/* Log the calculation */		  
	blscan_log_calc_bln(num_har, bln);
	
	/* Done */
	return;			
}

/*	**************************************************************	*/
/*
 *	blscan_get_bln_ave
 *	This function performs several VT vs X scans.  For each scan,
 *	it calculates the BL harmonic strengths.  It then finds the
 *	average and rms variation of the measurements.
 *
 *	Input:
 *	config, wire or coil configuration used to make to measurements
 *	num_har, number of harmonics to report (<= 10)
 *
 *	Output:
 *	bln_ave[0 to num_har], bln_ave[0] = 0, bln_ave[1 to n] are the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *	bln_rms[0 to num_har], bln_rms[0] = 0, bln_rms[1 to n] are the rms variation of BLn measurements (Tm)
 *
 *	Zachary Wolf
 *	3/17/00
 */

void blscan_get_bln_ave(enum blscan_config config, int num_har, double bln_ave[], double bln_rms[])
{
	/* Declare variables */
	int err;
	int i, j;
	int vt_fit_order = 0;
	double x[BLSCAN_MAX_NUM_SAMP + 1];
	double x_norm[BLSCAN_MAX_NUM_SAMP + 1];
	double vt[BLSCAN_MAX_NUM_SAMP + 1];
	double vt_fit[BLSCAN_MAX_NUM_SAMP + 1];
	double vt_coeff[BLSCAN_MAX_NUM_HAR + 1];  
	double bln[BLSCAN_MAX_NUM_HAR + 1];
	double vt_coeff_meas[BLSCAN_MAX_NUM_HAR + 1][BLSCAN_MAX_NUM_MEAS + 1];
	double bln_meas[BLSCAN_MAX_NUM_HAR + 1][BLSCAN_MAX_NUM_MEAS + 1];
	double meas_values[BLSCAN_MAX_NUM_MEAS + 1];
	double vtn_ave[BLSCAN_MAX_NUM_HAR + 1];
	double vtn_rms[BLSCAN_MAX_NUM_HAR + 1];
	
	/* Check input parameters */
	if (num_har < 3 || num_har > BLSCAN_MAX_NUM_HAR)
	{
		blscan_error("blscan_get_bln_ave: improper num_har");
		return;
	}
	
	/* Check module level parameters */
	err = blscan_check_param();
	if (err != 0)
	{
		blscan_error("Problem with parameter values");
		return;
	}
		
	/* When using a single coil, shift the origin so the coil is centered on the scan range */
	if (config == BLSCAN_SINGLE_COIL)
	{
		printf("\nShifting the origin so the single coil is centered...\n");
		vtscan_move_wire_abs(-(blscan_param.outer_coil_pos + blscan_param.inner_coil_pos) / 2.);
		vtscan_zero_pos();
	}
		
	/* When using a single wire or wire bundle, shift the origin so the wire is centered on the scan range */
	if (config == BLSCAN_WIRE)
	{
		printf("\nShifting the origin so the wire is centered...\n");
		vtscan_move_wire_abs(-blscan_param.wire_pos);
		vtscan_zero_pos();
	}
	
	/* Message */
	printf("\nPerforming measurement #");
	
	/* Loop over the measurements, save the results */
	for (j = 1; j <= blscan_param.num_meas_ave; j++)
	{
		printf("%2i", j);												   
		if (config == BLSCAN_DOUBLE_COIL)
		{
			vtscan_close_mux_chan(blscan_param.mux_chan_double_coil);
			vtscan_set_gain(blscan_param.gain_double_coil);		
			Delay(.5);
			blscan_get_vt(blscan_param.scan_start_pos, blscan_param.scan_dist, blscan_param.scan_num_samp, x, vt);
			vtscan_protect();
			vt_fit_order = 8;
			blscan_calc_vt_coeff(blscan_param.scan_num_samp, x, vt, vt_fit_order, x_norm, vt_coeff, vt_fit);
			blscan_calc_bln_double_coil(vt_fit_order, vt_coeff, num_har, bln);
		}
		else if (config == BLSCAN_SINGLE_COIL)
		{
			vtscan_close_mux_chan(blscan_param.mux_chan_single_coil);
			vtscan_set_gain(blscan_param.gain_single_coil);
			Delay(.5);
			blscan_get_vt(blscan_param.scan_start_pos, blscan_param.scan_dist, blscan_param.scan_num_samp, x, vt);
			vtscan_protect();
			vt_fit_order = 9;
			blscan_calc_vt_coeff(blscan_param.scan_num_samp, x, vt, vt_fit_order, x_norm, vt_coeff, vt_fit);
			blscan_calc_bln_single_coil(vt_fit_order, vt_coeff, num_har, bln);
		}
		else if (config == BLSCAN_WIRE)
		{
			vtscan_close_mux_chan(blscan_param.mux_chan_wire);
			vtscan_set_gain(blscan_param.gain_wire);
			Delay(.5);
			blscan_get_vt(blscan_param.scan_start_pos, blscan_param.scan_dist, blscan_param.scan_num_samp, x, vt);
			vtscan_protect();
			vt_fit_order = 10;
			blscan_calc_vt_coeff(blscan_param.scan_num_samp, x, vt, vt_fit_order, x_norm, vt_coeff, vt_fit);
			blscan_calc_bln_wire(vt_fit_order, vt_coeff, num_har, bln);
		}
		else blscan_error("Unknown coil configuration");
		for (i = 1; i <= vt_fit_order; i++) vt_coeff_meas[i][j] = vt_coeff[i];
		for (i = 1; i <= num_har; i++) bln_meas[i][j] = bln[i];
	}
	
	/* End messages */
	printf("\n");
	
	/* When using a single coil, shift the origin back to its original position */
	if (config == BLSCAN_SINGLE_COIL)
	{
		printf("\nShifting the origin back to its original position...\n");
		vtscan_move_wire_abs((blscan_param.outer_coil_pos + blscan_param.inner_coil_pos) / 2.);
		vtscan_zero_pos();
	}
	
	/* When using a wire, shift the origin back to its original position */
	if (config == BLSCAN_WIRE)
	{
		printf("\nShifting the origin back to its original position...\n");
		vtscan_move_wire_abs(blscan_param.wire_pos);
		vtscan_zero_pos();
	}
	
	/* Calculate the average and rms variation of the measurements */
	vtn_ave[0] = 0.;
	vtn_rms[0] = 0.;
	for (i = 1; i <= vt_fit_order; i++)
	{
		for (j = 1; j <= blscan_param.num_meas_ave; j++) meas_values[j-1] = vt_coeff_meas[i][j];
		StdDev(meas_values, blscan_param.num_meas_ave, &vtn_ave[i], &vtn_rms[i]);
	}
	bln_ave[0] = 0.;
	bln_rms[0] = 0.;
	for (i = 1; i <= num_har; i++)
	{
		for (j = 1; j <= blscan_param.num_meas_ave; j++) meas_values[j-1] = bln_meas[i][j];
		StdDev(meas_values, blscan_param.num_meas_ave, &bln_ave[i], &bln_rms[i]);
	}
	
	/* Log the result */
	blscan_log_get_bln_ave(vt_fit_order, vtn_ave, vtn_rms, num_har, bln_ave, bln_rms);
	
	/* Done */
	return;
}
						
/*	**************************************************************	*/
/*
 *	blscan_dat_bln_config_ave
 *	This function writes the results of a magnet strength measurement
 *	to the data file.
 *
 *	Input:
 *	config, coil or wire configuration used to make the measurements
 *	imag_ave, average of magnet current measurements (A)
 *	imag_rms, rms variation of magnet current measurements (A)
 *	num_har, number of harmonics to report (<= 10)
 *	bln_ave[0 to num_har], bln_ave[0] = 0, bln_ave[1 to n] are the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *	bln_rms[0 to num_har], bln_rms[0] = 0, bln_rms[1 to n] are the rms variation of BLn measurements (Tm)
 *
 *	Zachary Wolf
 *	3/22/00
 */

void blscan_dat_bln_config_ave(enum blscan_config config, double imag_ave, double imag_rms, int num_har, double bln_ave[], double bln_rms[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the dat file */
	file_ptr = fopen(dat_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_dat_bln_config_ave: Unable to open data file\n");
		return;
	}
	
	/* Write a header */
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Integrated Field Strength Harmonics\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Magnet Current = %12.4f +- %12.4f A\n", imag_ave, imag_rms);
	fprintf(file_ptr, "\n");
	if (config == BLSCAN_DOUBLE_COIL) fprintf(file_ptr, "Coil Configuration: Double Coil\n");
	if (config == BLSCAN_SINGLE_COIL) fprintf(file_ptr, "Coil Configuration: Single Coil\n");
	if (config == BLSCAN_WIRE)        fprintf(file_ptr, "Coil Configuration: Wire\n");
	fprintf(file_ptr, "BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)\n");
	fprintf(file_ptr, "R0 = %12.6f m\n", blscan_param.norm_radius);
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "            BLn_ave      BLn_rms  \n");
	fprintf(file_ptr, "   n         (Tm)         (Tm)    \n");
	fprintf(file_ptr, "-------  ------------+------------\n");
	for (i = 1; i <= num_har; i++)
		fprintf(file_ptr, "%4i     %12.6f %12.6f\n", i, bln_ave[i], bln_rms[i]);
	
	/* Close the dat file */
	fclose(file_ptr);
	
	/* Done */
	return;
}
						
/*	**************************************************************	*/
/*
 *	blscan_dat_bln_ave
 *	This function writes the results of a magnet strength measurement
 *	to the data file.
 *
 *	Input:
 *	imag_ave, average of magnet current measurements (A)
 *	imag_rms, rms variation of magnet current measurements (A)
 *	num_har, number of harmonics to report (<= 10)
 *	bln_ave[0 to num_har], bln_ave[0] = 0, bln_ave[1 to n] are the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *	bln_rms[0 to num_har], bln_rms[0] = 0, bln_rms[1 to n] are the rms variation of BLn measurements (Tm)
 *
 *	Zachary Wolf
 *	3/22/00
 */

void blscan_dat_bln_ave(double imag_ave, double imag_rms, int num_har, double bln_ave[], double bln_rms[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the dat file */
	file_ptr = fopen(dat_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_dat_bln_ave: Unable to open data file\n");
		return;
	}
	
	/* Write a header */
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Integrated Field Strength Harmonics Summary\n");
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "Magnet Current = %12.4f +- %12.4f A\n", imag_ave, imag_rms);
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)\n");
	fprintf(file_ptr, "R0 = %12.6f m\n", blscan_param.norm_radius);
	fprintf(file_ptr, "\n");
	fprintf(file_ptr, "            BLn_ave      BLn_rms      BLn/BL1      BLn/BL1   \n");
	fprintf(file_ptr, "   n         (Tm)         (Tm)          ave          rms     \n");
	fprintf(file_ptr, "-------  ------------+------------  ------------+------------\n");
	for (i = 1; i <= num_har; i++)
		fprintf(file_ptr, "%4i     %12.6f %12.6f  %12.8f %12.8f\n", i, bln_ave[i], bln_rms[i], bln_ave[i] / bln_ave[1], sqrt(pow(bln_rms[i], 2) / pow(bln_ave[1], 2) + pow(bln_rms[1], 2) * pow(bln_ave[i], 2) / pow(bln_ave[1], 4)));
	
	/* Close the dat file */
	fclose(file_ptr);
	
	/* Done */
	return;
}
						
/*	**************************************************************	*/
/*
 *	blscan_plt_bln_ave
 *	This function writes the results of a magnet strength measurement
 *	to the data file.
 *
 *	Input:
 *	imag_ave, average of magnet current measurements (A)
 *	imag_rms, rms variation of magnet current measurements (A)
 *	num_har, number of harmonics to report (<= 10)
 *	bln_ave[0 to num_har], bln_ave[0] = 0, bln_ave[1 to n] are the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *	bln_rms[0 to num_har], bln_rms[0] = 0, bln_rms[1 to n] are the rms variation of BLn measurements (Tm)
 *
 *	Zachary Wolf
 *	3/22/00
 */

void blscan_plt_bln_ave(double imag_ave, double imag_rms, int num_har, double bln_ave[], double bln_rms[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the dat file */
	file_ptr = fopen(plt_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_plt_bln_ave: Unable to open plot file\n");
		return;
	}
	
	/* Write a header */
	fprintf(file_ptr, "%%n imag_ave imag_rms bln_ave bln_rms\n");
	
	/* Write the data */
	for (i = 1; i <= num_har; i++)
		fprintf(file_ptr, "%5i  %12.5f  %12.5f  %12.6f  %12.6f\n", i, imag_ave, imag_rms, bln_ave[i], bln_rms[i]);
	
	/* Close the dat file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_save_vt_samp
 *	This function writes the VT vs X samples to a file.
 *
 *	Input:
 *	file_name, name of file to store the samples in
 *	num_samp, number of VT samples to take as the wire is moved
 *	x[0 to n], x[0] is the starting position, x[1 to n] are the positions of the samples (m)
 *	vt[0 to n], vt[0] = 0, vt[1 to n] = n measured integrated voltage samples (Vs)
 *
 *	Zachary Wolf
 *	3/23/00
 */

void blscan_save_vt_samp(char file_name[], int num_samp, double x[], double vt[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(file_name, "w");
	if (file_ptr == NULL)
	{
		printf("blscan_save_vt_samp: Unable to open save file\n");
		return;
	}
	
	/* Write the data to the file */
	fprintf(file_ptr, "%% i, x, vt\n");
	for (i = 0; i <= num_samp; i++)
		fprintf(file_ptr, "%5i  %12.5f  %12.5f\n", i, x[i], vt[i]);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

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

/*	INTERNAL FUNCTIONS  */

/*	**************************************************************	*/
/*
 *	blscan_error
 *	This function handles errors for the IMAG module.
 *
 *	Input:
 *	message, string to display in standard I/O
 *
 *	Zachary Wolf
 *	10/11/98
 */
 
void blscan_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);
}

/*	**************************************************************	*/
/*
 *	blscan_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 blscan_check_param()
{
	/* Declare variables */
	int i;
	
	/* log_file */
	if (CompareStrings(log_file, 0, "", 0, 0) == 0)
	{
		blscan_error("blscan_check_param: Log file not defined.");
		return -1;
	}
	
	/* dat_file */
	if (CompareStrings(dat_file, 0, "", 0, 0) == 0)
	{
		blscan_error("blscan_check_param: Magnet strength dat file not defined.");
		return -1;
	}

	/* plt_file */
	if (CompareStrings(plt_file, 0, "", 0, 0) == 0)
	{
		blscan_error("blscan_check_param: Magnet strength plt file not defined.");
		return -1;
	}
	
	/* outer_coil_pos */
	if (blscan_param.outer_coil_pos <= 0 || blscan_param.outer_coil_pos > .03)
	{
		blscan_error("blscan_check_param: outer_coil_pos is invalid");
		return -1;
	}
	
	/* inner_coil_pos */
	if (blscan_param.inner_coil_pos <= 0 || blscan_param.inner_coil_pos > blscan_param.outer_coil_pos)
	{
		blscan_error("blscan_check_param: inner_coil_pos is invalid");
		return -1;
	}

	/* num_turns_coil */
	if (blscan_param.num_turns_coil < 1 || blscan_param.num_turns_coil > 1000)
	{
		blscan_error("blscan_check_param: num_turns_coil is invalid");
		return -1;
	}

	/* scan_start_pos */
	if (fabs(blscan_param.scan_start_pos) > .1)
	{
		blscan_error("blscan_check_param: scan_start_pos is invalid");
		return -1;
	}

	/* scan_dist */
	if (fabs(blscan_param.scan_dist) > .2)
	{
		blscan_error("blscan_check_param: scan_dist is invalid");
		return -1;
	}

	/* scan_num_samp */
	if (blscan_param.scan_num_samp < 1 || blscan_param.scan_num_samp > BLSCAN_MAX_NUM_SAMP)
	{
		blscan_error("blscan_check_param: scan_num_samp is invalid");
		return -1;
	}


	/* num_meas_ave */
	if (blscan_param.num_meas_ave < 1 || blscan_param.num_meas_ave > BLSCAN_MAX_NUM_MEAS)
	{
		blscan_error("blscan_check_param: num_meas_ave is invalid");
		return -1;
	}

	/* norm_radius */
	if (blscan_param.norm_radius < 0. || blscan_param.norm_radius > 10.)
	{
		blscan_error("blscan_check_param: norm_radius is invalid");
		return -1;
	}
	
	/* If we made it this far, all parameters have values */
	return 0;
}

/*	**************************************************************	*/
/*
 *	blscan_fact
 *	This function returns the factorial of the input number.
 *
 *	Input:
 *	n, an integer
 *
 *	Output:
 *	n_fact, n * (n-1) * (n-2) * ... * 1
 *
 *	Zachary Wolf
 *	8/20/00
 */
 
double blscan_fact(int n)
{
	/* Declare variables */
	double n_fact;
	int i;
	
	/* Compute n factorial */
	n_fact = 1.;
	for (i = n; i >= 1; i--) n_fact = n_fact * i;
	
	/* Done */
	return n_fact;
}

/*	**************************************************************	*/
/*
 *	blscan_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 blscan_log_measurement(void)
{
	/* Declare variables */
	FILE* file_ptr;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_measurement: Unable to open log file\n");
		return;
	}
	
	/* Write to the log file */
	fprintf(file_ptr, "%s  Magnetic Field Harmonics Measurement\n", TimeStr());
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_get_vt
 *	This function logs the result of an integrated voltage measurement.
 *
 *	Zachary Wolf
 *	3/17/00
 */
 
void blscan_log_get_vt(void)
{
	/* Declare variables */
	FILE* file_ptr;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_blscan_ave: Unable to open log file\n");
		return;
	}
	
	/* Write the measured values to the log file */
	fprintf(file_ptr, "%s  BLscan_get_VT just finished a forward and reverse measurement\n", TimeStr());
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_calc_vt_coeff
 *	This function logs the result of a polynomial fit calculation.
 *
 *	Input:
 *	order, order of the polynomial fit
 *	vt_coeff[0 to order], the VTi in VT(x) = VT0 + VT1*(x/R0) + VT2*(x/R0)^2 + ... + VTorder*(x/R0)^order (Vs)
 *
 *	Zachary Wolf
 *	12/3/99
 */

void blscan_log_calc_vt_coeff(int order, double coeff[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_calc_coeff: Unable to open log file\n");
		return;
	}
	
	/* Write the data to the log file */
	fprintf(file_ptr, "%s  Polynomial Coefficient Calculation\n", TimeStr());
	fprintf(file_ptr, "          Polynomial Order = %i\n", order);
	fprintf(file_ptr, "                       Coeff    \n");
	fprintf(file_ptr, "              #         (Vs)    \n");
	fprintf(file_ptr, "           -------  ------------\n");
	for (i = 0; i <= order; i++)
		fprintf(file_ptr, "           %4i     %12.5f\n", i, coeff[i]);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_calc_bln
 *	This function logs the result of a magnet strength calculation.
 *
 *	Input:
 *	num_har, number of harmonics to record
 *	bln[0 to num_har], the coeff in BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)
 *
 *	Zachary Wolf
 *	12/3/99
 */

void blscan_log_calc_bln(int num_har, double bln[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_calc_bln: Unable to open log file\n");
		return;
	}
	
	/* Write the data to the log file */
	fprintf(file_ptr, "%s  BL vs X Polynomial Coefficient Calculation\n", TimeStr());
	fprintf(file_ptr, "          BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)\n");
	fprintf(file_ptr, "          R0 = %12.6f m\n", blscan_param.norm_radius);
	fprintf(file_ptr, "          Number of harmonics = %i\n", num_har);
	fprintf(file_ptr, "          Nturns per coil = %i\n", blscan_param.num_turns_coil);
	fprintf(file_ptr, "          D1 = %f\n", blscan_param.outer_coil_pos);
	fprintf(file_ptr, "          D0 = %f\n", blscan_param.inner_coil_pos);
	fprintf(file_ptr, "                        BLn     \n");
	fprintf(file_ptr, "              #         (Tm)    \n");
	fprintf(file_ptr, "           -------  ------------\n");
	for (i = 1; i <= num_har; i++)
		fprintf(file_ptr, "           %4i     %12.5f\n", i, bln[i]);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_M_double_coil
 *	This function logs the M matrix.
 *
 *	Input:
 *	M[8][8], M matrix in the BLn calculation
 *	invM[8][8], inverse of the M matrix
 *
 *	Zachary Wolf
 *	3/24/00
 */

void blscan_log_M_double_coil(double M[8][8], double invM[8][8])
{
	/* Declare variables */
	FILE* file_ptr;
	int i, j;
	
	/* Open the log file */									 
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_calc_coeff: Unable to open log file\n");
		return;
	}
	
	/* Write the data to the log file */
	fprintf(file_ptr, "%s  M matrix for the double coil\n", TimeStr());
	for (i = 0; i < 8; i++)
	{
		fprintf(file_ptr, "          ");
		for (j = 0; j < 8; j++) fprintf(file_ptr, "%11.3e ", M[i][j]);
		fprintf(file_ptr, "\n");
	}
	fprintf(file_ptr, "%s  invM matrix\n", TimeStr());
	for (i = 0; i < 8; i++)
	{
		fprintf(file_ptr, "          ");
		for (j = 0; j < 8; j++) fprintf(file_ptr, "%11.3e ", invM[i][j]);
		fprintf(file_ptr, "\n");
	}
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_M_single_coil
 *	This function logs the M matrix.
 *
 *	Input:
 *	M[9][9], M matrix in the BLn calculation
 *	invM[9][9], inverse of the M matrix
 *
 *	Zachary Wolf
 *	3/24/00
 */

void blscan_log_M_single_coil(double M[9][9], double invM[9][9])
{
	/* Declare variables */
	FILE* file_ptr;
	int i, j;
	
	/* Open the log file */									 
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_calc_coeff: Unable to open log file\n");
		return;
	}
	
	/* Write the data to the log file */
	fprintf(file_ptr, "%s  M matrix for the single coil\n", TimeStr());
	for (i = 0; i < 9; i++)
	{
		fprintf(file_ptr, "          ");
		for (j = 0; j < 9; j++) fprintf(file_ptr, "%11.3e ", M[i][j]);
		fprintf(file_ptr, "\n");
	}
	fprintf(file_ptr, "%s  invM matrix\n", TimeStr());
	for (i = 0; i < 9; i++)
	{
		fprintf(file_ptr, "          ");
		for (j = 0; j < 9; j++) fprintf(file_ptr, "%11.3e ", invM[i][j]);
		fprintf(file_ptr, "\n");
	}
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}

/*	**************************************************************	*/
/*
 *	blscan_log_get_bln_ave
 *	This function logs the result of a magnet strength calculation.
 *
 *	Input:
 *	order, order of the polynomial fit
 *	vtn_ave[0 to order], average of VTn coefficients in fit
 *	vtn_rms[0 to order], rms variation of VTn coefficients over several measurements
 *	num_har, number of harmonics to record
 *	bln_ave[0 to num_har], average of the BL expansion coefficients for each multipole
 *	bln_rms[0 to num_har], rms variation of the BL coefficients over several measurements
 *
 *	Zachary Wolf
 *	3/22/00
 */

void blscan_log_get_bln_ave(int order, double vtn_ave[], double vtn_rms[], int num_har, double bln_ave[], double bln_rms[])
{
	/* Declare variables */
	FILE* file_ptr;
	int i;
	
	/* Open the log file */
	file_ptr = fopen(log_file, "a");
	if (file_ptr == NULL)
	{
		printf("blscan_log_get_bln_ave: Unable to open log file\n");
		return;
	}
	
	/* Write the VTn coefficients to the log file */
	fprintf(file_ptr, "%s  Average Polynomial Coefficients Over Several Measurements\n", TimeStr());
	fprintf(file_ptr, "          VT(x) = VT0 + VT1*(x/R0) + ... + VTn*(x/R0)^n (Vs)\n");
	fprintf(file_ptr, "          R0 = %12.6f m\n", blscan_param.norm_radius);
	fprintf(file_ptr, "          Polynomial Order = %i\n", order);
	fprintf(file_ptr, "                      VTn_ave       VTn_rms   \n");
	fprintf(file_ptr, "              #       (Vs/m^#)      (Vs/m^#)  \n");
	fprintf(file_ptr, "           -------  ------------  ------------\n");
	for (i = 0; i <= order; i++)
		fprintf(file_ptr, "           %4i     %12.6f  %12.6f\n", i, vtn_ave[i], vtn_rms[i]);
	
	/* Write the BLn coefficients to the log file */
	fprintf(file_ptr, "%s  Average of BLn Polynomial Coefficients Over Several Measurements\n", TimeStr());
	fprintf(file_ptr, "          BL(x) = BL1 + BL2*(x/R0) + ... + BLn*(x/R0)^(n-1) (Tm)\n");
	fprintf(file_ptr, "          R0 = %12.6f m\n", blscan_param.norm_radius);
	fprintf(file_ptr, "          Number of harmonics = %i\n", num_har);
	fprintf(file_ptr, "                      BLn_ave       BLn_rms   \n");
	fprintf(file_ptr, "              #         (Tm)          (Tm)    \n");
	fprintf(file_ptr, "           -------  ------------  ------------\n");
	for (i = 1; i <= num_har; i++)
		fprintf(file_ptr, "           %4i     %12.6f  %12.6f\n", i, bln_ave[i], bln_rms[i]);
	
	/* Close the log file */
	fclose(file_ptr);
	
	/* Done */
	return;
}