extern struct{
	double xc, yc, zc, bxc, byc, bzc;
} rzforc_;

struct FileHeader {
	int i_max;
	long nR, nZ, nSheets, magic;
	double x_min, x_max, y_min, y_max, z_min, z_max, innerRadius, outerRadius, length, dR, dZ, error;
}Seg01FileHeader, Seg02FileHeader, CapSol01CStruct, CapSol01aCStruct, CapSol01bCStruct, CapSol01cCStruct, CapSol02CStruct, TRSSol01CStruct; //TRSSol01bCStruct, TRSSol01cCStruct, TRSSol01dCStruct;

float *Seg01mapX;
float *Seg01mapY;
float *Seg01mapZ;
float *Seg01mapBx;
float *Seg01mapBy;
float *Seg01mapBz;
struct FileHeader *Seg01FileHeaderPtr = &Seg01FileHeader;

float *Seg02mapX;
float *Seg02mapY;
float *Seg02mapZ;
float *Seg02mapBx;
float *Seg02mapBy;
float *Seg02mapBz;
struct FileHeader *Seg02FileHeaderPtr = &Seg02FileHeader;
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//const char map_fname[] = "./fields/seg_1.map";
const char map_fname_list[] = "../field_map_list.txt";
//const char map_fname1[] = "/var/tmp/sam/fieldmaps/seg_1.map";
//const char map_fname2[] = "/var/tmp/sam/fieldmaps/seg_2.map";

char map_fname1[1024];
char map_fname2[1024];

void readFile(const char * filename, struct FileHeader* fh, float* mapX, float* mapY, float* mapZ, float* mapBx, float* mapBy, float* mapBz) {
#ifdef WIN32
	FILE *in = fopen(filename,"rb");
#else
	FILE *in = fopen(filename,"r");
#endif
	if (!in){
		printf("Error opening %s\n",filename);
		exit(1);
	}

	int i = 0;
	char str_temp[140];

	// Read in first 4 lines (this is a header that just gives definitions of the variables)
	for (i = 0; i <=3; i++)
	{
		fgets(str_temp, 140, in);
	}

	// Read in each element to the necessary array
	i=0;
	while(feof(in) != 1)
	{
		fscanf(in, "%f", mapX + i);
		fscanf(in, "%f", mapY + i);
		fscanf(in, "%f", mapZ + i);
		fscanf(in, "%f", mapBx + i);
		fscanf(in, "%f", mapBy + i);
		fscanf(in, "%f", mapBz + i);
		*(mapBz + i) = *(mapBz + i) * -1; // need to flip Bz to get correct direction
		i++;
	}

	// Save the minimum x, y, z values (the first element in each array)
	fh->x_min = *(mapX + 0); 
	fh->y_min = *(mapY + 0);
	fh->z_min = *(mapZ + fh->i_max -2); // z directions in reverse order - z_min is the last element (last line is blank and want the (n-1)th element
	// Save the maximum x, y, z values (the final element in each array)
	fh->x_max = *(mapX + fh->i_max -2);
	fh->y_max = *(mapY + fh->i_max -2);
	fh->z_max = *(mapZ + 0); // z directions in reverse order - z_max is the first element

	fclose(in);
	in = 0;
	
	// Rpint out first 10 lines to make sure I;ve done it right
	/*for(i = 0; i <10; i++){
		printf("x = %f", *(mapX + i));
		printf(" y = %f", *(mapY + i));
		printf(" z = %f", *(mapZ + i));
		printf(" Bx = %f", *(mapBx + i));
		printf(" By = %f", *(mapBy + i));
		printf(" Bz = %f\n", *(mapBz + i));
	}*/
	//printf("file was opened \n");
}

// Read the number of entries in the file (for malloc)
void ReadFileSize(const char * filename, struct FileHeader * fh)
{
#ifdef WIN32
	FILE *in = fopen(filename,"rb");
#else
	FILE *in = fopen(filename,"r");
#endif
	if (!in){
		printf("Error opening '%s'\n",filename);
		exit(1);
	}
	int i = 0;
	char str_temp[200];

	// Read in first 4 lines (this is a header that just gives definitions of the variables)
	for (i = 0; i <=3; i++)
	{
		fgets(str_temp, 200, in);
	}

	fh->i_max = 0;
	i = 0;

	while(feof(in) != 1)
	{
		fgets(str_temp, 200, in);
		i++;
	}
	fh->i_max = i;
	//printf("i_max = %d", i);

	fclose(in);
	in = 0;
}


void findFieldAtPoint(float x, float y, float z, float *fieldPtr, float *mapX, float *mapY, float *mapZ, float *mapBx, float *mapBy, float *mapBz, struct FileHeader* fh){
	
	double localX = x;
	double localY = y;
	double localZ = z;
	
	// Deal with -ve y positions since we only get given magnetic field values for +ve y
	double sign = 1.0;
	if(localY < 0.0) {
		sign = -sign;
		localY = -localY;
	}

	double dX = 10; // 10mm (from file - each point is separated by 10mm in each direction)
	double dY = 10;
	double dZ = 10;
	long nX = (fh->x_max - fh->x_min)/dX +1; // +1 so that we include the minimum value
	long nY = (fh->y_max - fh->y_min)/dY +1;
	long nZ = (fh->z_max - fh->z_min)/dZ +1;

	// Find the right elements for the given position
	// Find integer part of steps (want to round down) before multiplying by nX, nY etc (otherwise fractional parts might multiple up to integer parts and get the wrong element)
	int x_steps = ((localX - fh->x_min)/dX);
	int y_steps = ((localY - fh->y_min)/dY);
	int z_steps = ((fh->z_max - localZ)/dZ); // NB z is in reverse (start at 3000)
	int i = y_steps + x_steps*nY + z_steps*nY*nX;
	/*for (i = 0; i < fh->i_max; i++)
	{
		//printf("localX=%f ", localX);
		if ( *(mapX + i) <= localX && (*(mapX+i)+dX) > localX && *(mapY + i) <= localY && (*(mapY+i)+dY) > localY && *(mapZ + i) <= localZ && (*(mapZ+i)+dZ) > localZ)
			break;
	}*/
	//printf("i = %d\n", i);
	
	// if we have found the correct position (i.e. the position is in this segment)
	if(i < fh->i_max -1){
		float fx = 1.0 - (localX - *(mapX+i)) / dX;
		float fy = 1.0 - (localY - *(mapY+i)) / dY;
		float fz = 1.0 - (localZ - *(mapZ+i)) / dZ;
		//printf("localX=%f localY=%f localZ=%f dX=%lf dY=%lf dZ=%lf, nZ=%d fx=%f fy=%f fz=%f \n", localX, localY, localZ, dX, dY, dZ, nZ, fx, fy, fz);
		//printf("(Before): Bx = %f, By = %f, Bz = %f\n", *(mapBx + i), *(mapBy + i), *(mapBz + i));
		float Bx = *(mapBx + i)*fx*fy*fz + // point (i, j, k)
					*(mapBx + i + nX*nY)*fx*fy*(1.0-fz) + // point (i, j, k+1)
					*(mapBx + i + 1)*fx*(1.0-fy)*fz + // point (i, j+1, k)
					*(mapBx + i + nY)*(1.0-fx)*fy*fz + // point (i+1, j, k)
					*(mapBx + i + nY + nX*nY)*(1.0-fx)*fy*(1.0-fz) + // point (i+1, j, k+1)
					*(mapBx + i + 1 + nX*nY)*fx*(1.0-fy)*(1.0-fz) + // point (i, j+1, k+1)
					*(mapBx + i + nY + 1)*(1.0-fx)*(1.0-fy)*fz + // point (i+1, j+1, k)
					*(mapBx + i + nY + 1 + nY*nX)*(1.0-fx)*(1.0-fy)*(1.0-fz); // point (i+1, j+1, k+1)
		float By = *(mapBy + i)*fx*fy*fz + // point (i, j, k)
					*(mapBy + i + nX*nY)*fx*fy*(1.0-fz) + // point (i, j, k+1)
					*(mapBy + i + 1)*fx*(1.0-fy)*fz + // point (i, j+1, k)
					*(mapBy + i + nY)*(1.0-fx)*fy*fz + // point (i+1, j, k)
					*(mapBy + i + nY + nX*nY)*(1.0-fx)*fy*(1.0-fz) + // point (i+1, j, k+1)
					*(mapBy + i + 1 + nX*nY)*fx*(1.0-fy)*(1.0-fz) + // point (i, j+1, k+1)
					*(mapBy + i + nY + 1)*(1.0-fx)*(1.0-fy)*fz + // point (i+1, j+1, k)
					*(mapBy + i + nY + 1 + nY*nX)*(1.0-fx)*(1.0-fy)*(1.0-fz); // point (i+1, j+1, k+1)
		float Bz = *(mapBz + i)*fx*fy*fz + // point (i, j, k)
					*(mapBz + i + nX*nY)*fx*fy*(1.0-fz) + // point (i, j, k+1)
					*(mapBz + i + 1)*fx*(1.0-fy)*fz + // point (i, j+1, k)
					*(mapBz + i + nY)*(1.0-fx)*fy*fz + // point (i+1, j, k)
					*(mapBz + i + nY + nX*nY)*(1.0-fx)*fy*(1.0-fz) + // point (i+1, j, k+1)
					*(mapBz + i + 1 + nX*nY)*fx*(1.0-fy)*(1.0-fz) + // point (i, j+1, k+1)
					*(mapBz + i + nY + 1)*(1.0-fx)*(1.0-fy)*fz + // point (i+1, j+1, k)
					*(mapBz + i + nY + 1 + nY*nX)*(1.0-fx)*(1.0-fy)*(1.0-fz); // point (i+1, j+1, k+1)
		//printf("(After): Bx = %f, By = %f, Bz = %f\n", Bx, By, Bz);

		// Assign the magnetic field values at this point to the array
		*(fieldPtr) = Bx; // there will be a sign change
		*(fieldPtr + 1) = sign * By;
		*(fieldPtr + 2) = Bz;
		//printf("(After): Bx=%lf By=%lf Bz= %lf \n", *(fieldPtr), *(fieldPtr+1), *(fieldPtr+2));
	}
	else { printf("ERROR: Position not in this segment\n"); }
}

void readfieldfrommap_() {
	float field[3] = {0., 0., 0.};
	float *fieldPtr = field;
	//int printFieldToFile = 0; // if = 0 don't print to .txt file, if = 1 then print field to .txt file
	
	FILE *map_fname_list_fh = fopen(map_fname_list, "r");
	if (!map_fname_list_fh)
		{fprintf(stderr, "ERROR: failed to open %s\n",map_fname_list);exit(1);}

	if ( fgets(map_fname1, sizeof(map_fname1),map_fname_list_fh) == NULL )
		{fprintf(stderr, "ERROR: reading map_fname1 from %s\n",map_fname_list);exit(1);}
	if ( fgets(map_fname2, sizeof(map_fname2),map_fname_list_fh) == NULL )
		{fprintf(stderr, "ERROR: reading map_fname2 from %s\n",map_fname_list);exit(1);}
	fclose(map_fname_list_fh);
	// strip newline char
	int new_line = strlen(map_fname1) -1;
	if (map_fname1[new_line] == '\n') {map_fname1[new_line] = '\0';}
	new_line = strlen(map_fname2) -1;
	if (map_fname2[new_line] == '\n') {map_fname2[new_line] = '\0';}


	
	// Read in the file if it hasn't been read before
	if (Seg01mapX == 0 || Seg01mapY == 0 || Seg01mapZ == 0 || Seg01mapBx == 0 || Seg01mapBy == 0 || Seg01mapBz == 0)
	{
		// Find out how many lines of data there are
		ReadFileSize(map_fname1, Seg01FileHeaderPtr);

		// Allocate enough memory for all the files
		Seg01mapX = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));
		Seg01mapY = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));
		Seg01mapZ = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));
		Seg01mapBx = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));
		Seg01mapBy = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));
		Seg01mapBz = malloc(Seg01FileHeaderPtr->i_max * sizeof (float));

		// Read the data in
		readFile(map_fname1, Seg01FileHeaderPtr, Seg01mapX, Seg01mapY, Seg01mapZ, Seg01mapBx, Seg01mapBy, Seg01mapBz);
	}

	// Read in the file if it hasn't been read before
	if (Seg02mapX == 0 || Seg02mapY == 0 || Seg02mapZ == 0 || Seg02mapBx == 0 || Seg02mapBy == 0 || Seg02mapBz == 0)
	{
		// Find out how many lines of data there are
		ReadFileSize(map_fname2, Seg02FileHeaderPtr);

		// Allocate enough memory for all the files
		Seg02mapX = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));
		Seg02mapY = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));
		Seg02mapZ = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));
		Seg02mapBx = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));
		Seg02mapBy = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));
		Seg02mapBz = malloc(Seg02FileHeaderPtr->i_max * sizeof (float));

		// Read in the data
		readFile(map_fname2, Seg02FileHeaderPtr, Seg02mapX, Seg02mapY, Seg02mapZ, Seg02mapBx, Seg02mapBy, Seg02mapBz);
	}
	
	
	// Only try and find field at that point if it is in the segment
	if (rzforc_.xc>=Seg01FileHeaderPtr->x_min && rzforc_.xc<=Seg01FileHeaderPtr->x_max && rzforc_.yc>=(-Seg01FileHeaderPtr->y_max) && rzforc_.yc<=Seg01FileHeaderPtr->y_max && rzforc_.zc>=Seg01FileHeaderPtr->z_min && rzforc_.zc<=Seg01FileHeaderPtr->z_max)
	{
		findFieldAtPoint(rzforc_.xc, rzforc_.yc, rzforc_.zc, fieldPtr, Seg01mapX, Seg01mapY, Seg01mapZ, Seg01mapBx, Seg01mapBy, Seg01mapBz, Seg01FileHeaderPtr);
	}
	else if (rzforc_.xc>=Seg02FileHeaderPtr->x_min && rzforc_.xc<=Seg02FileHeaderPtr->x_max && rzforc_.yc>=(-Seg02FileHeaderPtr->y_max) && rzforc_.yc<=Seg02FileHeaderPtr->y_max && rzforc_.zc>=Seg02FileHeaderPtr->z_min && rzforc_.zc<=Seg02FileHeaderPtr->z_max)
	{
		findFieldAtPoint(rzforc_.xc, rzforc_.yc, rzforc_.zc, fieldPtr, Seg02mapX, Seg02mapY, Seg02mapZ, Seg02mapBx, Seg02mapBy, Seg02mapBz, Seg02FileHeaderPtr);
	}
	else
	{
	//	printf("WARNING: Coords not in magnetic field map region: %f %f %f\n", rzforc_.xc, rzforc_.yc, rzforc_.zc);
	}
	
	//printf("R=%lf Z=%lf Br=%lf Bz= %lf \n", rzforc_.rc, rzforc_.zc, field[0], field[1]);	
	rzforc_.bxc = field[0];
	rzforc_.byc = field[1];
	rzforc_.bzc = field[2];
}