/*
Description: 
Evaluate adjacency matrix using the principle of intersecting circles of intersection.

The checks that are implemented here are:

The adjacency matrices should be input by calling this function and then using < name of adj mat file.  Where the adj mat file is in human readable nauty output.

This version differs from 'EvalAdjMats.c' as follows:
1) round_numi, where i is a number, has been entered in various functions as an easier way to control the round-off numbers.  Some of the round-off numbers have been changed so that less error is introduced as a result of these numbers.  This affects the number of packings found at n = 10, and accounts for around 2 more packings being found.

2) In rare cases, for some of the triangular bipyramids, there are dihedral angles that are not defined.  This occurs when >=3 points lie on a line.  This occurs in certain packings that contain octahedra.  In these cases, the global consistency check that is normally performed and checks for the sum and differences of dihedral angles is not valid.  In 'EvalAdjMats.c' running this check when the dihedral angle was not defined caused some physical packings to be erroneously ruled out as globally inconsistent and thus unphysical.  Now, the global consistency check first checks to see that all dihedral angles are well defined - if they are, the check proceeds as normal, if they're not then a special global consistency check for 3 points lying in a line is performed.

3) All functions that check for 9 particle physical new seeds now first check to see that the distance between 2 particles, r_ij, is unknown before writing in the distance corresponding to the appropriate r_ij of that seed.  (This prevents any already known distances being overwritten).
*/

#include <stdio.h>
#include <math.h>

#define MAXLEN 85000
#define MAXLEN2 100
#define MAXLEN3 50
#define MAXLEN4 10
#define MAXLEN5 1000
#define DUMMY_VAR 1000000
//#define PI = 3.141592653589793238462643
//Pi to 1024 decimal places: //for some reason the if statements are registering errors in them when use PI defined as a global variable, so instead are defining them locally in the functions - in spherical_triangle_check and in global_consistency_check
//#define PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327887;

/*Declare local functions:*/
int get_mx_start(const int, const int);

int row_col_2_idx(const int, const int, const int);

int idx2row(const int, const int);

int idx2col(const int, const int, const int);

double law_of_cosines_num(const double, const double, const double, const int);

double check_round(const double, const double);

double integer_simplify(const double, const double, int *, int *);

void simp_from_intsimp(double *, double *, const double, const double, const double, const int, const int);

void simp_nums_denoms_from_intsimp(double [], double [], double [], double [], const int, const int);

void unique_ring_check(int (*)[MAXLEN5], int [], const int, int *);

int find_4rings_pot_closed(int [][MAXLEN2], int (*)[MAXLEN5], const int);

int find_4rings(int [][MAXLEN2], int (*)[MAXLEN5], const int);

int find_5rings(int [][MAXLEN2], int (*)[MAXLEN5], const int);

int find_6rings(int [][MAXLEN2], int (*)[MAXLEN5], const int);

int check_connections(int [][MAXLEN2], int [12], int *, int *, const int, const int, const int);

int consistency_check(double [], const int, const int, const int, const double);

double spherical_triangle_check(const double, const double, const double, const double, const double, const double);

void recalc_dihed_ang_0cancelling(double [], int *, const double, const double, const double, const double, const double, const double);

int calc_rel_dist_num(double (*)[MAXLEN2], double (*)[MAXLEN2], const int, const int, int *, const int, const int, const int, const int, const int, const int, const int);

char* calc_rel_dist_analytic(int [][MAXLEN2], int [][MAXLEN2], char **, const int, const int);

int check_num_ints_of_2_int_circs(int [][MAXLEN2], const int, const int, const int, int *);

int check_num_points_on_int_circ(int [][MAXLEN2], int [][MAXLEN2], char **, double (*)[MAXLEN2], const int, const int, const int, int *);

int check_num_ints_of_m_int_circs(int [][MAXLEN2], int [][MAXLEN2], char **, double (*)[MAXLEN2], const int, const int, const int, int *);

int check_closed_5_rings(int [][MAXLEN2], const int, const int, int *);

int check_opp_points_of_4_ring_dont_touch(int [][MAXLEN2], const int, const int, int *);

void new_seed_8_particles(int [][MAXLEN2], char **, double (*)[MAXLEN2], const int);

int connected_4ring_check1(int [][MAXLEN2], int (*)[MAXLEN5], const int, const int);

void new_seed_9part_graph8901(int [][MAXLEN2], int (*)[MAXLEN5], char **, double (*)[MAXLEN2], const int, const int);

void new_seed_9part_graph13380(int [][MAXLEN2], int (*)[MAXLEN5], char **, double (*)[MAXLEN2], const int, const int);

void new_seed_9part_graph8926(int [][MAXLEN2], int [][MAXLEN2], int (*)[MAXLEN5], char **, double (*)[MAXLEN2], const int, const int);

void new_seed_9part_graph12811(int [][MAXLEN2], int (*)[MAXLEN5], char **, double (*)[MAXLEN2], const int, const int);

int new_seed_9part_graph8974(int [][MAXLEN2], int (*)[MAXLEN5], const int, const int);

int pot_closed_4ring_checks(int [][MAXLEN2], char **, double (*)[MAXLEN2], const int, const int, int *);

void new_seed_4ring_checks(int [][MAXLEN2], char **, double (*)[MAXLEN2], const int);

void new_seed_5ring_checks(int [][MAXLEN2], int [][MAXLEN2], char **, double (*)[MAXLEN2], const int);

int new_seed_6ring_checks(int [][MAXLEN2], char **, double (*)[MAXLEN2], const int, const int, int *);

int calc_new_rel_dists(int [][MAXLEN2], int [][MAXLEN2], char **, double (*)[MAXLEN2], int *, int *, const int, const int, const int, int *);

int global_consistency_check(double [], const int, const int, int *);

/*Start of main function:*/
main() {

	int c, flag, row_count, col_count, first, iter_count, max_col, i, j, k, num_int_circs, int_col_count, violated, mat_count, phys_count, unphys_count, only_print_stats, solved, row, col, print_adj_mats_with_sols, only_print_sols, print_all_adj_mats, num_violations, print_unphys_mats, mult_num_sol_count, valid_sols_count, contact_count, unknown, only_record_known_sols;
	int a[MAXLEN2][MAXLEN2], int_mat[MAXLEN2][MAXLEN2], violations[MAXLEN2];  //without using pointers here, 9 particles is the maximum number that this program can evaluate.  Later, I'll rewrite with pointers to evaluate > 9.
	int phys_mats[MAXLEN], unphys_mats[MAXLEN]; //once again, the maximum length here is for now dictated by 9 particles being the maximum amount that we'll evaluate.
	double num_sol[MAXLEN5][MAXLEN2]; //this way, we're only saving 100 potential solutions in num_sol, I think this should be enough, but theoretically there can be 2^(n(n-1)/2 - (3n-6)) solutions to try.
	char *sol[MAXLEN2]; //if I add in the analytical calculation of the formulas, will then presumably need to add another dimension to this, which might be tricky - I'm not quite sure how to do that
	
	/*Initialize local variabels*/
	flag = 0;
	first = 0;
	row_count = 0;
	col_count = 0;
	max_col = -1;
	iter_count = 0;
	int_col_count = 0;
	mat_count = 0;
	phys_count = 0;
	unphys_count = 0;
	num_violations = 11;
	only_print_stats = 1; //set to 1 if only want to print all answers, physical and unphysical, without any of the matrices, set to 0 if want everything to print out
	print_adj_mats_with_sols = 1; //set to 1 if want the adjacency matrices associated with solutions to be printed out, set to 0 if not
	print_all_adj_mats = 0; //set to 1 if want all adjacency matrices to be printed, set to 0 if don't
	only_print_sols = 1; //set to 1 if only want information about matrices with solutions to print out, otherwise set to 0
	print_unphys_mats = 0; //set to 1 if you want the numbers corresponding to unphysical matrices printed out, set to 0 if not 
	only_record_known_sols = 0; //set to 1 if only want to record solutions where all relative distances are known (i.e. not including ones that haven't violated any rule but also don't have all distances solved)
	mult_num_sol_count = 1;
	valid_sols_count = 0;
	for (i = 0; i < num_violations; ++i) {
		violations[i] = 0;
	}
	
	while ((c = getchar()) != EOF) {
		if (c == 'G') { //this is the start of a new adjacency matrix (specifically the header for a new adjacency matrix)
			if (only_print_sols == 0)
				printf("\n");
			++first;
			flag = 1;
			++mat_count;
		}
		else if (c == '.') {
			if (flag == 1) {
				if (only_print_sols == 0)
					printf(".\n");
			}
			flag = 0;
		}
		if (flag == 1) {//print header of each adjacency matrix
			if (only_print_sols == 0) //only print the header if we want to print information about every adjacency matrix
				putchar(c);
		}
		else if (flag == 0) { //save the current adjacency matrix read in the input file to the temporary adjacency matrix, a, which is an array.
			if (c == '1' || c == '0') {
				if (c == 49)
					a[row_count][col_count] = 1;
				if (c == 48)
					a[row_count][col_count] = 0;
				++col_count;
				iter_count = 1;
			}
			if (c == '\n') {
				if (iter_count == 1) {
					++row_count;
					max_col = col_count;
					col_count = 0; //re-initialize
				}
			}
		}
		if (row_count == max_col) { //if we've finished copying the entire adjacency matrix, we can print it out or evaluate it.
			/*Print the adjacency matrix out, if desired*/
			if (print_all_adj_mats != 0) {
				for (i = 0; i < max_col; ++i) {
					for (j = 0; j < max_col; ++j) {
						if (j < max_col - 1)
							printf("%d",a[i][j]);
						else
							printf("%d\n",a[i][j]);
					}
				}
			}
			/*Construct the intersection circle matrix*/
			for (i = 0; i < max_col; ++i) {
				for (j = i+1; j < max_col; ++j) {
					if (a[i][j] == 1) {
						for (k = 0; k < max_col; k++) {
							if (k != i && k != j) { //don't bother to check whether the adj mat entries are 1 if we're on the same row or column as the int circ
								if (a[k][i] == 1 && a[k][j] == 1)
									int_mat[k][int_col_count] = 1;
								else
									int_mat[k][int_col_count] = 0;
							}
						}
					}
					else
						for (k = 0; k < max_col; ++k) {
							int_mat[k][int_col_count] = 0;
						}
					++int_col_count;
				}
			}
			/*Print the intersection circle matrix, if desired*/
			num_int_circs = max_col*(max_col-1)/2;
			if (only_print_stats == 0) {
				for (i = 0; i < max_col; ++i) {
					for (j = 0; j < num_int_circs; ++j) {
						if (j < num_int_circs - 1)
							printf("%d",int_mat[i][j]);
						else
							printf("%d\n",int_mat[i][j]);
					}
				}
			}
			
			/*Check physicality of intersection circle matrix by checking for violations of the rules constructed thus far.*/
			violated = 0; //initialize
			/*Initialize the solution vectors.*/
			for (i = 0; i < num_int_circs; ++i) {
				row = idx2row(i,max_col)-1;
				col = idx2col(i,row+1,max_col)-1;
				if (a[row][col] == 1) {
					sol[i] = "R";
					num_sol[0][i] = 1;
				}
				else {
					sol[i] = "Unknown";
					num_sol[0][i] = DUMMY_VAR;
				}
			}
			
			/*CHECK 1*/
			/*Check that no 2 intersection circles intersect at > 2 points*/
			violated = check_num_ints_of_2_int_circs(int_mat, max_col, num_int_circs, only_print_sols, violations);
			
			/*CHECK 2*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				/*Check that no more than 5 points lie on any 1 intersection circle.  And if 3-5 points lie on an intersection circle, check that all of these points don't touch eachother.*/
				violated = check_num_points_on_int_circ(a, int_mat, sol, num_sol, max_col, num_int_circs, only_print_sols, violations); //note don't have to put & in front of the pointers to the arrays, as an array 																																			already implies the address of that object, whereas an int doesn't (and if I had 																																			an int that was a pointer here, it would have to appear as &int - likewise, 																																		passing an int pointer into the local function requiers calling it as *int, but 																																		this is not necessary for an array).
			}
			
			/*CHECK 3*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				/*Check that 6 or more intersection circles do not intersect at more than 1 point.  Also, check for 2 intersections of 4 or 5 intersection circles, and impose the relevant criteria for them, check 					whether they're violated, and if not, record the appropriate solutions to the unknown relative distances.*/
				violated = check_num_ints_of_m_int_circs(a, int_mat, sol, num_sol, max_col, num_int_circs, only_print_sols, violations);
			}
			
			/*CHECK 4*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				/*Check that a closed 5 ring does not surround a circle of intersection.*/
				//printf("\nGraph %d, closed 5 ring check:\n",mat_count);
				violated = check_closed_5_rings(a, max_col, only_print_sols, violations);
			}
			
			/*CHECK 5*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				/*Check that 2 points connected to opposite sides of a closed 4 ring by 2 bonds do not touch.*/
				violated = check_opp_points_of_4_ring_dont_touch(a, max_col, only_print_sols, violations);
			}
			
			/*CHECK 6*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if ((violated == 0) && (max_col >= 8)) {
				/*Check for 8 particle new seeds.*/
				new_seed_8_particles(a, sol, num_sol, max_col);
			}
			
			/*CHECK 7*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				//printf("\nGraph %d, connected 4 ring check:\n",mat_count);
				/*Check that no violations involving potentially closed off 4 rings occur (i.e. 4 rings whose crossterms can = R).*/
				violated = pot_closed_4ring_checks(a, sol, num_sol, max_col, only_print_sols, violations);
			}
			
			/*CHECKs 8*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if ((violated == 0) && (max_col >= 9)) {
				//printf("\nGraph %d, n >= 9 new seed checks:\n",mat_count);
				/*Check for new seeds of >= 9 particles involving closed 6 rings.*/
				violated = new_seed_6ring_checks(a, sol, num_sol, max_col, only_print_sols, violations); //I've put the 6 ring check first, as it can violate a matrix, and thus make the other checks unecessary
				if (violated == 0) {
					/*Check for new seeds of >= 9 particles involving closed 5 rings.*/
					new_seed_5ring_checks(a, int_mat, sol, num_sol, max_col);
					/*Check for new seeds of >= 9 particles involving closed 4 rings.*/
					new_seed_4ring_checks(a, sol, num_sol, max_col);
				}
			}
			
			/*CHECK 9*/
			/*If we haven't already violated one of the rules, then check another rule*/
			if (violated == 0) {
				/*Calculate all new relative distances, if possible - and check to see that no distance is < R.*/
				//printf("\nGraph %d:\n",mat_count);
				violated = calc_new_rel_dists(a, int_mat, sol, num_sol, &mult_num_sol_count, &valid_sols_count, max_col, num_int_circs, only_print_sols, violations);
			}
			
			/*Document which adjacency matrices violate the rules, and are thus unphysical*/
			if (violated == 1) {
				if (print_unphys_mats != 0)
					unphys_mats[unphys_count] = mat_count;
				++unphys_count;
			}
	
			/*If graph did not violate any existing rules, then document that and fill in all known relative distance solutions*/
			if (violated == 0) {
				unknown = 0; //initialize, in case want to record all solutions - this will be only place it is set
				if (only_record_known_sols == 1) {
					//for (i = 0; i < valid_sols_count; ++i) { //for now will leave this out
						unknown = 0; //initialize
						for (j = 0; j < num_int_circs; ++j) { //Determine if a non-violated solutuion has unknown distances (if a distance is unknown it will be initialized to 10^6)
							if (num_sol[0][j] == DUMMY_VAR) {
								unknown = 1; //set to 1 if any distance is unknown
								break;
							}
						}
					//}
				}
				if (unknown == 0) { //only record if all distances are known or if want to record all potential solutions (including ones with unkown distances)
					if (only_print_stats != 0 && only_print_sols != 0) { //print header if we haven't already printed it
						printf("\n\n\nGraph %d, order %d.\n",mat_count,max_col);
					}
					if (print_adj_mats_with_sols != 0) { //print adjacency matrix for those that have physical solutions, if desired
						for (i = 0; i < max_col; ++i) {
							for (j = 0; j < max_col; ++j) {
								if (print_all_adj_mats == 0) { //if the intersection circle matrix is being printed out also, then print out a new line with an adjacency matrix header
									if (i == 0 && j == 0)
										printf("\nAdjacency Matrix:\n");
								printf("%d",a[i][j]);
								if (j == max_col - 1)
									printf("\n");
								}
							}
						}
					}
					printf("\nGraph did not violate any of the established rules.\n");
					for (j = 0; j < num_int_circs; ++j) { //print relative distances, both known and unknown (if a distance is unknown it will be written as "r[i][j] = unknown")
						if (j == 0)
							printf("\n");
						row = idx2row(j,max_col);
						col = idx2col(j,row,max_col);
						printf("r[%d][%d] = ",row,col);
						printf("%s\n",sol[j]);
					}
					for (i = 0; i < valid_sols_count; ++i) {
						if (i == 0)
							printf("\nThe Numerical Approximations of the Analytical Solutions:\n(Here, just the coefficients are written, and all solutions should be taken to be *R)\n");
						contact_count = 0; //initialize contact count
						for (j = 0; j < num_int_circs; ++j) { //print the numerical approximations of the analytical relative distances, both known and unknown (if a distance is unknown it will be initialized to 10^6)
							if ((j == 0) && (valid_sols_count > 1))
								printf("\nSolution %d:\n",i+1);
							row = idx2row(j,max_col);
							col = idx2col(j,row,max_col);
							printf("r[%d][%d] = %f\n",row,col,num_sol[i][j]);
							if (num_sol[i][j] == 1) {
								++contact_count; //document the number of contacts the structure makes
							}
						}
						//printf("\nNumber of contacts = %d\n",contact_count);
						if (contact_count > 3*max_col-6)
							printf("\nThis structure has greater than 3n-6 contacts: it has %d contacts\n",contact_count);
					}
					phys_mats[phys_count] = mat_count;
					++phys_count;
				}
			}
			row_count = 0; //re-initialize
			col_count = 0;
			iter_count = 0;
			int_col_count = 0;
			mult_num_sol_count = 1;
			valid_sols_count = 0;
		}
	}
	
	/*Print how many matrices where physical and unphysical, and which fell under each category*/
	printf("\nFinal Statistics:\n\n%d/%d matrices were unphysical.\n%d/%d matrices were physical.\n",unphys_count,mat_count,phys_count,mat_count);
	if (print_unphys_mats != 0) {
		printf("\nUnphysical Matrices: ");
		for (i = 0; i < unphys_count; ++i) {
			if (i < unphys_count - 1)
				printf("%d, ",unphys_mats[i]);
			else
				printf("%d.\n",unphys_mats[i]);
		}
	}
	printf("\nPhysical Matrices: ");
	for (i = 0; i < phys_count; ++i) {
		if (i < phys_count - 1)
			printf("%d, ",phys_mats[i]);
		else
			printf("%d.\n",phys_mats[i]);
	}
	
	/*Print how many times each rule was violated*/
	printf("\n%d matrices were unphysical because they implied that 2 or more intersecting circles intersect > 2 times.\n",violations[0]);
	printf("%d matrices were unphysical because they implied that > 5 points lie on an intersection circle.\n",violations[1]);
	printf("%d matrices were unphysical because they implied that there were > 2 points on an intersection circle that all touched eachother.\n",violations[2]);
	printf("%d matrices were unphysical because they implied that 4 intersection circles intersect at 2 points a distance R apart.\n",violations[3]);
	printf("%d matrices were unphysical because they implied that 5 intersection circles intersect at 2 points a distance R apart.\n",violations[4]);
	printf("%d matrices were unphysical because they implied that > 5 intersection circles intersect at 2 points.\n",violations[5]);
	printf("%d matrices were unphysical because they implied that a closed 5 ring surrounds a circle of intersection.\n",violations[6]);
	printf("%d matrices were unphysical because they implied that 2 points on opposite sides of a closed 4 ring touch.\n",violations[7]);
	printf("%d matrices were unphysical because they had no solution or they implied that a distance was < R.\n",violations[8]);
	printf("%d matrices were unphysical because they contained the 9 particle unphysical new seed (graph 8974), which implies a distance < R.\n",violations[9]);
	printf("%d matrices were unphysical because points in connected 4 rings touched eachother that implied a distance < R.\n",violations[10]);
	

} //end of main

/*Local functions:*/

/*
Description: 
Gives the starting index for row i of a symmetric nxn adjacency matrix.  Currently does this for a 0 based counting system.  I.e., the matrix element numbers are 0 based, going from 0 to n*(n-1)/2 - 1, and fill up only the upper off-diagonal triangle of the nxn matrix.  The row number, i, is 1 based, going from 1 to n.  So, inputing the 1 based row number, i, this returns the 1st element number in that row, where the element number is 0 based and only the upper off-diagonal triangle of the symmetric matrix is numbered.

Output: the 1st index number of a row
*/
int get_mx_start(const int i, const int n) {
	return (i-1)*n-(i*(i-1)/2); //add +1 here if want this to be for a 1 based counting system
}

/*
Description:
Gives the corresponding matrix index of a symmetric nxn matrix.  The matrix elements are 0 based, going from 0 to n*(n-1)/2 - 1, and filling up only the off-diagonal upper triangle of the matrix (as this is a symmetric adjacency matrix, and so the diagonal elements = 0 and the bottom triangle elements = top triangle elements).  It accepts the corresponding row and column of the nxn matrix, which are both 1 based, and it accepts the number of particles, n (also 1 based).  Thus, it accepts the 1 based row and column and returns the 0 based symmetric matrix element.

Note that at the time I put in the 9 particle new seed rules - I finally updated this function to check itself that the row < col, and switch if it isn't - thus, from this point on, we no longer have to input row < col when calling this function (and I can CHANGE all parts of the code that do this) - I can simply enter in any order row, col into the function - and it will parse out here row < col, and always only return the upper 0ff-diagonal 0 based matrix element (accepting a 1 based row and column element).

Output: the 0-based symmetric off-diagonal matrix index
*/
int row_col_2_idx(const int row_enter, const int col_enter, const int n) {

	//declare local variables
	int idx, row, col, start_row;
	
	//enter in below with row < col
	if (row_enter < col_enter) {
		row = row_enter;
		col = col_enter;
	}
	else {
		row = col_enter;
		col = row_enter;
	}
	
	start_row = get_mx_start(row,n);
	idx = start_row + col - (row + 1);
	return idx;
}

/*
Description: 
Gives the corresponding row of a given matrix index, i, of a symmetric nxn adjacency matrix.
Currently does this for a 0 based counting system, as C uses.  So, this accepts a zero based symmetric matrix indice and returns a 1 based row indice.

Output: the row
*/
int idx2row(const int i, const int n) {

	//declare and initialize local variables
	int found, row, start_row, start_next_row;
	found = 0;
	row = 0;

	while (found == 0) {
		++row;
		if (row == 1)
			start_row = get_mx_start(row,n);
		start_next_row = get_mx_start(row+1,n);
		if (i >= start_row && i < start_next_row)
			found = 1;
		else
			start_row = start_next_row;
	}
	
	return row;
}

/*
Description: 
Gives the corresponding column of a given a symmetric matrix index, i, whose corresponding row is j, and whose total dimensions are nxn (but whose elements are only numbered in the off-diagonal upper triangle).  
Currently does this for a 0 based counting system, as C uses.  So, this accepts a zero based matrix indice and a 1 based row indice and returns a 1 based col indice.

Output: the column
*/
int idx2col(const int i, const int j, const int n) {

	//declare local variable
	int plus_num;
	
	plus_num = i - get_mx_start(j, n);
	return plus_num + j + 1;
}

/*
Description:
Accepts 3 relative distances that comprise a triangle, and returns the angle opposite one of those relative distances.  The last input, indicates which angle should be returned - i.e. 1 would imply that the angle opposite the first relative distance should be returned, and so on for 2 and 3.

Output: The angle opposite one of the relative distances in a triangle.
*/
double law_of_cosines_num(const double r12, const double r13, const double r23, const int pick) {
	
	/*declare local variables*/
	double angle, angle_arg, numerator;
	
	if (pick == 1) {
		numerator = pow(r13,2) + pow(r23,2) - pow(r12,2);
		angle_arg = numerator/(2*r13*r23); //the angle opposite r12
	}
	if (pick == 2) {
		numerator = pow(r12,2) + pow(r23,2) - pow(r13,2);
		angle_arg = numerator/(2*r12*r23); //the angle opposite r13
	}
	if (pick == 3) {
		numerator = pow(r12,2) + pow(r13,2) - pow(r23,2);
		angle_arg = numerator/(2*r12*r13); //the angle opposite r23
	}
	
	angle_arg = check_round(angle_arg,numerator);
	
	angle = acos(angle_arg);
	
	return angle;
}

/*
Description: Accepts a double as an input, and if it's above 1 in the 7th decimal place or higher, it will round the number to 1 exactly.

Output: 1, if the number has been rounded, the input number arg, if it has not been rounded.
*/
double check_round(const double arg, const double numerator) {
	//printf("\nCheck round was called!\n");
	/*Declare and initialize local variables*/
	int i, round, prec1, prec2;
	double num, rem;
	char *string_num;
	round = 1;
	num = arg;
	//printf("num = %f\n",arg);
	prec1 = -10; //Testing different values, it seems as long as this is -5 or higher (I've tested up to -16), it doesn't matter.  Originally set to -5, this is the precision to which we check if the numerator = 0, if it is > 10^(prec1) - i.e. greater than 0 to a precision of -prec1, then we don't round the quotient to 0
	prec2 = -10; //Testing different values, it seems that -10 is the optimal value.  Originally set to -6, this is the precision to which we check if the remainder = 0, if it is > 10^(prec2) - i.e. greater than 0 to a precision of -prec2, then we don't round
	
	//printf("numerator rem = %f\n",numerator);
	if (fabs(numerator) > pow(10,prec1)) { //if the numerator is not within 10^prec1 of 0, then don't round the quotient to 0
		round = 0;
		//printf("\nI'm greater than 10^(-5) within 0\n");
	}
	if (round == 1) {
		//printf("\nRound to 0 was called!\n");
		num = 0;
	}
	if (round == 0) { //if haven't already rounded, then check for one rounding (on 6/02/09, changed this function - switched ordering, and rounded to zero in all cases...)
		round = 1; //re-initialize
		if (isinf(fabs(arg)) == 0) { //if the argument is not infinite, and we're in the scenario where we might have to 'round' to 1, do this check
			rem = fabs(fabs(num) - 1);
			//printf("rem = %f\n",rem);
			if (fabs(rem) > pow(10,prec2)) //don't round in this case
				round = 0;
			if (round == 1) {
				//printf("\nRound to 1 was called!\n");
				if (arg > 0)
					num = 1;
				else
					num = -1;
			}
			/*if (round == 0) { //check for any integer rounding -- i.e. if num = 2.00004, then round to 2
				i = 0; //initialize
				while (fabs(arg) > (double)i-1) {
					//printf("\nfabs(arg) = %f, i = %d\n",fabs(arg),i);
					round = 1; //re-initialize
					rem = fabs(fabs(num) - i);
					if (fabs(rem) > pow(10,-5)) //don't round in this case
						round = 0;
					if  (round == 1) {
						if (arg >= 0)
							num = i;
						else
							num = -i;
					}
					++i; //increment i
				}
			}*/
		}
	}
	
	//printf("round = %d, num = %f\n",round,num);
	
	return num;

}

/*
Description
This function accepts a numerator and a denominator, and checks whether the numerator is an integer multiple of the denominator.  If it is, it returns the common integer value between the
numerator and the denominator -- so that the simplified(numerator/denominator) would then by (numerator/common_value)/(denominator/common_value).  If there is no common value, it returns 1.

Inputs:
	the numerator
	the denominator
	'multiple' -- this is returned as 1 if the numerator is an integer multiple of the denominator, or vice-versa, and 0 otherwise.  If *multiple is returned back as 1, then do
		numerator/denominator = common_val or numerator/denominator = 1/common_val (depending on whether num > denom or num < denom, respectively).  If *multiple is returned back as 0,
		then just do numerator/denominator = (numerator/common_val)/(denominator/common_val)
	'num_greater' -- this is returned as 1 if the numerator is greater than the denominator, as 0 if the numerator is less than the denominator, and as -1 if the numerator is = the denominator.

Ouptut: The common integer value between the numerator and the denominator (defaulting to 1).
*/
double integer_simplify(const double num, const double denom, int *multiple, int *num_greater) {

	/*Declare and Initialize Local Variables*/
	int i;
	double common_val, temp;
	common_val = 1;
	*multiple = 0;
	*num_greater = -1;
	
	/*Check if numerator = integer multiple of denominator or vice-versa*/
	if (fabs(round(num*pow(10,5))) == fabs(round(denom*pow(10,5)))) { //I won't return this as *multiple = 1, because here I will want to do both num/common_val and denom/common_val
		common_val = num;
	}
	else if (fabs(round(num*pow(10,5))) > fabs(round(denom*pow(10,5)))) {
		temp = num/denom;
		//printf("\nIn Simp Func2:\nnum/denom = %f,floor(num/denom) = %f,abs of rounded difference = %f\n", temp,floor(round(temp*pow(10,5)))/pow(10,5),fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))));
		if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if the num/denom is an integer -- if it is, then num/denom = floor(num/denom), otherwise num/denom > floor(num/denom)
			common_val = temp;
			*multiple = 1;
			*num_greater = 1;
		}
	}
	else if (fabs(round(num*pow(10,5))) < fabs(round(denom*pow(10,5)))) {
		temp = denom/num;
		if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if the denom/num is an integer -- if it is, then denom/num = floor(denom/num), otherwise denom/num > floor(denom/num)
			common_val = temp;
			*multiple = 1;
			*num_greater = 0;
		}
	}
	
	/*Check if numerator and denominator share a common integer divisor*/
	if (common_val == 1) { //only do this check if numerator is not an integer multiple of the denominator, or vice versa.
		if (fabs(round(num*pow(10,5))) > fabs(round(denom*pow(10,5)))) {
			*num_greater = 1;
			for (i = floor(round(num*pow(10,5)))/pow(10,5); i >= 2; --i) { //start with integer nearest to, but not > numerator and increment downwards to 2, looking for an integer value in common -- I determine the floor in this way, by rounding *10^5 and then dividing by 10^5 to avoid it giving a floor of 3 if it = 3.99999999, for example.
				temp = num/i;
				if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if the numerator/i is an integer -- if it is, then numerator/i = floor(numerator/i), otherwise numerator/i > floor(numerator/i)
					if (i <= denom) //only bother to check if i is a common divisor of denom if it is a common divisor of num
						temp = denom/i;
					else
						temp = i/denom;
					if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if i is a common divisor of the denominator by checking if denominator/i or the inverse (as appropriate) is an integer
						common_val = i;
						break; //as I'm working back from the greatest possible common divisor, once I hit a common divisor, I'm guarranteed it's the largest one, and can stop searching for it
					}
				}
			}
		}
		else if (fabs(round(num*pow(10,5))) < fabs(round(denom*pow(10,5)))) {
			*num_greater = 0;
			for (i = floor(round(denom*pow(10,5)))/pow(10,5); i >= 2; --i) { //start with integer nearest to, but not > denominator and increment downwards to 2, looking for an integer value in common -- I determine the floor in this way, by rounding *10^5 and then dividing by 10^5 to avoid it giving a floor of 3 if it = 3.99999999, for example.
				temp = denom/i;
				if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if the denominator/i is an integer -- if it is, then denominator/i = floor(denominator/i), otherwise denominator/i > floor(denominator/i)
					if (i <= num) //only bother to check if i is a common divisor of num if it is a common divisor of denom
						temp = num/i;
					else
						temp = i/num;
					if (fabs(round(temp*pow(10,5)) - floor(round(temp*pow(10,5)))) < 1) { //check if i is a common divisor of the numerator by checking if numerator/i or the inverse (as appropriate) is an integer
						common_val = i;
						break; //as I'm working back from the greatest possible common divisor, once I hit a common divisor, I'm guarranteed it's the largest one, and can stop searching for it
					}
				}
			}
		}
	}
		
	return common_val;
	
}

/*
Description:
This function takes the output of integer_simplify and produces the simplified numerator and denominator.

Ouptut:
	Technically, this function has no output - but the following two values are the results of this function, and they are both updated as a pointer to a double:
	simp_num: the simplified numerator.
	simp_denom: the simplified denominator.
*/
void simp_from_intsimp(double *simp_num, double *simp_denom, const double num, const double denom, const double common_val, const int mult, const int num_greater) {
	
	
	if (mult == 1) { //if numerator is an integer multiple of denominator or vice-versa, then num/denom = common_val or num/denom = 1/common_val depending on if num > or < denom, respectively.
		if (num_greater == 1) {
			*simp_num = common_val;
			*simp_denom = 1;
		}
		else if (num_greater == 0) {
			*simp_num = 1;
			*simp_denom = common_val;
		}
	}
	else { //if mult = 0 (or simply if not = 1) then simp_num = num/common_val and simp_denom = denom/common_val
		*simp_num = num/common_val; //note if common_val = 1, this will simply return the numerator and denominator unaltered
		*simp_denom = denom/common_val;
	}
	
}

/*
Description:
This function accepts an array of numerators and denominators and returns arrays of simplified numerators and denominators as simplified by integer_simplify (both in straight integer and
sqrt(integer) modes).
As each term in entered for the numerator and denominator is simplified individually, for n numerator terms and m denominator terms, we conduct
for i = 1:n
	for j = 1:m
		simplify(num_i, denom_j);
		update num_i and denom_j values to their respective simplified values;
	end
end

Inputs:
	...
	n_nums: the number of numerators that are individually simplifying.
	n_denoms: the number of denominators that are individually simplifying.

Output: Technically, there is no output, but the simplified numerators and denominators are updated in the simp_nums and simp_denoms arrays.
*/
void simp_nums_denoms_from_intsimp(double simp_nums[], double simp_denoms[], double nums[], double denoms[], const int n_nums, const int n_denoms) {

	/*Declare and Initialize Local Variables*/
	int i, j, mult, num_greater;
	double comm_val;
	mult = 0;
	num_greater = 0;
	for (i = 0; i < n_nums; ++i) {
		simp_nums[i] = nums[i];
	}
	for (i = 0; i < n_denoms; ++i) {
		simp_denoms[i] = denoms[i];
	}
	
	/*Main Loop*/
	for (i = 0; i < n_nums; ++i) {
		for (j = 0; j < n_denoms; ++j) {
			comm_val = integer_simplify(simp_nums[i], simp_denoms[j], &mult, &num_greater); //First do integer simplification
			//printf("\nIn Simplify Func1:\ncomm_val = %f, in_num = %f, in_denom = %f, multiple = %d, num_greater = %d\n",comm_val,simp_nums[i],simp_denoms[j],mult,num_greater);
			simp_from_intsimp(&simp_nums[i], &simp_denoms[j], simp_nums[i], simp_denoms[j], comm_val, mult, num_greater);
			comm_val = integer_simplify(pow(simp_nums[i],2), pow(simp_denoms[j],2), &mult, &num_greater); //Then do sqrt(integer) simplification
			simp_from_intsimp(&simp_nums[i], &simp_denoms[j], simp_nums[i], simp_denoms[j], sqrt(comm_val), mult, num_greater); //if using the sqrt(integer) simplification then should feed in sqrt(common_val) as common_val
		}
	}
	

}

/*
Description:
This function accepts an array of m rings (where m is 4, 5, or 6 - the array is thus of dimensions [4][MAXLEN5], [5][MAXLEN5], or [6][MAXLEN5]), and tests whether the 'test_ring' is already contained in the array (potentially in a different order).  If it is already contained, it is not recorded into the array; if it isn't, then it is documented.  (Note that I only wrote this function when it became apparent that I was going to use it alot; therefore, some functions that could use this, don't, and instead have their own check - such as the 8 particle new seed check, for example - I can go through later an dreplace these to use this function, if desired).

Output: No output, but the array containing all of the m rings is updated if the test ring is unique, and the integer documenting the number of rings is similarly updated.
*/
void unique_ring_check(int (*rings)[MAXLEN5], int test_ring[], const int num_points, int *num_rings) {
	
	/*Declare local variables*/
	int i, j, k, document, count;
	document = 1;
	
	for (i = 0; i < *num_rings; ++i) { //loop through all rings documented up until now
		count = 0; //initialize
		for (j = 0; j < num_points; ++j) { //loop through all points in the rings
			for (k = 0; k < num_points; ++k) { //loop through all points in the test ring that we are checking for uniqueness
				if (rings[j][i] == test_ring[k])
					++count;
			}
		}
		if (count == num_points) { //if all points in the test ring are equal to all points in the ith ring, albeit in a potentially different order, then it is equal to the ith ring
			document = 0;
			break;
		}
	}
	
	if (document == 1) { //if the test ring is unique, document it
		for (i = 0; i < num_points; ++i)
			rings[i][*num_rings] = test_ring[i];
		++*num_rings;
	}

}

/*
Description:
Given an adjacency matrix, a, find all unique 4 rings that are contained within it; where here, the crossterms in the 4 rings can also potentially be in contact (i.e. a_ij = 1, a_jk = 1, a_kp = 1, a_pi = 1; with no equirements on a_ik and a_jp, i.e. the crossterms can touch or not touch, it doesn't matter).

Output: The number of unique 4 rings.  The array containing all of the 4 rings is also updated.
*/
int find_4rings_pot_closed(int a[][MAXLEN2], int (*rings4)[MAXLEN5], const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, num_rings;
	int test_ring[4];
	num_rings = 0;
	
	for (i = 0; i < max_col; ++i) {
		for (j = 0; j < max_col; ++j) {
			if (a[i][j] == 1) {
				for (k = 0; k < max_col; ++k) {
					if (k != i) {
						if (a[j][k] == 1) {
							for (p = 0; p < max_col; ++p) {
								if (p != j) {
									if ((a[k][p] == 1) && (a[p][i] == 1)) {
										test_ring[0] = i;
										test_ring[1] = j;
										test_ring[2] = k;
										test_ring[3] = p;
										unique_ring_check(rings4, test_ring, 4, &num_rings);
									}
								}
							}
						}
					}
				}
			}
		}
	}
	
	return num_rings;

}

/*
Description:
Given an adjacency matrix, a, find all unique closed 4 rings that are contained within it (i.e. a_ij = 1, a_jk = 1, a_kp = 1, a_pi = 1; a_ik = a_jp = 0 - crossterms don't touch).

Output: The number of unique closed 4 rings.  The array containing all of the closed 4 rings is also updated.
*/
int find_4rings(int a[][MAXLEN2], int (*closed_4rings)[MAXLEN5], const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, num_rings;
	int test_ring[4];
	num_rings = 0;
	
	for (i = 0; i < max_col; ++i) {
		for (j = 0; j < max_col; ++j) {
			if (a[i][j] == 1) {
				for (k = 0; k < max_col; ++k) {
					if (k != i) {
						if ((a[j][k] == 1) && (a[i][k] == 0)) {
							for (p = 0; p < max_col; ++p) {
								if (p != j) {
									if ((a[k][p] == 1) && (a[p][i] == 1) && (a[p][j] == 0)) {
										test_ring[0] = i;
										test_ring[1] = j;
										test_ring[2] = k;
										test_ring[3] = p;
										unique_ring_check(closed_4rings, test_ring, 4, &num_rings);
									}
								}
							}
						}
					}
				}
			}
		}
	}
	
	return num_rings;

}

/*
Description:
Given an adjacency matrix, a, find all unique closed 5 rings that are contained within it (i.e. a_ij = a_jk = a_kp = a_pq = a_qi = 1; a_ik = a_jp = a_pi = a_jq = a_kq = 0 - crossterms don't touch).

(Note that I only wrote this function when it became apparent that I was going to use it alot; therefore, some functions that could use this, don't, and instead have their own check - such as the check for whether a closed 5 ring surrounds a circle of intersection - I can go through later an dreplace these to use this function, if desired).

Output: The number of unique closed 5 rings.  The array containing all of the closed 5 rings is also updated.
*/
int find_5rings(int a[][MAXLEN2], int (*closed_5rings)[MAXLEN5], const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, num_rings;
	int test_ring[5];
	num_rings = 0;
	
	for (i = 0; i < max_col; ++i) {
		for (j = 0; j < max_col; ++j) {
			if (a[i][j] == 1) {
				for (k = 0; k < max_col; ++k) {
					if (k != i) {
						if ((a[j][k] == 1) && (a[i][k] == 0)) {
							for (p = 0; p < max_col; ++p) {
								if ((p != i) && (p != j)) {
									if ((a[k][p] == 1) && (a[p][i] == 0) && (a[p][j] == 0)) {
										for (q = 0; q < max_col; ++q) {
											if ((q != j) && (q != k)) {
												if ((a[p][q] == 1) && (a[q][i] == 1) && (a[q][j] == 0) && (a[q][k] == 0)) { //then we have a closed 5 ring, so check to see whether it is unique
													test_ring[0] = i;
													test_ring[1] = j;
													test_ring[2] = k;
													test_ring[3] = p;
													test_ring[4] = q;
													unique_ring_check(closed_5rings, test_ring, 5, &num_rings);
												}
											}
										}
									}
								}
							}
						}
					}
				}
			}
		}
	}
	
	return num_rings;

}

/*
Description:
Given an adjacency matrix, a, find all unique closed 6 rings that are contained within it (i.e. a_ij = a_jk = a_kp = a_pq = a_qz = a_zi = 1; a_ik = a_ip = a_iq = a_jp = a_jq = a_jz = a_kq = a_kz = a_pz = 0 - crossterms don't touch).

Output: The number of unique closed 6 rings.  The array containing all of the closed 5 rings is also updated.
*/
int find_6rings(int a[][MAXLEN2], int (*closed_6rings)[MAXLEN5], const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, z, num_rings;
	int test_ring[6];
	num_rings = 0;
	
	for (i = 0; i < max_col; ++i) {
		for (j = 0; j < max_col; ++j) {
			if (a[i][j] == 1) {
				for (k = 0; k < max_col; ++k) {
					if (k != i) {
						if ((a[j][k] == 1) && (a[i][k] == 0)) {
							for (p = 0; p < max_col; ++p) {
								if ((p != i) && (p != j)) {
									if ((a[k][p] == 1) && (a[p][i] == 0) && (a[p][j] == 0)) {
										for (q = 0; q < max_col; ++q) {
											if ((q != i) && (q != j) && (q != k)) {
												if ((a[p][q] == 1) && (a[q][i] == 0) && (a[q][j] == 0) && (a[q][k] == 0)) {
													for (z = 0; z < max_col; ++z) {
														if ((z != j) && (z != k) && (z != p)) {
															if ((a[q][z] == 1) && (a[z][i] == 1) && (a[z][j] == 0) && (a[z][k] == 0) && (a[z][p] == 0)) { //then we have a closed 6 ring, so check to see whether it is unique
																test_ring[0] = i;
																test_ring[1] = j;
																test_ring[2] = k;
																test_ring[3] = p;
																test_ring[4] = q;
																test_ring[5] = z;
																unique_ring_check(closed_6rings, test_ring, 6, &num_rings);
															}
														}
													}
												}
											}
										}
									}
								}
							}
						}
					}
				}
			}
		}
	}
	
	return num_rings;

}

/*
Description:
Accepts an array of points, and checks whether m out of the n total points are touched by a given particle (where the particle is required to touch m of these points).  If m < n, it also documents which points were not touched, by updating an array containing these points.  'row' indicates the row of the adjacency matrix, and thus the particle being tested.  (12 is the maximum of the points array, as a particle can only touch a maximum of 12 points).

Output: 1 if the number of required points are touched, 0 if they are not.  The array containing which of these points the particle does not touch is also updated.
*/
int check_connections(int a[][MAXLEN2], int points[12], int *not_touched_points, int *not_touched_idx, const int num_points, const int num_need_touch, const int row) {

	/*Declare and initialize local variables*/
	int i, count, num_not_touched;
	count = 0;
	num_not_touched = 0;
	
	for (i = 0; i < num_points; ++i) {
		if (a[row][points[i]] == 1)
			++count;
		else {
			not_touched_points[num_not_touched] = points[i]; //record the points the particle doesn't touch
			not_touched_idx[num_not_touched] = i;
			++num_not_touched;
		}
	}
	if (count == num_need_touch)
		return 1;
	else
		return 0;

}

/*
Description:
Check all dihedral angle combinations for all spherical triangles that involve line segment rij, and see if any imply that rij has been calculated incorrectly.

Note that i and j are entered in 1 based.

Output: Whether all possible spherical triangles are consistent or not - if at least 1 is found not consistent, return 0, if all are found consistent, return 1.
*/
int consistency_check(double test_sol[], const int max_col, const int i, const int j, const double dist) {

	/*Declare local variables*/
	int idx_ij, idx_ik, idx_ip, idx_jk, idx_jp, idx_pk, k, p, consistent;
	consistent = 1;
	
	/*For each relative distance rij, there will be 2 other points that will specify a spherical triange with that line segment.  Thus, go through all these possibilities, and check for an inconsistency in any spherical triangle with enough information to calculate the dihedral angles, etc.*/
	for (k = 1; k <= max_col; ++k) {
		for (p = k+1; p <= max_col; ++p) {
			if (k != i && k != j && p != i && p != j) {
				if (i < k)
					idx_ik = row_col_2_idx(i,k,max_col);
				else
					idx_ik = row_col_2_idx(k,i,max_col);
				if (i < p)
					idx_ip = row_col_2_idx(i,p,max_col);
				else
					idx_ip = row_col_2_idx(p,i,max_col);
				if (j < k)
					idx_jk = row_col_2_idx(j,k,max_col);
				else
					idx_jk = row_col_2_idx(k,j,max_col);
				if (j < p)
					idx_jp = row_col_2_idx(j,p,max_col);
				else
					idx_jp = row_col_2_idx(p,j,max_col);
				idx_pk = row_col_2_idx(k,p,max_col); //k is always < p from this for loop
				if (test_sol[idx_ik] != DUMMY_VAR && test_sol[idx_ip] != DUMMY_VAR && test_sol[idx_jk] != DUMMY_VAR && test_sol[idx_jp] != DUMMY_VAR && test_sol[idx_pk] != DUMMY_VAR) { //only check the spherical trangle for 																																														consistency if all of its relative 																																														distances are known
					consistent = spherical_triangle_check(dist, test_sol[idx_ik], test_sol[idx_ip], test_sol[idx_jk], test_sol[idx_jp], test_sol[idx_pk]);
					//printf("Spherical Triangle = %d,%d,%d,%d\nConsistent = %d\n",i,j,k,p,consistent);
					//printf("The distances for the spherical triangle are r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f\n",i,j,dist,i,k,test_sol[idx_ik],i,p,test_sol[idx_ip],j,k,test_sol[idx_jk],j,p,test_sol[idx_jp],p,k,test_sol[idx_pk]);
				}
				if (consistent == 0)
					break;
			}
		}
		if (consistent == 0)
			break;
	}
	
	return consistent;

}

/*
Description:
Check the spherical triangle that involves an unknown distance which was solved for consistency by seeing whether the spherical law of cosines is satisfied for all of its dihedral angles.  If it can not be satisfied, I believe nan should arise.

Output: 1 if the spherical triangle is consistent, 0 if not.
*/
double spherical_triangle_check(const double r_ij, const double r_ik, const double r_ip, const double r_jk, const double r_jp, const double r_pk) {

	/*Declare and initialize local variables*/
	int consistent, round_num;
	double a, b, c, A, B, C, arg_A, arg_B, arg_C, num_A, num_B, num_C, PI;
	consistent = 1;
	round_num = 5; //originally set to 5 - this is the precision to which the dihedral angle = 0 or pi is checked.  Changing to 3 didn't make a difference.
	
	//Pi to 1024 decimal places:
	PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327887;
	
	/*Calculate all of the known angles from the relative distances*/
	/*a = law_of_cosines_num(r_ip,r_jp,r_ij,2);
	b = law_of_cosines_num(r_ik,r_pk,r_ip,2);
	c = law_of_cosines_num(r_ik,r_jk,r_ij,2);*/ //had it like this until 6/02/09 -- this is incorrect, I don't know how it worked as well as it did with this...
	a = law_of_cosines_num(r_ip,r_jp,r_ij,3);
	b = law_of_cosines_num(r_jk,r_pk,r_jp,1);
	c = law_of_cosines_num(r_ik,r_pk,r_ip,1);
	//printf("r_ij = %f, r_ik = %f, r_ip = %f, r_jk = %f, r_jp = %f, r_pk = %f\n",r_ij,r_ik,r_ip,r_jk,r_jp,r_pk);
	//printf("a = %f, b = %f, c = %f\n",a,b,c);
	
	/*if any of the angles = pi (180 degrees) or 0, then we have a line, not a triangle, and the dihedral angle is not well definted, and thus this check can not be performed - 
	and nothing can be invalidated by this check in this case, so return consistent = 1 by default in this case.  Thus, only calculate the dihedral angles if they are well defined.*/
	if (isnan(a) == 1 || isnan(b) == 1 || isnan(c) == 1) { //first check that all of the angles are well defined, if not return inconsistent, if well-defined proceed with check
		consistent = 0;
	}
	//may as well change this to >0 instead of >=1, but will save that for later.
	else if ((fabs(round(a*pow(10,round_num)) - round(PI*pow(10,round_num))) >= 1) && (fabs(round(b*pow(10,round_num)) - round(PI*pow(10,round_num))) >= 1) && (fabs(round(c*pow(10,round_num)) - round(PI*pow(10,round_num))) >= 1) && (round(a*pow(10,round_num)) >= 1) && (round(b*pow(10,round_num)) >= 1) && (round(c*pow(10,round_num)) >= 1)) { //check that none of the angles = pi or 0 to round_num decimal places
		/*Calculate the dihedral angles associated with this spherical triangle*/
		num_A = cos(a) - cos(b)*cos(c);
		arg_A = num_A/(sin(c)*sin(b));
		//printf("arg_A = %f\n",arg_A);
		arg_A = check_round(arg_A,num_A);
		A = acos(arg_A);
		num_B = cos(b) - cos(a)*cos(c);
		arg_B = num_B/(sin(c)*sin(a));
		//printf("\nnumB = %f, cos(b) = %f, cos(a) = %f, cos(c) = %f, sin(a) = %f, sin(c) = %f, denomB = %f\n",num_B,cos(b),cos(a),cos(c),sin(a),sin(c),sin(a)*sin(c));
		//printf("\narg = %f, arg = %f\n",(cos(b) - cos(a)*cos(c))/(sin(a)*sin(c)),num_b/(sin(a)*sin(c)));
		//printf("arg_B = %f\n",arg_B);
		arg_B = check_round(arg_B,num_B);
		B = acos(arg_B);
		num_C = cos(c) - cos(a)*cos(b);
		arg_C = num_C/(sin(a)*sin(b));
		//printf("arg_C = %f\n",arg_C);
		arg_C = check_round(arg_C,num_C);
		C = acos(arg_C);
		//printf("A = %f, B = %f, C = %f\n",A,B,C);
		//printf("arg_A = %f, arg_B = %f, arg_C = %f\n",arg_A,arg_B,arg_C);
		//printf("abs(arg_A) = %f, abs(arg_B) = %f, abs(arg_C) = %f\nconsistent = %d\n",fabs(arg_A),fabs(arg_B),fabs(arg_C),consistent);
		
		//if (fabs(arg_A) > 1 || isnan(arg_A) == 1 || fabs(arg_B) > 1 || isnan(arg_B) == 1 || fabs(arg_C) > 1 || isnan(arg_C) == 1) {
		if (isnan(A) == 1 || isnan(B) == 1 || isnan(C) == 1) {
			consistent = 0;
			//printf("\nIn Sphere Triang Check:\nabs(arg_A) = %f, abs(arg_B) = %f, abs(arg_C) = %f\nconsistent = %d\n",fabs(arg_A),fabs(arg_B),fabs(arg_C),consistent);
		}
	}
		
	return consistent;

}

/*
Description:
Recalculate the dihedral angles when have zeros cancelling (as example in beginning of Self-Assembly notebook 13).  This function performs the zero cancelling.
Here we will calculate the zero cancelling in the dihedral angle as follows:
The formula for the dihedral angle, A, is given by
cos(a) - cos(b)*cos(c) = sin(b)*sin(c)*cos(A)
This is called when either sin(b) or sin(c) = 0 and cos(a) - cos(b)*cos(c) = 0.  In this case, the non-zeros sin function, let this be sin(b), for example, will multiply out the cosines, 
and the zeros will cancel:
i.e. sin(b)*cos(a) - sin(b)*cos(b)*cos(c) = sin(c)*cos(A), where we simplify the LHS by cancelling out complementary numerators and denominators btwn sin(b) and the cos terms.  Then, 
the simplified cos(a) - cos(b)*cos(c) cancels out with sin(c), and we have A = arccos(simplified version of sin(b)).

If sin(b) and sin(c) = 0, and cos(a) - cos(b)*cos(c) = 0 then we will simly try A = arccos(1) -- although I'm not sure I should be doing this...but for now we will try this...

What we will be using for these cos and sin terms are as follows:
a = arccos(f(dists)), and the same for b and c.  The dihedral angle formula thus has cos(arccos(a)) and sin(arccos(a)), and so on.  cos(arccos(x)) = x, and sin(arccos(x)) = sqrt(1-x^2), so
we will directly evaluate those terms as opposed to the cos and sin terms.  I.e. we will evaluate
f(dists) - f(dists)*f(dists) = sqrt(1-f(dists)^2)*sqrt(1 - f(dists)^2)*cos(A)
This is
(r_ab^2 + r_ac^2 - r_a^2)/(2*r_ab*r_ac) - (r_bc^2 + r_ab^2 - r_b^2)/(2*r_bc*r_ab) * (r_ac^2 + r_bc^2 - r_c^2)/(2*r_ac*r_bc) = sqrt(1 - ((r_bc^2 + r_ab^2 - r_b^2)/(2*r_bc*r_ab))^2)*sqrt(1 - (r_ac^2 + r_bc^2 - r_c^2)^2/(4*r_ac^2*r_bc^2))*cos(A)
And I want to evaluate numerator and denominator cancelling separately...
I'll combine the numerator into one term, keep the denominator as separate terms...and do the cancelling that way...

Inputs: The distances of the tetrahedron for which the dihedral angle, A, is calculated.
		r_bc: the distance 'OA'(from my general spherical triangle) -- this is the side distance used in the calculation of angles b and c -- i.e. this distance separates angle b from angle c.
		r_ab: the distance 'OC' -- this is the side distance used in the calculation of angles a and b -- i.e. this distance separates angle a from angle b.
		r_ac: the distance 'OB' -- this is the side distance used in the calculation of angles a and c -- i.e. this distance separates angle a from angle c.
		r_b: the distance 'AC' -- this is the distance opposite angle b, it is only used in the calculation of angle b.
		r_a: the distance 'BC' -- this is the distance opposite angle a, it is only used in the calculation of angle a.
		r_c: the distance 'AB' -- this is the distance opposite angle c, it is only used in the calculation of angle c.

Output: A, the dihedral angle.
*/
void recalc_dihed_ang_0cancelling(double A[], int *numsA, const double r_bc, const double r_ab, const double r_ac, const double r_b, const double r_a, const double r_c) {

	/*Declare and Initialize Local Variables*/
	int prec, sin_b_0, sin_c_0, cos_a_0, cos_b_0, cos_c_0, mult, num_greater, simp;
	double r_ab_sq, r_bc_sq, r_ac_sq, num_cos_a, num_cos_b, num_cos_c, num_sin_b, num_sin_c, num_cos_b_sq, num_cos_c_sq, denom_a, denom_b, denom_c, sin_b, sin_c, cos_a, cos_b, cos_c, comm_val, num_arg, denom_arg, arg_A;
	double nums_cosa[1], simp_nums_cosa[1], denoms_a[3], simp_denoms_a[3], nums_cosb[1], simp_nums_cosb[1], denoms_b[3], simp_denoms_b[3], nums_cosc[1], simp_nums_cosc[1], denoms_c[3], simp_denoms_c[3], nums_sinb[1], simp_nums_sinb[1], simp_denoms_sinb[3], nums_sinc[1], simp_nums_sinc[1], simp_denoms_sinc[3], num_arg1[3], denom_arg1[3], num_arg2[1], denom_arg2[1];
	prec = 5; //originally set to 5 - this is the precision to which each angle = 0 is checked.  Changing to 3 or 7 makes no difference, changing to 10 makes 1 packing be missed out.
	sin_b_0 = 0;
	sin_c_0 = 0;
	cos_a_0 = 0;
	cos_b_0 = 0;
	cos_c_0 = 0;
	*numsA = 2;
	simp = 1;
	
	r_ab_sq = pow(r_ab,2);
	r_bc_sq = pow(r_bc,2);
	r_ac_sq = pow(r_ac,2);
	num_cos_a = r_ab_sq + r_ac_sq - pow(r_a,2);
	num_cos_b = r_bc_sq + r_ab_sq - pow(r_b,2);
	num_cos_c = r_ac_sq + r_bc_sq - pow(r_c,2);
	num_cos_b_sq = pow(num_cos_b,2);
	num_cos_c_sq = pow(num_cos_c,2);
	denom_a = 2*r_ab*r_ac;
	denom_b = 2*r_bc*r_ab;
	denom_c = 2*r_ac*r_bc;
	num_sin_b = sqrt(4*r_bc_sq*r_ab_sq - num_cos_b_sq);
	num_sin_c = sqrt(4*r_ac_sq*r_bc_sq - num_cos_c_sq);
	sin_b = num_sin_b/denom_b;
	sin_c = num_sin_c/denom_c;
	cos_a = num_cos_a/denom_a;
	cos_b = num_cos_b/denom_b;
	cos_c = num_cos_c/denom_c;
	
	
	if (sin_b*pow(10,prec) <= 1)
		sin_b_0 = 1;
	if (sin_c*pow(10,prec) <= 1)
		sin_c_0 = 1;
	if (cos_a*pow(10,prec) <= 1)
		cos_a_0 = 1;
	if (cos_b*pow(10,prec) <= 1)
		cos_b_0 = 1;
	if (cos_c*pow(10,prec) <= 1)
		cos_c_0 = 1;
	
	if ((sin_b_0 == 1) && (sin_c_0 == 1)) { //if both sinb and sinc = 0
		A[0] = acos(1);
		A[1] = acos(-1); //the negative solution
	}
	else {
		/*Simplify expressions for the cosines of each of the angles, by simplifying the numerator individually with each term of the denominator.*/
		denoms_a[0] = 2.0;
		denoms_a[1] = r_ab;
		denoms_a[2] = r_ac;
		nums_cosa[0] = num_cos_a;
		simp_nums_denoms_from_intsimp(simp_nums_cosa, simp_denoms_a, nums_cosa, denoms_a, 1, 3);
		//printf("\nSimplifying A:\nsimp_anum = %f, denom: simpa1 = %f, simpa2 = %f, simpa3 = %f\nanum = %f, denom: a1 = %f, a2 = %f, a3 = %f.\n",simp_nums_cosa[0],simp_denoms_a[0],simp_denoms_a[1],simp_denoms_a[2],nums_cosa[0],denoms_a[0],denoms_a[1],denoms_a[2]);
		denoms_b[0] = 2.0;
		denoms_b[1] = r_bc;
		denoms_b[2] = r_ab;
		nums_cosb[0] = num_cos_b;
		simp_nums_denoms_from_intsimp(simp_nums_cosb, simp_denoms_b, nums_cosb, denoms_b, 1, 3);
		denoms_c[0] = 2.0;
		denoms_c[1] = r_ac;
		denoms_c[2] = r_bc;
		nums_cosc[0] = num_cos_a;
		simp_nums_denoms_from_intsimp(simp_nums_cosc, simp_denoms_c, nums_cosc, denoms_c, 1, 3);
		nums_sinb[0] = num_sin_b;
		simp_nums_denoms_from_intsimp(simp_nums_sinb, simp_denoms_sinb, nums_sinb, denoms_b, 1, 3);
		nums_sinc[0] = num_sin_c;
		simp_nums_denoms_from_intsimp(simp_nums_sinc, simp_denoms_sinc, nums_sinc, denoms_c, 1, 3);
		if (sin_b_0 == 1) { //if only sin b = 0 (i.e. sin(c) is non-zero)
			simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinc, simp_denoms_a, 3, 3);
			simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosa, simp_nums_sinc, 1, 1);
			if (denom_arg2[0] == simp_nums_sinc[0])
				simp = 0;
		}
		else if (sin_c_0 == 1) { //if only sin c = 0 (i.e. sin(b) is non-zero)
			simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinb, simp_denoms_a, 3, 3);
			simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosa, simp_nums_sinb, 1, 1);
			if (denom_arg2[0] == simp_nums_sinb[0])
				simp = 0;
		}
		if ((num_arg2[0] == simp_nums_cosa[0]) && (simp == 0)) { //if no simplification was made for (simplified cos a num) / (simplified sin c num)
			num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2]; //numerator is simp denom sin c or sin b
		}
		else { //if a simplification was made
			num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2]*num_arg2[0]; //numerator is simp denom sin c (or b) * simp num cos a
		}
		denom_arg = denom_arg1[0]*denom_arg1[1]*denom_arg1[2]*denom_arg2[0];
		arg_A = num_arg/denom_arg;
		arg_A = check_round(arg_A,num_arg);
		A[0] = acos(arg_A);
		A[1] = acos(-arg_A); //the negative zero cancelling solution
		//printf("\nIn recalc of dihedral angle:\nNum_A = %f, denom_A = %f, arg_A = %f, A = %f\nsimp_nums_cosa = %f, simp_nums_sinb = %f, simp_denoms_sinb1 = %f, simp_denoms_sinb2 = %f, simp_denoms_sinb3 = %f, simp_denoms_a1 = %f, simp_denoms_a2 = %f, simp_denoms_a3 = %f\n",num_arg,denom_arg,arg_A,A[0],simp_nums_cosa[0],simp_nums_sinb[0],simp_denoms_sinb[0],simp_denoms_sinb[1],simp_denoms_sinb[2],simp_denoms_a[0],simp_denoms_a[1],simp_denoms_a[2]);
		//printf("\nsimplified: num_arg1[0] = %f, num_arg1[1] = %f, num_arg1[2] = %f, num_arg2[0] = %f, denom_arg1[0] = %f, denom_arg1[1] = %f, denom_arg1[2] = %f, denom_arg2[0] = %f.\n",num_arg1[0], num_arg1[1], num_arg1[2], num_arg2[0], denom_arg1[0], denom_arg1[1], denom_arg1[2], denom_arg2[0]);
		if (cos_a_0 == 1) { //do the alternate zero cancelling solution also if cos(a) = 0
			*numsA = 4;
			if ((cos_b_0 == 1) && (cos_c_0 == 1)) {
				A[2] = acos(1);
				A[3] = acos(-1);
			}
			else { //if either cos(b) or cos(c) = 0
				if ((cos_b_0 == 1) && (cos_c_0 == 0)) { //if cos(c) is non-zero
					if (sin_b_0 == 1) { //if only sin b = 0 (i.e. sin(c) is non-zero)
						simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinc, simp_denoms_c, 3, 3);
						simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosc, simp_nums_sinc, 1, 1);
						if (denom_arg2[0] == simp_nums_sinc[0])
							simp = 0;
					}
					else if (sin_c_0 == 1) { //if only sin c = 0 (i.e. sin(b) is non-zero)
						simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinb, simp_denoms_c, 3, 3);
						simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosc, simp_nums_sinb, 1, 1);
						if (denom_arg2[0] == simp_nums_sinb[0])
							simp = 0;
					}
					if ((num_arg2[0] == simp_nums_cosc[0]) && (simp == 0)) { //if no simplification was made for (simplified cos a num) / (simplified sin c num)
						num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2];
					}
					else { //if a simplification was made
						num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2]*num_arg2[0];
					}
				}
				else if ((cos_b_0 == 0) && (cos_c_0 == 1)) { //if cos(b) is non-zero
					if (sin_b_0 == 1) { //if only sin b = 0 (i.e. sin(c) is non-zero)
						simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinc, simp_denoms_b, 3, 3);
						simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosb, simp_nums_sinc, 1, 1);
						if (denom_arg2[0] == simp_nums_sinc[0])
							simp = 0;
					}
					else if (sin_c_0 == 1) { //if only sin c = 0 (i.e. sin(b) is non-zero)
						simp_nums_denoms_from_intsimp(num_arg1, denom_arg1, simp_denoms_sinb, simp_denoms_b, 3, 3);
						simp_nums_denoms_from_intsimp(num_arg2, denom_arg2, simp_nums_cosb, simp_nums_sinb, 1, 1);
						if (denom_arg2[0] == simp_nums_sinb[0])
							simp = 0;
					}
					if ((num_arg2[0] == simp_nums_cosb[0]) && (simp == 0)) { //if no simplification was made for (simplified cos a num) / (simplified sin c num)
						num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2];
					}
					else { //if a simplification was made
						num_arg = num_arg1[0]*num_arg1[1]*num_arg1[2]*num_arg2[0];
					}
				}
				denom_arg = denom_arg1[0]*denom_arg1[1]*denom_arg1[2]*denom_arg2[0];
				arg_A = num_arg/denom_arg;
				arg_A = check_round(arg_A,num_arg);
				A[2] = acos(arg_A);
				A[3] = acos(-arg_A); //the negative zero cancelling solution
			}
		}
	}

}

/*
Description: 
Calculates the relative distance between two points, when those 2 points can be decomposed into one polygon added to another (where all the relative distances amongst those 2 polygons are known).  This function simply uses the standard math.h library to give the analytically calculated answer numerically.  Note that i, j, k, p, q are entered in 1 based.

col is the potential solution column, row is the unknown relative distance index.

Also note that without loss of generality, we have set pk = 1 for this general rule.

Output: The number of consistent solutions, given a partial solution vector as input.  The relative distance between the 2 points is also outputted numerically (as a double) by having the num_sol array updated.
*/
int calc_rel_dist_num(double (*mult_num_sols_write)[MAXLEN2], double (*mult_num_sols_old)[MAXLEN2], const int col, const int row, int *num_new_sols, const int max_col, const int num_int_circs, const int i, const int j, const int k, const int p, const int q) {

	/*Declare local variables*/
	int idx_ik, idx_ip, idx_iq, idx_jk, idx_jp, idx_jq, idx_pk, idx_pq, idx_kq, orientation, num_cons_sols, z, s, l, recalc_A2, recalc_A3, nums_A2, nums_A3, consistent, round_num1, round_num2, round_num3, bla1, bla2;
	double a1, b1, c1, a2, b2, c2, a3, b3, c3, A1, A2, A3, arg_A2, arg_A3, arg_a1, num_A2, num_A3, dist_old, dist_old2, dist_old3, dist, PI;
	double test_sol[num_int_circs], A2_array[4], A3_array[4];
	dist = DUMMY_VAR;
	dist_old = DUMMY_VAR;
	consistent = 0;
	num_cons_sols = 0;
	round_num1 = 5; //originally set to 5 - this is the precision to which we determine whether we have a new distance and thus whether we should check it.  Changing this to 3 caused a packing to be missed out, and changing it to 10 made no difference.
	round_num2 = 5; //originally set to 5 - this is the precision to which we check whether the numerator and denominator in zero cancelling are both zero.  Changing this to 3, 7, 10, and 16 made no difference in the output.
	round_num3 = 16; //set at 5, 8, 16 with the anlge (b1, a2, etc.) check, 1 packing is missed out at n=9.  Set at 5 for dihedral angle not = 0, Pi, 2Pi check it only registers 3 packings at n = 9.
	
	
	/*printf("\nIn Calc Rel Dist:\n");
	printf("col = %d\n",col);
	for (z = 0; z < num_int_circs; ++z) {
		bla1 = idx2row(z,max_col);
		bla2 = idx2col(z,bla1,max_col);
		printf("mn_sol_in[%d][%d] = %f\n",bla1,bla2,mult_num_sols_old[col][z]);
	}
	
	printf("\nLocal Solution:\n");
	printf("\ni = %d, j = %d, k = %d, p = %d, q = %d\n",i,j,k,p,q);*/
	
	/*Get the indices associated with all the relative distances.*/
	if (i < k)
		idx_ik = row_col_2_idx(i,k,max_col);
	else
		idx_ik = row_col_2_idx(k,i,max_col);
	if (i < p)
		idx_ip = row_col_2_idx(i,p,max_col);
	else
		idx_ip = row_col_2_idx(p,i,max_col);
	if (i < q)
		idx_iq = row_col_2_idx(i,q,max_col);
	else
		idx_iq = row_col_2_idx(q,i,max_col);
	if (j < k)
		idx_jk = row_col_2_idx(j,k,max_col);
	else
		idx_jk = row_col_2_idx(k,j,max_col);
	if (j < p)
		idx_jp = row_col_2_idx(j,p,max_col);
	else
		idx_jp = row_col_2_idx(p,j,max_col);
	if (j < q)
		idx_jq = row_col_2_idx(j,q,max_col);
	else
		idx_jq = row_col_2_idx(q,j,max_col);
	if (p < k)
		idx_pk = row_col_2_idx(p,k,max_col);
	else
		idx_pk = row_col_2_idx(k,p,max_col);
	if (p < q)
		idx_pq = row_col_2_idx(p,q,max_col);
	else
		idx_pq = row_col_2_idx(q,p,max_col);
	if (k < q)
		idx_kq = row_col_2_idx(k,q,max_col);
	else
		idx_kq = row_col_2_idx(q,k,max_col);
		
	//printf("\nidx_ik = %d, idx_ip = %d, idx_iq = %d, idx_jk = %d, idx_jp = %d, idx_jq = %d, idx_pk = %d, idx_pq = %d, idx_kq = %d\n",idx_ik,idx_ip,idx_iq,idx_jk,idx_jp,idx_jq,idx_pk,idx_pq,idx_kq);
	//printf("\nnum_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f, num_sol[%d] = %f\n",idx_ik,num_sol[idx_ik],idx_ip,num_sol[idx_ip],idx_iq,num_sol[idx_iq],idx_jk,num_sol[idx_jk],idx_jp,num_sol[idx_jp],idx_jq,num_sol[idx_jq],idx_pk,num_sol[idx_pk],idx_pq,num_sol[idx_pq],idx_kq,num_sol[idx_kq]);
	
	/*Calculate all of the known angles from the relative distances*/
	b1 = law_of_cosines_num(mult_num_sols_old[col][idx_jp],mult_num_sols_old[col][idx_jk],mult_num_sols_old[col][idx_pk],2);
	c1 = law_of_cosines_num(mult_num_sols_old[col][idx_ip],mult_num_sols_old[col][idx_ik],mult_num_sols_old[col][idx_pk],2);
	a2 = law_of_cosines_num(mult_num_sols_old[col][idx_ip],mult_num_sols_old[col][idx_iq],mult_num_sols_old[col][idx_pq],2);
	b2 = law_of_cosines_num(mult_num_sols_old[col][idx_kq],mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_pq],1);
	c2 = law_of_cosines_num(mult_num_sols_old[col][idx_ip],mult_num_sols_old[col][idx_ik],mult_num_sols_old[col][idx_pk],2);
	a3 = law_of_cosines_num(mult_num_sols_old[col][idx_jp],mult_num_sols_old[col][idx_jq],mult_num_sols_old[col][idx_pq],2);
	b3 = law_of_cosines_num(mult_num_sols_old[col][idx_jk],mult_num_sols_old[col][idx_jp],mult_num_sols_old[col][idx_pk],1);
	c3 = law_of_cosines_num(mult_num_sols_old[col][idx_kq],mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_pq],1);
	
	//printf("\nb1 = %f, c1 = %f, a2 = %f, b2 = %f, c2 = %f, a3 = %f, b3 = %f, c3 = %f\n",b1,c1,a2,b2,c2,a3,b3,c3); //uncomment this when want to check
	//printf("r_jp = %f, r_jq = %f, r_pq = %f\n",num_sol[idx_jp],num_sol[idx_jq],num_sol[idx_pq]);
	//printf("r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f\n",j,p,mult_num_sols_old[col][idx_jp],j,q,mult_num_sols_old[col][idx_jq],p,q,mult_num_sols_old[col][idx_pq]);
	//printf("r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f, r_%d_%d = %f\n",j,p,mult_num_sols_old[col][idx_jp],j,q,mult_num_sols_old[col][idx_jq],p,q,mult_num_sols_old[col][idx_pq],i,p,mult_num_sols_old[col][idx_ip],i,q,mult_num_sols_old[col][idx_iq],i,k,mult_num_sols_old[col][idx_ik],j,p,mult_num_sols_old[col][idx_jp],j,q,mult_num_sols_old[col][idx_jq],j,k,mult_num_sols_old[col][idx_jk]); //uncomment this when want to check
	
	/*Calculate the relative distance*/
	num_A2 = cos(a2) - cos(b2)*cos(c2);
	arg_A2 = num_A2/(sin(c2)*sin(b2));
	num_A3 = cos(a3) - cos(b3)*cos(c3);
	arg_A3 = num_A3/(sin(c3)*sin(b3));
	arg_A2 = check_round(arg_A2,num_A2);
	arg_A3 = check_round(arg_A3,num_A3);
	A2 = acos(arg_A2); //note we don't need to put pow(num_sol[idx_pk],2) in any of the expressions for this general rule, as the spherical triangles can always be normalized s.t. this = 1
	A3 = acos(arg_A3);
	//printf("A2 = acos(%f), A3 = acos(%f)\n",arg_A2,arg_A3);
	//printf("A2 = %f, A3 = %f\n",A2,A3);
	
	
	/*Try the summing dihedral angles option*/
	A1 = A2 + A3; //first try summing the angles
	arg_a1 = sin(c1)*sin(b1)*cos(A1) + cos(b1)*cos(c1);
	arg_a1 = check_round(arg_a1,1); //just put 1 here as the numerator as a1 isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 in the numerator will ensure it doesn't round to 1
	a1 = acos(arg_a1);
	if (isnan(a1) == 0) {
		dist = sqrt(pow(mult_num_sols_old[col][idx_ip],2) + pow(mult_num_sols_old[col][idx_jp],2) - 2*mult_num_sols_old[col][idx_ip]*mult_num_sols_old[col][idx_jp]*cos(a1));
		dist = check_round(dist, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
		dist_old = dist;
		//printf("Old dist = %f\n",dist_old);
		//printf("A1 = A2 + A3 gives dist = %f\n",dist);
		if (dist >= 1) {
			for (z = 0; z < num_int_circs; ++z)
				test_sol[z] = mult_num_sols_old[col][z];
			test_sol[row] = dist;
			consistent = consistency_check(test_sol, max_col, i, j, dist);
			if (consistent == 1) {
				for (z = 0; z < num_int_circs; ++z)
					mult_num_sols_write[*num_new_sols][z] = test_sol[z];
				//printf("\nI was called 1\n");
				//printf("\nSolution Recorded.\n");
				++*num_new_sols;
				++num_cons_sols;
				//printf("\nSumming A2 and A3 found a consistent solution:\nA1 = %f, a1 = acos(%f), a1 = %f, Implies dist = %f\n",A1,arg_a1,a1,dist);
			}
			//else
				//printf("\nLocal Distance is Inconsistent!\n");
		}
	}
	/*Try the subtracting dihedral angles option*/
	if (A2 < A3) {
		A1 = A3 - A2;
		//printf("A1 = A3 - A2\n");
	}
	else {
		A1 = A2 - A3;
		//printf("A1 = A2 - A3\n");
	}
	arg_a1 = sin(c1)*sin(b1)*cos(A1) + cos(b1)*cos(c1);
	arg_a1 = check_round(arg_a1, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
	a1 = acos(arg_a1);
	dist = sqrt(pow(mult_num_sols_old[col][idx_ip],2) + pow(mult_num_sols_old[col][idx_jp],2) - 2*mult_num_sols_old[col][idx_ip]*mult_num_sols_old[col][idx_jp]*cos(a1));
	dist = check_round(dist, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
	dist_old2 = dist;
	//printf("New dist = %f\n",dist);
	//printf("New dist = %f, old dist = %f\n",dist,dist_old);
	//printf("Subtracting A2 and A3 gives dist = %f\n",dist);
	//if ((dist >= 1) && (dist != dist_old)) { //only record or test unique distances
	if ((dist >= 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old*pow(10,round_num1))) > 1)) { //only record or test unique distances - only check that they're equivalent out to round_num1 decimal places...
		for (z = 0; z < num_int_circs; ++z)
			test_sol[z] = mult_num_sols_old[col][z];
		test_sol[row] = dist;
		consistent = consistency_check(test_sol, max_col, i, j, dist);
		if (consistent == 1) {
			for (z = 0; z < num_int_circs; ++z)
				mult_num_sols_write[*num_new_sols][z] = test_sol[z];
			//printf("\nI was called 2\n");
			//printf("\nSolution Recorded.\n");
			++*num_new_sols;
			++num_cons_sols;
			//printf("\nSubtracting A2 and A3 found a consistent solution:\nA1 = %f, a1 = acos(%f), a1 = %f, Implies dist = %f\n",A1,arg_a1,a1,dist);
		}
		//else
			//printf("\nLocal Distance is Inconsistent!\n");
	} //end of calculating distances directly	
	
	recalc_A2 = 0; //initialize
	recalc_A3 = 0;
	A2_array[0] = A2;
	A3_array[0] = A3;
	nums_A2 = 1;
	nums_A3 = 1;
	/*Check if any zeros cancel, and recalc dihedral angles and distance if any zeros cancel*/
	if ((fabs(round(num_A2*pow(10,round_num2))) < 1) && ((fabs(round(sin(c2)*pow(10,round_num2))) < 1) || (fabs(round(sin(b2)*pow(10,round_num2))) < 1))) {
		recalc_A2 = 1;
		//A2 = recalc_dihed_ang_0cancelling(sol_num[idx_pk],sol_num[idx_pq],sol_num[idx_ip],sol_num[idx_kq],sol_num[idx_iq],sol_num[idx_ik]);
		//A2 = recalc_dihed_ang_0cancelling(mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_pq],mult_num_sols_old[col][idx_ip],mult_num_sols_old[col][idx_kq],mult_num_sols_old[col][idx_iq],mult_num_sols_old[col][idx_ik]);
		recalc_dihed_ang_0cancelling(A2_array, &nums_A2, mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_pq],mult_num_sols_old[col][idx_ip],mult_num_sols_old[col][idx_kq],mult_num_sols_old[col][idx_iq],mult_num_sols_old[col][idx_ik]);
	}
	if ((fabs(round(num_A3*pow(10,round_num2))) < 1) && ((fabs(round(sin(c3)*pow(10,round_num2))) < 1) || (fabs(round(sin(b3)*pow(10,round_num2))) < 1))) {
		recalc_A3 = 1;
		//A3 = recalc_dihed_ang_0cancelling(sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_pq],sol_num[idx_jk],sol_num[idx_jq],sol_num[idx_kq]);
		//A3 = recalc_dihed_ang_0cancelling(mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_jp],mult_num_sols_old[col][idx_pq],mult_num_sols_old[col][idx_jk],mult_num_sols_old[col][idx_jq],mult_num_sols_old[col][idx_kq]);
		recalc_dihed_ang_0cancelling(A3_array, &nums_A3, mult_num_sols_old[col][idx_pk],mult_num_sols_old[col][idx_jp],mult_num_sols_old[col][idx_pq],mult_num_sols_old[col][idx_jk],mult_num_sols_old[col][idx_jq],mult_num_sols_old[col][idx_kq]);
	}
	if ((recalc_A2 == 1) || (recalc_A3 == 1)) { //if any zero canceling occurs, then recalculate the distances using the recalculated dihedral angles
		//printf("\nZero Cancelling Solution:\n");
		for (s = 0; s < nums_A2; ++s) {
			A2 = A2_array[s];
			for (l = 0; l < nums_A3; ++l) {
				A3 = A3_array[l];
				/*Try the summing dihedral angles option*/
				 A1 = A2 + A3; //first try summing the angles
				 arg_a1 = sin(c1)*sin(b1)*cos(A1) + cos(b1)*cos(c1);
				 arg_a1 = check_round(arg_a1,1); //just put 1 here as the numerator as a1 isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 in the numerator will ensure it doesn't round to 1
				 a1 = acos(arg_a1);
				 if (isnan(a1) == 0) {
					 dist = sqrt(pow(mult_num_sols_old[col][idx_ip],2) + pow(mult_num_sols_old[col][idx_jp],2) - 2*mult_num_sols_old[col][idx_ip]*mult_num_sols_old[col][idx_jp]*cos(a1));
					 dist = check_round(dist, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
					 dist_old3 = dist;
					 //printf("Old dist = %f\n",dist_old);
					 //printf("A1 = A2 + A3 gives dist = %f\n",dist);
					 if ((dist >= 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old*pow(10,round_num1))) > 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old2*pow(10,round_num1))) > 1)) { //only record or test unique distances - only check that they're equivalent out to round_num1 decimal places...
						 for (z = 0; z < num_int_circs; ++z)
							 test_sol[z] = mult_num_sols_old[col][z];
						 test_sol[row] = dist;
						 consistent = consistency_check(test_sol, max_col, i, j, dist);
						 if (consistent == 1) {
							 for (z = 0; z < num_int_circs; ++z)
								 mult_num_sols_write[*num_new_sols][z] = test_sol[z];
							 //printf("\nI was called 1\n");
							 //printf("\nSolution Recorded.\n");
							 ++*num_new_sols;
							 ++num_cons_sols;
							 //printf("\nFrom Zero Canceling: Summing A2 and A3 found a consistent solution:\nA1 = %f, a1 = acos(%f), a1 = %f, Implies dist = %f\n",A1,arg_a1,a1,dist);
						 }
						 //else
							//printf("\nLocal Distance is Inconsistent!\n");
					 }
				 }
				 /*Try the subtracting dihedral angles option*/
				 if (A2 < A3) {
					 A1 = A3 - A2;
					 //printf("A1 = A3 - A2\n");
				 }
				 else {
					 A1 = A2 - A3;
					 //printf("A1 = A2 - A3\n");
				 }
				 arg_a1 = sin(c1)*sin(b1)*cos(A1) + cos(b1)*cos(c1);
				 arg_a1 = check_round(arg_a1, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
				 a1 = acos(arg_a1);
				 dist = sqrt(pow(mult_num_sols_old[col][idx_ip],2) + pow(mult_num_sols_old[col][idx_jp],2) - 2*mult_num_sols_old[col][idx_ip]*mult_num_sols_old[col][idx_jp]*cos(a1));
				 dist = check_round(dist, 1); //just put 1 here as the numerator as this isn't a quotient, and so we won't have a scenario were the quotient is infinite but should really be 0 as it's 0/0 - and a 1 will ensure it doesn't round to 0
				 //printf("New dist = %f\n",dist);
				 //printf("New dist = %f, old dist = %f\n",dist,dist_old);
				 //printf("Subtracting A2 and A3 gives dist = %f\n",dist);
				 //if ((dist >= 1) && (dist != dist_old)) { //only record or test unique distances
				 if ((dist >= 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old*pow(10,round_num1))) > 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old2*pow(10,round_num1))) > 1) && (fabs(round(dist*pow(10,round_num1)) - round(dist_old3*pow(10,round_num1))) > 1)) { //only record or test unique distances - only check that they're equivalent out to round_num1 decimal places...
					 for (z = 0; z < num_int_circs; ++z)
						 test_sol[z] = mult_num_sols_old[col][z];
					 test_sol[row] = dist;
					 consistent = consistency_check(test_sol, max_col, i, j, dist);
					 if (consistent == 1) {
						 for (z = 0; z < num_int_circs; ++z)
							 mult_num_sols_write[*num_new_sols][z] = test_sol[z];
						 //printf("\nI was called 2\n");
						 //printf("\nSolution Recorded.\n");
						 ++*num_new_sols;
						 ++num_cons_sols;
						 //printf("\nFrom Zero Canceling: Subtracting A2 and A3 found a consistent solution:\nA1 = %f, a1 = acos(%f), a1 = %f, Implies dist = %f\n",A1,arg_a1,a1,dist);
					 }
					 //else
						//printf("\nLocal Distance is Inconsistent!\n");
				 }
		 	}
		 }
	} //end of calculating distances after zeros cancelling, if appropriate.
	
	//printf("\nNum locally consistent solutions written = %d.\n",*num_new_sols);
	
	/*for (z = 0; z < *num_new_sols; ++z) {
		printf("\nSol Num %d of %d:\n",z+1,*num_new_sols);
		for (s = 0; s < num_int_circs; ++s) {
			bla1 = idx2row(s,max_col);
			bla2 = idx2col(s,bla1,max_col);
			printf("mn_write[%d][%d] = %f\n",bla1,bla2,mult_num_sols_write[z][s]);
		}
	}*/
	
	//printf("\nA3 = acos(%f), A2 = acos(%f)\n",(cos(a3) - cos(b3)*cos(c3))/(sin(c3)*sin(b3)),(cos(a2) - cos(b2)*cos(c2))/(sin(c2)*sin(b2)));
	//printf("\nA1 = %f, a1 = %f, dist = %f\n",A1,a1,dist);
	
	return num_cons_sols;
	
}

/*
Description:
Calculates the relative distance between two points, when those 2 points can be decomposed into one polygon added to another (where all the relative distances amongst those 2 polygons are known).  This function writes out the analytical expression as a string.  Later, I may link it to Maxima or write my own program so that I can fully simplify the expression and give simplified, shorter strings - but for now, it will give long, horrible expressions.

Also note that without loss of generality, we have set pk = 1 for this general rule.

Output: The relative distance between 2 points, outputted as a string.
*/
char* calc_rel_dist_analytic(int a[][MAXLEN2], int int_mat[][MAXLEN2], char **sol, const int max_col, const int num_int_circs) {

	/*Declare local variables*/
	char *dist;
	
	/*Calculate the relative distance*/
	
	return dist;

}

/*
Description: Checks that no 2 circles of intersection intersect more than twice.

Output: Returns whether the check has been violated.  Returns 1 if it has been, and 0 if it hasn't been.
*/
int check_num_ints_of_2_int_circs(int int_mat[][MAXLEN2], const int max_col, const int num_int_circs, const int only_print_sols, int *violations) {
	
	/*declare and initialize local variables*/
	int i, j, k, violated, num_points, num_intersections;
	int num_int_check[MAXLEN];
	
	violated = 0;
	num_intersections = 0;

	/*perform the check*/
	for (i = 0; i < num_int_circs; ++i) {
		for (j = 0; j < max_col; ++j) { //initialize the matrix that helps document how many intersections occur btwn 2 int circs
			num_int_check[j] = 0;
		}
		num_points = 0; //initialize
		for (j = 0; j < max_col; ++j) {
			if (int_mat[j][i] == 1) {
				num_int_check[num_points] = j;
				++num_points;
			}
		}
		for (j = i+1; j < num_int_circs; ++j) { //loop through the rest of the columns, to see if any of them share > 2 intersections
			num_intersections = 0; //initialize
			for (k = 0; k < num_points; ++k) {
				if (int_mat[num_int_check[k]][j] == 1)
					++num_intersections;
				if (num_intersections > 2) {
					if (only_print_sols == 0)
						printf("\nGraph is unphysical because it implies that 2 or more intersecting circles intersect > 2 times.\n");
					++violations[0];
					violated = 1;
					break;
				}
			}
			if (num_intersections > 2)
				break;
		}
		if (num_intersections > 2)
			break;
	}
	
	return violated;
	
}

/*
Description: 
This function checks for the following things:
1) That no more than 5 points lie on any one intersection circle.
2) Not all points on an intersection circle touch eachother (as long as there are > 2 points).
3) Records the known relative distance between the non-touching points.

Output: Returns whether the check has been violated.  Returns 1 if it has been, and 0 if it hasn't been.
*/
int check_num_points_on_int_circ(int a[][MAXLEN2], int int_mat[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col, const int num_int_circs, const int only_print_sols, int *violations) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, num_points_on_int_circ, p_int_circ_count, new_indice, pot_conn_points_count, pot_conn_points_count2, pot_conn_points_count3, deg_sep, violated, idx, row, col;
	int points_on_circs[5], pot_connected_points[13], pot_connected_points_2[13], pot_connected_points_3[13];
	
	violated = 0;
	
	/*Perform the checks*/
	for (i = 0; i < num_int_circs; ++i) {
		num_points_on_int_circ = 0; //initialize
		p_int_circ_count = 0;
		for (j = 0; j < max_col; ++j) {
			if (int_mat[j][i] == 1) { //record which points lie on each intersection circle
				points_on_circs[num_points_on_int_circ] = j;
				++num_points_on_int_circ;
			}
			if (num_points_on_int_circ > 5) {
				violated = 1;
				if (only_print_sols == 0)
					printf("\nGraph is unphysical because it implies that > 5 points lie on an intersection circle.\n");
				++violations[1];
				break;
			}
		}
		if (num_points_on_int_circ > 5)
			break;
		else { //check to see that not all points on an intersection circle are touching eachother, as this is physically impossible.
			if (num_points_on_int_circ > 2) { //Don't check if there are only 2 points on an int circ, as then it's just a tetramer, and of course the 1 relative distance btwn those 2 points can = 1
				for (j = 0; j < num_points_on_int_circ; ++j) {
					for (k = j+1; k < num_points_on_int_circ; ++k) {
						if (a[points_on_circs[j]][points_on_circs[k]] == 1)
							++p_int_circ_count;
						else {
							new_indice = 1; //initialize
							idx = row_col_2_idx(points_on_circs[j]+1,points_on_circs[k]+1,max_col);
							if (num_sol[0][idx] != DUMMY_VAR) // if the distance is already known, then it's not a new indice
								new_indice = 0;
							if (new_indice == 1) {
								/*Document the which relative distance solution code to use*/
								if (num_points_on_int_circ == 3) {
									pot_conn_points_count = 0; //initialize
									/*Check to see whether the 3 points lying on the intersection circle are connected - if they are, assign the tetrahedral distance, as it forms a tetramer.*/
									for (p = 0; p < num_points_on_int_circ; ++p) {
										if (points_on_circs[p] != points_on_circs[j] && points_on_circs[p] != points_on_circs[k]) { //don't bother checking if it corresponds to the 2 points we're already analyzing
											if (a[points_on_circs[j]][points_on_circs[p]] == 1) {
												pot_connected_points[pot_conn_points_count] = points_on_circs[p];
												++pot_conn_points_count;
											}
										}
									}
									for (p = 0; p < pot_conn_points_count; ++p) {
										if (a[pot_connected_points[p]][points_on_circs[k]] == 1) { //if one of the points connecting to the 1st point also connects to the 2nd point, then the 2 points are separated by 1 																										particle/point, and the degree of separation is 1, and the distance between them is the tetrahedral distance.
											sol[idx] = "2*sqrt(2/3)*R"; //the tetrahedral relative distance solution
											num_sol[0][idx] = 2*sqrt(2.0/3.0); //NOTE that we have to put 2.0/3.0 to get the correct floating point answer - if just have 2/3, because it's int/int, it gives something crazy
											break;
										}
									}
								} //end of assigning relative distances to the non-touching 3 points lying on an intersection circle
								else if (num_points_on_int_circ == 4) {
									pot_conn_points_count = 0; //initialize
									/*Find the degree of separation btwn the 2 points we're analyzing (i.e. how many points lie btwn them along the intersection ring on which they lie)*/
									for (p = 0; p < num_points_on_int_circ; ++p) { //find all points that connect to the first of the 2 points we're analyzing
										if (points_on_circs[p] != points_on_circs[j] && points_on_circs[p] != points_on_circs[k]) { //don't bother checking if it corresponds to the 2 points we're already analyzing
											if (a[points_on_circs[j]][points_on_circs[p]] == 1) {
												pot_connected_points[pot_conn_points_count] = points_on_circs[p];
												++pot_conn_points_count;
											}
										}
									}
									pot_conn_points_count2 = 0; //initialize
									deg_sep = 2;
									for (p = 0; p < pot_conn_points_count; ++p) {
										if (a[pot_connected_points[p]][points_on_circs[k]] == 1) { //if one of the points connecting to the 1st point also connects to the 2nd point, then the 2 points are separated by 1 																										particle/point, and the degree of separation is 1, and the distance between them is the tetrahedral distance.
											sol[idx] = "2*sqrt(2/3)*R"; //the tetrahedral relative distance solution
											num_sol[0][idx] = 2*sqrt(2.0/3.0);
											deg_sep = 1;
											break;
										}
										else {
											for (q = 0; q < num_points_on_int_circ; ++q) {
												if (points_on_circs[q] != points_on_circs[j] && points_on_circs[q] != points_on_circs[k] && points_on_circs[q] != pot_connected_points[p]) { //don't bother checking if it 																																																 corresponds to one of the 																																																 points we're already 																																																 analyzing.
													if (a[pot_connected_points[p]][points_on_circs[q]] == 1) {
														pot_connected_points_2[pot_conn_points_count2] = points_on_circs[q];
														++pot_conn_points_count2;
													}
												}
											}
										}
									} //end of first degree of separation analyses (i.e. which points are connected to the 1st point by 1 degree of separation)
									if (deg_sep > 1) { //if the degree of separation is > 1 over 4 points lying on an intersection circle, and those 4 points are connected, then the degree of separation must be 2...
										for (p = 0; p < pot_conn_points_count2; ++p) { //check to see whether those 4 points are connected
											if (a[pot_connected_points_2[p]][points_on_circs[k]] == 1) { //if one of the points connecting to a point that connects to the 1st point also connects to the 2nd point, then 																											   the 2 points are separated by 2 particles/points, and the degree of separation is 2, and the distance between 																											   them is the longest distance btwn 2 non-touching points of 4 points lying on an intersection circle.
												sol[idx] = "5/3*R"; //the relative distance btwn the 2 farthest non-touching points of 4 points on an int circ
												num_sol[0][idx] = 5.0/3.0;
												break;
											}
										}
									}
								} //end of assigning relative distances to the non-touching 4 points lying on an intersection circle
								else if (num_points_on_int_circ == 5) {
									pot_conn_points_count = 0; //initialize
									/*Find the degree of separation btwn the 2 points we're analyzing (i.e. how many points lie btwn them along the intersection ring on which they lie)*/
									for (p = 0; p < num_points_on_int_circ; ++p) { //find all points that connect to the first of the 2 points we're analyzing
										if (points_on_circs[p] != points_on_circs[j] && points_on_circs[p] != points_on_circs[k]) { //don't bother checking if it corresponds to the 2 points we're already analyzing
											if (a[points_on_circs[j]][points_on_circs[p]] == 1) {
												pot_connected_points[pot_conn_points_count] = points_on_circs[p];
												++pot_conn_points_count;
											}
										}
									}
									pot_conn_points_count2 = 0; //initialize
									deg_sep = 3;
									for (p = 0; p < pot_conn_points_count; ++p) {
										if (a[pot_connected_points[p]][points_on_circs[k]] == 1) { //if one of the points connecting to the 1st point also connects to the 2nd point, then the 2 points are separated by 1 																										particle/point, and the degree of separation is 1, and the distance between them is the tetrahedral distance.
											sol[idx] = "2*sqrt(2/3)*R"; //the tetrahedral relative distance solution
											num_sol[0][idx] = 2*sqrt(2.0/3.0);
											deg_sep = 1;
											break;
										}
										else {
											for (q = 0; q < num_points_on_int_circ; ++q) {
												if (points_on_circs[q] != points_on_circs[j] && points_on_circs[q] != points_on_circs[k] && points_on_circs[q] != pot_connected_points[p]) { //don't bother checking if it 																																																 corresponds to one of the 																																																 points we're already 																																																 analyzing.
													if (a[pot_connected_points[p]][points_on_circs[q]] == 1) {
														pot_connected_points_2[pot_conn_points_count2] = points_on_circs[q];
														++pot_conn_points_count2;
													}
												}
											}
										}
									} //end of first degree of separation analysis (i.e. which points are connected to the 1st point by 1 degree of separation)
									pot_conn_points_count3 = 0; //initialize
									if (deg_sep > 1) {
										for (p = 0; p < pot_conn_points_count2; ++p) {
											if (a[pot_connected_points_2[p]][points_on_circs[k]] == 1) { //if one of the points connecting to a point that connects to the 1st point also connects to the 2nd point, then 																											   the 2 points are separated by 2 particles/points, and the degree of separation is 2, and the distance between 																											   them is the longest distance btwn 2 non-touching points of 4 points lying on an intersection circle.
												sol[idx] = "5/3*R"; //the relative distance btwn the 2 farthest non-touching points of 4 points on an int circ
												num_sol[0][idx] = 5.0/3.0;
												deg_sep = 2;
												break;
											}
											else {
												for (q = 0; q < num_points_on_int_circ; ++q) {
													if (points_on_circs[q] != points_on_circs[j] && points_on_circs[q] != points_on_circs[k] && points_on_circs[q] != pot_connected_points_2[p]) { //don't bother checking 																																																	   if it corresponds to 																																																   one of the points we're 																																																	   already analyzing - 																																																	   could put in check to 																																																	   see that it's also not 																																																	   equal to any of the 																																																	   pot_connected_points - 																																																	   but I don't think this 																																																	   is necessary, as this 																																																	   will automatically not 																																																	   be the case...
														if (a[pot_connected_points_2[p]][points_on_circs[q]] == 1) {
															pot_connected_points_3[pot_conn_points_count3] = points_on_circs[q];
															++pot_conn_points_count3;
														}
													}
												}
											}
										}
									} //end of 2nd degree of separation analysis (i.e. which points are connected to the 2nd point by 2 degrees of separation)
									if (deg_sep > 2) { //if the degree of separation is > 2 over 5 points lying on an intersection circle, and those points are all connected, then the degree of separation must be 3...
										for (p = 0; p < pot_conn_points_count3; ++p) { //check to see whether those 5 points are connected
											if (a[pot_connected_points_3[p]][points_on_circs[k]] == 1) { //if the 4 points connect in a line on an int circ to the 5th point (this being the connection between the 4th 																											and 5th points), then the 5 points are connected, and we can apply the appropriate relative distance...
												sol[idx] = "4*sqrt(6)/9*R"; //the relative distance btwn the 2 non-touching points of 5 points on an int circ that make it so ring of 5 points is not closed
												num_sol[0][idx] = 4*sqrt(6)/9;
												break;
											}
										}
									}
								} //end of assigning solutions to the non-touching relative distances of 5 points lying on an intersection circle
							}
						} //end of statement that assigns non-touching relative distance solutions to non-connected points of a string of connected points lying on an intersection circle
						if (p_int_circ_count == num_points_on_int_circ) {
							if (only_print_sols == 0)
								printf("\nGraph is unphysical because all %d points on an intersection circle touch eachother.\n",p_int_circ_count);
							++violations[2];
							violated = 1;
							break;
						}
					}
					if (p_int_circ_count == num_points_on_int_circ)
						break;
				}
			}
		}
		if (p_int_circ_count == num_points_on_int_circ && num_points_on_int_circ > 2)
			break;
	}
	
	return violated;

}

/*
Description:
This function checks for the following things:
1) That 6 or more intersection circles don't intersect at more than 1 point.
2) If 5 intersection circles intersect at 2 points, they must lie within the same plane, and we can record the appropriate relative distance btwn the non-touching points.
3) If 4 intersection circles are associated with 4 points (i.e. a closed 4 ring) and intersect at 2 points, they must lie within the same plane, and we can record the relative distance btwn the non-touching points.

Output: Returns whether the check has been violated.  Returns 1 if it has been, and 0 if it hasn't been.
*/
int check_num_ints_of_m_int_circs(int a[][MAXLEN2], int int_mat[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col, const int num_int_circs, const int only_print_sols, int *violations) {

	/*Declare and initialize local variables*/
	int i,j, k, p, q, num_already_encountered, skip, row, col, shared_int_circs_indice, num_shared_int_circs, unique_points_count, test_point, not_unique, violated, idx;
	int shared_int_circs[MAXLEN3], already_encountered[MAXLEN3], shared_int_circs_indices[MAXLEN3], unique_points[MAXLEN3];
	
	violated = 0;
	
	/*Perform the checks*/
	num_already_encountered = 0; //initialize
	for (i = 0; i < max_col; ++i) {
		skip = 0; //initialize
		for (j = 0; j < num_already_encountered; ++j) {
			if (i == already_encountered[j])
				skip = 1;
		}
		if (skip == 0) { //only check if we haven't already established that this is one of 2 points that lies on 4, 5, or 6 intersection circles (as if > 2 such points exist, the above check for > 2 						intersections of 2 or more intersection circles would have caught it).
			shared_int_circs_indice = 0; //initialize
			for (j = 0; j < num_int_circs; ++j) {
				if (int_mat[i][j] == 1) {
					shared_int_circs[shared_int_circs_indice] = j;
					++shared_int_circs_indice;
				}
			}
			for (j = i+1; j < max_col; ++j) { //iterate through the remaining rows, to see whether any other points share 4, 5, 6, or more of the same intersection circles
				num_shared_int_circs = 0; //initialize
				for (k = 0; k < shared_int_circs_indice; ++k) {
					if (int_mat[j][shared_int_circs[k]] == 1) {
						shared_int_circs_indices[num_shared_int_circs] = shared_int_circs[k];
						++num_shared_int_circs;
					}
				}
				if (num_shared_int_circs == 4) { //if any of these scenarios are encountered, we can stop searching for more points that share 4, 5, 6, or more int circs, as if they existed, the above check 													looking for > 2 intersections of 2 or more intersection circles would have found them.  Also, we don't check the scenario where 3 intersection circles 													intersect at 2 points, as this is simply a tetramer, and will have already been picked up in the check with 3, 4, or 5 points lying on 1 int circ.
					unique_points_count = 0; //initialize
					for (k = 0; k < 4; ++k) { //determine how many unique points are associated with the 4 intersection circles
						row = idx2row(shared_int_circs_indices[k],max_col); //note that these answers for row and col are the actual particle numbers (not the zeros based particle numbers)
						col = idx2col(shared_int_circs_indices[k],row,max_col);
						for (q = 0; q <= 1; ++q) {
							not_unique = 0; //initialize
							if (q == 0)
								test_point = row;
							else
								test_point = col;
							for (p = 0; p < unique_points_count; ++p) {
								if (test_point == unique_points[p])
									not_unique = 1;
							}
							if (not_unique == 0) {
								unique_points[unique_points_count] = test_point;
								++unique_points_count;
							}
						}
					}
					if (unique_points_count == 4) { //the following checks are only valid if the 4 intersection circles are associated with only 4 unique points (if associated with > 4, this does not apply).
						if (a[i][j] == 1) { //if 2 points lie at the intersections of >= 2 circles of intersections, then the distance btwn them can not be 1.
							if (only_print_sols == 0)
								printf("\nGraph is unphysical because it implies that 4 intersection circles intersect at 2 points a distance R apart.\n");
							++violations[3];
							violated = 1;
						}
						else { //record the solution
							idx = row_col_2_idx(i+1,j+1,max_col);
							sol[idx] = "sqrt(2)*R"; //the relative distance btwn 2 intersection points of 4 mutually intersecting circles
							num_sol[0][idx] = sqrt(2);
						}
						break;
					}
				}
				else if (num_shared_int_circs == 5) {
					if (a[i][j] == 1) { //if 2 points lie at the intersections of >= 2 circles of intersections, then the distance btwn them can not be 1.
						if (only_print_sols == 0)
							printf("\nGraph is unphysical because it implies that 5 intersection circles intersect at 2 points a distance R apart.\n");
						++violations[4];
						violated = 1;
					}
					else { //record the solution
						idx = row_col_2_idx(i+1,j+1,max_col);
						sol[idx] = "sqrt(2 - 2/sqrt(5))*R";
						num_sol[0][idx] = sqrt(2 - 2/sqrt(5));
					}
					break;
				}
				else if (num_shared_int_circs >= 6) {
					if (only_print_sols == 0)
						printf("\nGraph is unphysical because it implies that %d intersection circles intersect at 2 points.\n",num_shared_int_circs);
					++violations[5];
					violated = 1;
					break;
				}
			}
		}
		if (violated == 1)
			break;
	}
	
	return violated;
	
}

/*
Description: 
Tests whether a closed 5 ring surrounds an intersection circle (if it does, this is a physical violation).  (This is different from the check where we tested for no more than 5 points on an intersection circle, as here the points need not lie on the intersection circle.)

Output: Returns whether the check has been violated.  Returns 1 if it has been, and 0 if it hasn't been.
*/
int check_closed_5_rings(int a[][MAXLEN2], const int max_col, const int only_print_sols, int *violations) {
	
	/*Declare and initialize local variables*/
	int i, j, k, l, p, q, num_already_encountered, num_closed_5_rings, num_potential_points_in_closed_ring_1, num_potential_points_in_closed_ring_2, not_encountered, check, violated;
	int points_in_closed_rings[5][MAXLEN2], encountered_points[MAXLEN4], potential_points_in_closed_ring_1[MAXLEN4], potential_points_in_closed_ring_2[MAXLEN4];
	
	violated = 0;
	
	/*Perform the check*/
	num_already_encountered = 0; //initialize
	num_closed_5_rings = 0;
	for (i = 0; i < max_col; ++i) { //first check to see whether a closed 5 ring exists
		for (j = i+1; j < max_col; ++j) { //these outer 2 nested for loops find where a_ij = 1, then we go on below and check to see if that 1 connection is part of a closed 5 ring...
			num_potential_points_in_closed_ring_1 = 0; //initialize
			num_potential_points_in_closed_ring_2 = 0;
			if (a[i][j] == 1) { //check whether there's a closed 5 ring that includes points i and j
				encountered_points[num_already_encountered] = i; //record which points have been tested
				encountered_points[num_already_encountered+1] = j;
				num_already_encountered = num_already_encountered + 2;
				for (k = j+1; k < max_col; ++k) { //iterate through remaining columns of the current row (row i) to look for a connected a_ik = 1.
					if (a[i][k] == 1) {
						potential_points_in_closed_ring_1[num_potential_points_in_closed_ring_1] = k;
						++num_potential_points_in_closed_ring_1;
					}
				}
				for (k = 0; k < max_col; ++k) { //iterate through row j to look for a connected a_jk = 1.
					if (k != i) {
						if (a[j][k] == 1) {
							potential_points_in_closed_ring_2[num_potential_points_in_closed_ring_2] = k;
							++num_potential_points_in_closed_ring_2;
						}
					}
				}
				for (k = 0; k < num_potential_points_in_closed_ring_1; ++k) { //check the newer points in the potential 5 ring, to see if they close to form the 5 ring (these would be points 3 and 4 in the 																					ring), i.e. this loop iterates through the point connected to i.
					for (l = 0; l < num_potential_points_in_closed_ring_2; ++l) { //This loop iterates through the point connected to j.
						if (potential_points_in_closed_ring_1[k] != potential_points_in_closed_ring_2[l] && a[j][potential_points_in_closed_ring_1[k]] != 1 && a[i][potential_points_in_closed_ring_2[l]] != 1) { //if the i and j are not connected to the same point (as this would imply a closed 3 ring).
							for (p = 0; p < max_col; ++p) {
								//This was messing things up, and so without this, and with the other modifications, it's finding the 5 rings, but it has a list which is not unique - so will have to put in unique check later!
								/*not_encountered = 1; //initialize
								for (q = 0; q < num_already_encountered; ++q) {
									if (p == encountered_points[q])
										not_encountered = 0;
								}*/
								//if (not_encountered == 1 && potential_points_in_closed_ring_1[k] != p && potential_points_in_closed_ring_2[l] != p) {
								
								if (potential_points_in_closed_ring_1[k] != p && potential_points_in_closed_ring_2[l] != p) {
									if (a[potential_points_in_closed_ring_1[k]][p] == 1 && a[potential_points_in_closed_ring_2[l]][p] == 1 && a[i][p] != 1 && a[j][p] != 1) { //if there exists a 5th point that touches both the 3rd and 4th points (thus closing the 5 ring), and that does not touch the 1st and 2nd point (which would close off the 5 ring).
										points_in_closed_rings[0][num_closed_5_rings] = i; //the 1st point
										points_in_closed_rings[1][num_closed_5_rings] = j; //the 2nd point
										points_in_closed_rings[2][num_closed_5_rings] = potential_points_in_closed_ring_2[l]; //the 3rd point
										points_in_closed_rings[3][num_closed_5_rings] = p; //point 4 in closed ring
										points_in_closed_rings[4][num_closed_5_rings] = potential_points_in_closed_ring_1[k]; //the 5th point
										++num_closed_5_rings; //iterate the number of closed 5 rings that exist
									}
								}
							}
						}
					}
				} 
			} //end of check for 5 ring including points i and j (i.e. for the 1st and 2nd points in the closed 5 ring).
		}
	} //end of checking for closed 5 rings
	
	/*Print closed 5 rings*/
	/*printf("\nnum_closed_5_rings = %d\n",num_closed_5_rings);
	printf("\nClosed 5 rings:\n");
	for (i = 0; i < num_closed_5_rings; ++i) {
		for (j = 0; j < 5; ++j) {
			if (j < 4)
				printf("%d",points_in_closed_rings[j][i]+1);
			else
				printf("%d\n",points_in_closed_rings[j][i]+1);
		}
	}*/
	
	/*Check for the existence of intersection circles, and then see whether a closed 5 ring surrounds it.*/
	for (j = 0; j < max_col; ++j) {
		for (k = j+1; k < max_col; ++k) {
			if (a[j][k] == 1) { //check intersection circle Ijk
				for (i = 0; i < num_closed_5_rings; ++i) { //now check to see whether any of the closed 5 rings surround the intersection circle.
					check = 1; //initialize
					for (p = 0; p < 5; ++p) {
						if (j == points_in_closed_rings[p][i] || k == points_in_closed_rings[p][i]) {
							check = 0;
							break;
						}
					}
					if (check == 1) { //only test for a violation if the originating points of the intersection circle are different from the points in the closed 5 ring
						if (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) { //check whether the 3 required bonds on the 'bottom' exist, 																																									which are needed to constrain the 5 ring
							if ((a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) || (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1) || (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1)) { //check to see whether the corresponding required 2 'upper' bonds exist
								violated = 1;
							}
						}
						if (violated != 1) {
							if (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) {
								if ((a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) || (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) || (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1)) { //check to see whether the corresponding required 2 'upper' bonds exist
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) {
								if ((a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) || (a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) || (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1)) {
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) {
								if ((a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1) || (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) || (a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1)) {
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) {
								if ((a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) || (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) || (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1)) {
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) { //check the 'flip-side' of the bottom satisfied - i.e. 																																									either side can be top or bottom
								if ((a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) || (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1) || (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1)) { //check to see whether the corresponding required 2 'upper' bonds exist
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) {
								if ((a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) || (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) || (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1)) { //check to see whether the corresponding required 2 'upper' bonds exist
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1) {
								if ((a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) || (a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) || (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1)) {
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[k][points_in_closed_rings[1][i]] == 1 && a[k][points_in_closed_rings[3][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) {
								if ((a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[2][i]] == 1) || (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) || (a[j][points_in_closed_rings[2][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1)) {
									violated = 1;
								}
							}
						}
						if (violated != 1) {
							if (a[k][points_in_closed_rings[0][i]] == 1 && a[k][points_in_closed_rings[2][i]] == 1 && a[k][points_in_closed_rings[4][i]] == 1) {
								if ((a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1) || (a[j][points_in_closed_rings[1][i]] == 1 && a[j][points_in_closed_rings[4][i]] == 1) || (a[j][points_in_closed_rings[0][i]] == 1 && a[j][points_in_closed_rings[3][i]] == 1)) {
									violated = 1;
								}
							}
						}
					}
				}
			}
			if (violated == 1) {
				if (only_print_sols == 0)
					printf("\nGraph is unphysical because it implies that a closed 5 ring surrounds a circle of intersection.\n");
				++violations[6];
				break;
			}
		}
		if (violated == 1)
			break;
	}
	
	return violated;

}

/*
Description: 
Check that 2 points lying on different triangles of a partially formed tetramer (i.e. a closed 4 ring that can form a tetramer) do not touch.  Put another way - check that 2 poins attached (by 2 bonds) to opposite sides of a closed 4 rings do not touch.  This is different from the previous check that ensured that the 2 points of 4 intersection rings intersecting at 2 points did not touch, in that the 4 intersection rings of this closed 4 ring need not intersect at 2 points (i.e. we could for example have 2 points hanging off opposite faces of a tetramer).

Output: Returns whether the check has been violated.  Returns 1 if it has been, and 0 if it hasn't been.
*/
int check_opp_points_of_4_ring_dont_touch(int a[][MAXLEN2], const int max_col, const int only_print_sols, int *violations) {
	
	/*Declare and initialize local variables*/
	
	int i, j, k, l, p, m, violated;
	
	violated = 0;
	
	/*Perform check*/
	for (i = 0; i < max_col; ++i) { //loop through rows of the adjacency matrix
		for (j = i+1; j < max_col; ++j) { //loop through columns of the adjacency matrix (that we haven't already checked)
			if (a[i][j] == 1) { //search for neighboring particles, this would be the 1st 2 points in the potential closed 4 ring
				for (k = 0; k < max_col; ++k) {
					if (k != i && k != j) {
						if (a[j][k] == 1) { //this represents the 2nd point connecting to the 3rd point in the potential closed 4 ring
							for (l = 0; l < max_col; ++l) {
								if (l != i && l != j && l != k) {
									if (a[k][l] == 1 && a[l][i] == 1) { //if the 3rd point connects to the 4th point, and if the 4th point connects to the 1st point, respectively - then we have a closed 4 ring
										for (p = 0; p < max_col; ++p) {
											if (p != i && p != j && p != k && p != l) { //check that this is a new point, and not one of the points in the closed 4 ring
												if (a[p][l] == 1 && a[p][j] == 1) {
													for (m = 0; m < max_col; ++m) { //if we have a point connected to 1 side of the closed 4 ring by 2 bonds, then continue the check to see whether we have a 2nd point 																						connected to the opposite side
														if (m != i && m != j && m != k && m != l) {
															if (a[m][i] == 1 && a[m][k] == 1) { //if this 2nd point exists, then check whether the 2 points on opposite sides of the closed 4 ring touch...
																if (a[m][p] == 1) { //if these 2 points touch, then it is a physical violation
																	//printf("Closed 4 ring is: %d,%d,%d,%d.\nPoints on opposit sides are: %d,%d\n",i,j,k,l,p,m);
																	if (only_print_sols == 0)
																		printf("\nGraph is unphysical because it implies that 2 points on opposite sides of a closed 4 ring touch.\n");
																	++violations[7];
																	violated = 1;
																	break;
																}
															}
														}
													}
												}
											}
											if (violated == 1)
												break;
										} //end of loops that check to see whether 2 points on opposite sides of a closed 4 ring touch
									}
								}
								if (violated == 1)
									break;
							}
						}
					}
					if (violated == 1)
						break;
				}
			}
			if (violated == 1)
				break;
		}
		if (violated == 1)
			break;
	}
	
	return violated;

}

/*
Description: Check adjacency matrices for the 8 particle new seed characteristic.  If they have this characteristic, then insert the 8 particle new seed solution where appropriate.

Output: No output - but the solution arrays are updated to include the 8 particle new seed solution, if appropriate.
*/
void new_seed_8_particles(int a[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, z, num_4_rings, count, document, d1_idx1, d1_idx2, d3_idx1, d3_idx2, d3_idx3, d3_idx4;
	double d1, d3;
	int points_in_4_rings[MAXLEN2][4]; //this array documents all non-closed 4 rings
	/*Here I'm puting in the numerical approximations of the analytical solutions I calculated, as the analytical solutions were horribly long and complicated to look at.  But they are contained in my notebook.  For 	now, I'll let the general relative distance rule solve for the 3rd distance of the 8 particle new seed packing (which would = d4), as that's probably what I would use analytically to derive it anyway.*/
	d1 = 1.28916854644831; //the distance between the non-closed part of the 4 rings (i.e. the 2 end points of the 4 rings)
	d3 = 1.5129998501150983; //the cross distance of the 4 rings - i.e. the diagonals of the non-closed 4 rings
	num_4_rings = 0;
	
	/*Find all non-closed 4 rings in the adjacency matrix*/
	for (i = 0; i < max_col; ++i) {
		for (j = 0; j < max_col; ++j) { //could formally below check for i != j, but a_ii = 0 anyway, so we need not do this
			if (a[i][j] == 1) { //if i and j are connected, then continue to check whether i and j are part of a non-closed 4 ring
				for (k = 0; k < max_col; ++k) {
					if ((k != i) && (k != j) && (a[j][k] == 1) && (a[i][k] == 0)) { //if j connects to k and i does not touch k (i.e. the crossterm != 1), then this might be a 4 ring and continue the check (and indeed don't really need to check k != j, as a_jj = 0 automatically)
						for (p = 0; p < max_col; ++p) {
							if ((p != i) && (p != j) && (p != k) && (a[k][p] == 1) && (a[p][j] == 0) && (a[p][i] == 0)) { //if k connects to p, and the crossterm pj is = 0, and the 2 end points, p and i, don't connect, then we have a non-closed 4 ring - once again, don't need to check p != k, as a_kk = 0, but have it here anyway...
								document = 1; //initialize
								for (q = 0; q < num_4_rings; ++q) { //loop through the 4 rings
									count = 0; //initialize
									for (z = 0; z < 4; ++z) { //loop through the points in the 4 rings
										if ((points_in_4_rings[q][z] == i) || (points_in_4_rings[q][z] == j) || (points_in_4_rings[q][z] == k) || (points_in_4_rings[q][z] == p))
											++count;
									}
									if (count == 4) { //if the count is = 4, then the 4 ring is not unique, so don't document
										document = 0;
										break;
									}
								}
								if (document == 1) { //if the 4 ring is unique (or if it is the 1st 4 ring) then document it
									points_in_4_rings[num_4_rings][0] = i;
									points_in_4_rings[num_4_rings][1] = j;
									points_in_4_rings[num_4_rings][2] = k;
									points_in_4_rings[num_4_rings][3] = p;
									++num_4_rings;
								}
							}
						}
					}
				}
			}
		}
	}
	
	/*Check to see whether the endpoints of any of the 4 rings are int points of any of the other 4 rings - i.e. if the end points are attached to all points of another 4 ring*/
	for (i = 0; i < num_4_rings; ++i) {
		for (j = i+1; j < num_4_rings; ++j) {
			if ((a[points_in_4_rings[i][0]][points_in_4_rings[j][0]] == 1) && (a[points_in_4_rings[i][0]][points_in_4_rings[j][1]] == 1) && (a[points_in_4_rings[i][0]][points_in_4_rings[j][2]] == 1) && (a[points_in_4_rings[i][0]][points_in_4_rings[j][3]] == 1) && (a[points_in_4_rings[i][3]][points_in_4_rings[j][0]] == 1) && (a[points_in_4_rings[i][3]][points_in_4_rings[j][1]] == 1) && (a[points_in_4_rings[i][3]][points_in_4_rings[j][2]] == 1) && (a[points_in_4_rings[i][3]][points_in_4_rings[j][3]] == 1)) { //if the end points of one 4 ring touches all points of another 4 ring
				if ((a[points_in_4_rings[j][0]][points_in_4_rings[i][0]] == 1) && (a[points_in_4_rings[j][0]][points_in_4_rings[i][1]] == 1) && (a[points_in_4_rings[j][0]][points_in_4_rings[i][2]] == 1) && (a[points_in_4_rings[j][0]][points_in_4_rings[i][3]] == 1) && (a[points_in_4_rings[j][3]][points_in_4_rings[i][0]] == 1) && (a[points_in_4_rings[j][3]][points_in_4_rings[i][1]] == 1) && (a[points_in_4_rings[j][3]][points_in_4_rings[i][2]] == 1) && (a[points_in_4_rings[j][3]][points_in_4_rings[i][3]] == 1)) { //if the end points of that other 4 ring also touch all the points of the 1st 4 ring, then we have an 8 particle new seed
					if (points_in_4_rings[i][0] < points_in_4_rings[i][3])
						d1_idx1 = row_col_2_idx(points_in_4_rings[i][0] + 1, points_in_4_rings[i][3] + 1, max_col); //the row and column inputs here must be 1 based, which is why I added 1 - this is the distance between the end points of the non-closed 4 ring
					else
						d1_idx1 = row_col_2_idx(points_in_4_rings[i][3] + 1, points_in_4_rings[i][0] + 1, max_col);
					if (points_in_4_rings[j][0] < points_in_4_rings[j][3])
						d1_idx2 = row_col_2_idx(points_in_4_rings[j][0] + 1, points_in_4_rings[j][3] + 1, max_col);
					else
						d1_idx2 = row_col_2_idx(points_in_4_rings[j][3] + 1, points_in_4_rings[j][0] + 1, max_col);
					if (points_in_4_rings[i][0] < points_in_4_rings[i][2])
						d3_idx1 = row_col_2_idx(points_in_4_rings[i][0] + 1, points_in_4_rings[i][2] + 1, max_col); //these are the diagonal distances of the non-closed 4 rings
					else
						d3_idx1 = row_col_2_idx(points_in_4_rings[i][2] + 1, points_in_4_rings[i][0] + 1, max_col);
					if (points_in_4_rings[i][1] < points_in_4_rings[i][3])
						d3_idx2 = row_col_2_idx(points_in_4_rings[i][1] + 1, points_in_4_rings[i][3] + 1, max_col);
					else
						d3_idx2 = row_col_2_idx(points_in_4_rings[i][3] + 1, points_in_4_rings[i][1] + 1, max_col);
					if (points_in_4_rings[j][0] < points_in_4_rings[j][2])
						d3_idx3 = row_col_2_idx(points_in_4_rings[j][0] + 1, points_in_4_rings[j][2] + 1, max_col);
					else
						d3_idx3 = row_col_2_idx(points_in_4_rings[j][2] + 1, points_in_4_rings[j][0] + 1, max_col);
					if (points_in_4_rings[j][1] < points_in_4_rings[j][3])
						d3_idx4 = row_col_2_idx(points_in_4_rings[j][1] + 1, points_in_4_rings[j][3] + 1, max_col);
					else
						d3_idx4 = row_col_2_idx(points_in_4_rings[j][3] + 1, points_in_4_rings[j][1] + 1, max_col);
					num_sol[0][d1_idx1] = d1;
					num_sol[0][d1_idx2] = d1;
					num_sol[0][d3_idx1] = d3;
					num_sol[0][d3_idx2] = d3;
					num_sol[0][d3_idx3] = d3;
					num_sol[0][d3_idx4] = d3;
					sol[d1_idx1] = "1.28916854644831*R (I've written it out numerically, as the cubic solution has imaginary numbers that cancel, and it looks horrible)";
					sol[d1_idx2] = "1.28916854644831*R";
					sol[d3_idx1] = "1.5129998501150983*R (Again, I've written out the analytical solution numerically, as it's a long, ugly expression)";
					sol[d3_idx2] = "1.5129998501150983*R";
					sol[d3_idx3] = "1.5129998501150983*R";
					sol[d3_idx4] = "1.5129998501150983*R";
				}
			}
		}
	}

}

/*
Description: See whether 3 connected 4 rings (whose crossterms also are potentially closed) are connected in a manner that implies particle overlap.  This is the rule that became relevant for 10 particles in notebook 10.  The configuration is as follows: 4 rings m,j,k,l; j,k,p,n; k,n,q,j -- where the 1st and 2nd 4 ring are connected by bond jk, the 1st and 3rd 4 ring are connected by bond kl, and the 2nd and 3rd 4 ring are connected by bond kn (so the point in common with all three 4 rings is k).  Then, a point i touches points j, l, and n (i.e. it touches all three 4 rings by touching the point of the common bond that is not in common with all three 4 rings, i.e. that is not k).  Then, a point o touches the point in common with all three 4 rings, k, and the point opposite it (m, p, or q -- which are the points of each 4 ring that are not shared).

In this configuration the following bonds can not be made: mn, oi, pl, qj -- but I won't code these in, as they will already be found by the 2 points on opposite sides of closed 4 ring rule, and the following bonds which will be coded in, can not be made: on, qm, pm, pq.  If any of these last 4 bonds are made, it implies a distance < R, and thus a physical violation.

Output: 1 if the matrix violates this rule, 0 if it doesn't.
*/
int connected_4ring_check1(int a[][MAXLEN2], int (*pot_closed_4rings)[MAXLEN5], const int num_4rings, const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, z, s, check_1i, check_2i, check_1j, check_2j, count, count2, found_point1, found_point2, shared_point, shared_point_idx, not_shared_idx, idx, i_exist, i_exist2, o_exist, o1_exist, o2_exist, o3_exist, o1_point, o2_point, o3_point, found_o1, found_o2, found_o3, num_points_common, num_points_common2, already_encountered, violated;
	int shared_points[2][4], points_i_touch[3], points_i_touch2[4], next_point_1st_2nd_rings[2][2], next_point_3rd_ring[2][2], not_touched_points[4], not_touched_idx[4], opp_shared_point[3], o_touch_points[2], o1_touch_points[2], o2_touch_points[2], o3_touch_points[2];
	violated = 0;
	
	for (i = 0; i < num_4rings; ++i) {
		//printf("\n4 ring 1:");
		//for (j = 0; j < 4; ++j)
			//printf("%d,",pot_closed_4rings[j][i]+1);
		for (j = i+1; j < num_4rings; ++j) {
			//printf("\n4 ring 2:");
			//for (k = 0; k < 4; ++k)
				//printf("%d,",pot_closed_4rings[k][j]+1);
			count = 0; //initialize
			num_points_common = 0;
			for (k = 0; k < 4; ++k) {
				for (p = 0; p < 4; ++p) {
					if (pot_closed_4rings[k][i] == pot_closed_4rings[p][j])
						++num_points_common;
				}
			}
			if (num_points_common == 2) { //only look for a common bond between the 1st two 4 rings if they share only 2 points in common
				for (k = 0; k < 4; ++k) { //look for a common bond between the 1st two 4 rings
					check_1i = k;
					if (k == 3)
						check_2i = 0;
					else
						check_2i = k + 1;
					//printf("\ncheck_1i = %d, check_2i = %d\n",check_1i,check_2i);
					for (p = 0; p < 4; ++p) {
						found_point1 = 0; //initialize
						found_point2 = 0;
						check_1j = p;
						if (p == 3)
							check_2j = 0;
						else
							check_2j = p + 1;
						//printf("\ncheck_1j = %d, check_2j = %d\n",check_1j,check_2j);
						/*if (((pot_closed_4rings[check_1i][i] == pot_closed_4rings[check_1j][j]) || (pot_closed_4rings[check_1i][i] == pot_closed_4rings[check_2j][j])) && ((pot_closed_4rings[check_2i][i] == pot_closed_4rings[check_1j][j]) || (pot_closed_4rings[check_2i][i] == pot_closed_4rings[check_2j][j]))) {
							++count;
							shared_bond1[0] = pot_closed_4rings[check_1i][i];
							shared_bond1[1] = pot_closed_4rings[check_2i][i];
						}*/
						if (pot_closed_4rings[check_1i][i] == pot_closed_4rings[check_1j][j]) {
							shared_points[0][0] = pot_closed_4rings[check_1i][i];
							shared_points[0][1] = check_1i;
							shared_points[0][2] = check_1j;
							found_point1 = 1;
						}
						else if (pot_closed_4rings[check_1i][i] == pot_closed_4rings[check_2j][j]) {
							shared_points[0][0] = pot_closed_4rings[check_1i][i];
							shared_points[0][1] = check_1i;
							shared_points[0][2] = check_2j;
							found_point1 = 1;
						}
						if (pot_closed_4rings[check_2i][i] == pot_closed_4rings[check_1j][j]) {
							shared_points[1][0] = pot_closed_4rings[check_2i][i];
							shared_points[1][1] = check_2i;
							shared_points[1][2] = check_1j;
							found_point2 = 1;
						}
						else if (pot_closed_4rings[check_2i][i] == pot_closed_4rings[check_2j][j]) {
							shared_points[1][0] = pot_closed_4rings[check_2i][i];
							shared_points[1][1] = check_2i;
							shared_points[1][2] = check_2j;
							found_point2 = 1;
						}
						if ((found_point1 == 1) && (found_point2 == 1)) {
							++count;
							//printf("\ncount = %d, shared bond = %d,%d\n",count,shared_points[0][0]+1,shared_points[1][0]+1);
						}
					}
				} //end of looking for common bond between 1st two 4 rings
				//printf("\ntotal count = %d\n",count);
				if (count == 1) {
					for (k = j+1; k < num_4rings; ++k) { //look through the 3rd 4 ring
						num_points_common = 0; //re-initialize
						num_points_common2 = 0;
						for (p = 0; p < 4; ++p) {
							for (q = 0; q < 4; ++q) {
								//printf("\n4 ring 1[%d][%d] = %d, 4 ring 2[%d][%d] = %d, 4 ring 3[%d][%d] = %d\n",p+1,i+1,pot_closed_4rings[p][i],p+1,j+1,pot_closed_4rings[p][j],q+1,k+1,pot_closed_4rings[q][k]);
								if (pot_closed_4rings[p][i] == pot_closed_4rings[q][k])
									++num_points_common;
								if (pot_closed_4rings[p][j] == pot_closed_4rings[q][k])
									++num_points_common2;
							}
						}
						//printf("\n4 ring 3:");
						//for (p = 0; p < 4; ++p)
							//printf("%d,",pot_closed_4rings[p][k]+1);
						count2 = 0; //initialize
						//printf("\nnum_points_common = %d, num_points_common2 = %d\n",num_points_common, num_points_common2);
						if ((num_points_common == 2) && (num_points_common2 == 2)) { //only continue the search if the 3rd 4 ring shares exactly 2 points in common with the 1st and 2nd 4 rings
							for (p = 0; p < 4; ++p) {
								if (pot_closed_4rings[p][k] == shared_points[0][0]) {
									shared_point = shared_points[0][0];
									shared_point_idx = 0;
									not_shared_idx = 1;
									shared_points[0][3] = p;
									points_i_touch[2] = shared_points[1][0]; //shared point between 1st and 2nd 4 rings (which is not shared by all three 4 rings)
									++count2;
								}
								if (pot_closed_4rings[p][k] == shared_points[1][0]) {
									shared_point = shared_points[1][0];
									shared_point_idx = 1;
									not_shared_idx = 0;
									shared_points[1][3] = p;
									points_i_touch[2] = shared_points[0][0]; //shared point between 1st and 2nd 4 rings (which is not shared by all three 4 rings)
									++count2;
								}
							}
							//printf("\ncount2 = %d\n",count2);
							if (count2 == 1) {
								if (shared_points[shared_point_idx][1] == 3) { //record the point next to the shared point in 1st 4 ring
									if (pot_closed_4rings[0][i] == shared_points[not_shared_idx][0]) { //if the point next to it by adding 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by substracting 1
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[2][i]; //record the point - 0 column has the points, 0 row has for ring 1, 1 row has for ring 2
										next_point_1st_2nd_rings[0][1] = 2; //record the index - 1 column has the indices, 0 row has for ring 1, 1 row has for ring 2
									}
									else { //if it's not equal to the other point of the shared bond, then this is the adjacent point we want to record
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[0][i];
										next_point_1st_2nd_rings[0][1] = 0;
									}
									opp_shared_point[0] = pot_closed_4rings[1][i]; //record point opposite to the shared point for ring 1
								}
								else if (shared_points[shared_point_idx][1] == 0) {
									if (pot_closed_4rings[4][i] == shared_points[not_shared_idx][0]) { //if the point next to it by subtracting 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by adding 1
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[1][i];
										next_point_1st_2nd_rings[0][1] = 1;
									}
									else { //if not equal, it is that point
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[3][i];
										next_point_1st_2nd_rings[0][1] = 3;
									}
									opp_shared_point[0] = pot_closed_4rings[2][i]; //record point opposite to the shared point for ring 1
								}
								else {
									if (pot_closed_4rings[shared_points[shared_point_idx][1] + 1][i] == shared_points[not_shared_idx][0]) { //if the point next to it by adding 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by substracting 1
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[shared_points[shared_point_idx][1] - 1][i];
										next_point_1st_2nd_rings[0][1] = shared_points[shared_point_idx][1] - 1;
									}
									else { //if it's not equal to the other point of the shared bond, then this is the adjacent point we want to record
										next_point_1st_2nd_rings[0][0] = pot_closed_4rings[shared_points[shared_point_idx][1] + 1][i];
										next_point_1st_2nd_rings[0][1] = shared_points[shared_point_idx][1] + 1;
									}
									if (shared_points[shared_point_idx][1] == 1)
										opp_shared_point[0] = pot_closed_4rings[3][i]; //record point opposite to the shared point for ring 1
									else if (shared_points[shared_point_idx][1] == 2)
										opp_shared_point[0] = pot_closed_4rings[0][i]; //record point opposite to the shared point for ring 1
								}
								//printf("\nopp k point 4 ring 1: %d\n",opp_shared_point[0]+1);
								if (shared_points[shared_point_idx][2] == 3) { //record the point next to the shared point in 2nd 4 ring
									if (pot_closed_4rings[0][j] == shared_points[not_shared_idx][0]) { //if the point next to it by adding 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by substracting 1
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[2][j]; //record the point
										next_point_1st_2nd_rings[1][1] = 2; //record the index
									}
									else { //if it's not equal to the other point of the shared bond, then this is the adjacent point we want to record
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[0][j];
										next_point_1st_2nd_rings[1][1] = 0;
									}
									opp_shared_point[1] = pot_closed_4rings[1][j]; //record point opposite to the shared point for ring 2
								}
								else if (shared_points[shared_point_idx][2] == 0) {
									if (pot_closed_4rings[4][j] == shared_points[not_shared_idx][0]) { //if the point next to it by subtracting 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by adding 1
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[1][j];
										next_point_1st_2nd_rings[1][1] = 1;
									}
									else { //if not equal, it is that point
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[3][j];
										next_point_1st_2nd_rings[1][1] = 3;
									}
									opp_shared_point[1] = pot_closed_4rings[2][j]; //record point opposite to the shared point for ring 2
								}
								else {
									if (pot_closed_4rings[shared_points[shared_point_idx][2] + 1][j] == shared_points[not_shared_idx][0]) { //if the point next to it by adding 1 is = to the other point in the shared bond, then the adjacent point is the point next to it by substracting 1
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[shared_points[shared_point_idx][2] - 1][j];
										next_point_1st_2nd_rings[1][1] = shared_points[shared_point_idx][2] - 1;
									}
									else { //if it's not equal to the other point of the shared bond, then this is the adjacent point we want to record
										next_point_1st_2nd_rings[1][0] = pot_closed_4rings[shared_points[shared_point_idx][2] + 1][j];
										next_point_1st_2nd_rings[1][1] = shared_points[shared_point_idx][2] + 1;
									}
									if (shared_points[shared_point_idx][2] == 1)
										opp_shared_point[1] = pot_closed_4rings[3][j]; //record point opposite to the shared point for ring 2
									else if (shared_points[shared_point_idx][2] == 2)
										opp_shared_point[1] = pot_closed_4rings[0][j]; //record point opposite to the shared point for ring 2
									//printf("\nk idx in 2nd ring = %d, 2nd 4 ring[3] = %d, opp k point in ring 2 = %d\n",next_point_1st_2nd_rings[1][1],pot_closed_4rings[3][j]+1,opp_shared_point[1]+1);
								}
								//printf("\nopp k point 4 ring 2: %d\n",opp_shared_point[1]+1);
								if (shared_points[shared_point_idx][3] == 3) { //record the point next to the shared point in the 3rd 4 ring
									next_point_3rd_ring[0][0] = pot_closed_4rings[0][k]; //record the point obtained by adding one to the index
									next_point_3rd_ring[0][1] = 0; //record the index for this point
									next_point_3rd_ring[1][0] = pot_closed_4rings[2][k]; //record the point obtained by subtracting one to the index
									next_point_3rd_ring[1][1] = 2; //record the index for this point
									opp_shared_point[2] = pot_closed_4rings[1][k]; //record point opposite to the shared point for ring 3
								}
								else if (shared_points[shared_point_idx][3] == 0) {
									next_point_3rd_ring[0][0] = pot_closed_4rings[1][k]; //record the point obtained by adding one to the index
									next_point_3rd_ring[0][1] = 1; //record the index for this point
									next_point_3rd_ring[1][0] = pot_closed_4rings[3][k]; //record the point obtained by subtracting one to the index
									next_point_3rd_ring[1][1] = 3; //record the index for this point
									opp_shared_point[2] = pot_closed_4rings[2][k]; //record point opposite to the shared point for ring 3
								}
								else {
									next_point_3rd_ring[0][0] = pot_closed_4rings[shared_points[shared_point_idx][3] + 1][k]; //record the point obtained by adding one to the index
									next_point_3rd_ring[0][1] = shared_points[shared_point_idx][3] + 1; //record the index for this point
									next_point_3rd_ring[1][0] = pot_closed_4rings[shared_points[shared_point_idx][3] - 1][k]; //record the point obtained by subtracting one to the index
									next_point_3rd_ring[1][1] = shared_points[shared_point_idx][3] - 1; //record the index for this point
									if (shared_points[shared_point_idx][3] == 1)
										opp_shared_point[2] = pot_closed_4rings[3][k]; //record point opposite to the shared point for ring 3
									else if (shared_points[shared_point_idx][3] == 2)
										opp_shared_point[2] = pot_closed_4rings[0][k]; //record point opposite to the shared point for ring 3
									//printf("\nk idx in 3rd ring = %d, 3rd 4 ring[0] = %d, opp k point in ring 3 = %d\n",next_point_3rd_ring[1][1],pot_closed_4rings[0][k]+1,opp_shared_point[2]+1);
								}
								//printf("\nopp k point 4 ring 3: %d\n",opp_shared_point[2]+1);
								found_point1 = 0; //re-initialize
								found_point2 = 0;
								if (next_point_1st_2nd_rings[0][0] == next_point_3rd_ring[0][0]) { //see if the next point in the 1st 4 ring = one of the next points in the 3rd 4 ring
									points_i_touch[1] = next_point_1st_2nd_rings[0][0]; //shared point between 1st and 3rd 4 rings (which is not shared by all three 4 rings)
									idx = 1;
									found_point1 = 1;
								}
								else if (next_point_1st_2nd_rings[0][0] == next_point_3rd_ring[1][0]) { //see if the next point in the 1st 4 ring = one of the next points in the 3rd 4 ring
									points_i_touch[1] = next_point_1st_2nd_rings[0][0]; //shared point between 1st and 3rd 4 rings (which is not shared by all three 4 rings)
									idx = 0;
									found_point1 = 1;
								}
								if (found_point1 == 1) { //if one of these points is equal, then proceed to check further
									if (next_point_1st_2nd_rings[1][0] == next_point_3rd_ring[idx][0]) { //see if the next point in the 2nd 4 ring = the left over point in the 3rd 4 ring
										points_i_touch[0] = next_point_1st_2nd_rings[1][0]; //shared point between 2nd and 3rd 4 rings (which is not shared by all three 4 rings)
										found_point2 = 1;
									}
								}
								//printf("\nfound point1 = %d, found_point2 = %d\n");
								if ((found_point1 == 1) && (found_point2 == 1)) { //if the next point on the 1st 4 ring = one of the next points on the 3rd 4 ring, and the next point on the 2nd 4 ring also = one of the next points on the 3rd 4 ring, then proceed with the check
									//printf("\npoints i touch:");
									//for (p = 0; p < 3; ++p)
										//printf("%d,",points_i_touch[p]+1);
									for (p = 0; p < max_col; ++p) { //look for the existence of point i
										already_encountered = 0; //initialize
										for (z = 0; z < 4; ++z) {
											if ((p == pot_closed_4rings[z][i]) || (p == pot_closed_4rings[z][j]) || (p == pot_closed_4rings[z][k])) { //only check if this point isn't one that's already been encountered
												already_encountered = 1;
											}
										}
										if (already_encountered == 0) {
											i_exist = check_connections(a, points_i_touch, not_touched_points, not_touched_idx, 3, 3, p);
											if (i_exist == 1) { //if i exists, look for point o
												//printf("\ni_point = %d\n",p+1);
												for (q = 0; q < 3; ++q) {
													o_touch_points[0] = shared_point;
													o_touch_points[1] = opp_shared_point[q];
													//printf("\npoints o touch:%d,%d",shared_point+1,o_touch_points[1]+1);
													for (z = 0; z < max_col; ++z) {
														already_encountered = 0; //re-initialize
														for (s = 0; s < 4; ++s) {
															if ((z == pot_closed_4rings[s][i]) || (z == pot_closed_4rings[s][j]) || (z == pot_closed_4rings[s][k]) || (z == p)) { //only check if this point isn't one that's already been encountered
																already_encountered = 1;
															}
														}
														if (already_encountered == 0) { //only check if this point hasn't already been encountered
															o_exist = check_connections(a, o_touch_points, not_touched_points, not_touched_idx, 2, 2, z);
															if (o_exist == 1) { //check that o doesn't touch n, if it joins to 1st 4 ring; l, if it joins to 2nd 4 ring; and j, if it joins to 3rd 4 ring
																//printf("\no_point = %d\n",z+1);
																if (a[z][points_i_touch[q]] == 1) { //the way I've recorded points_i_touch, if o touches the 1st ring (i.e. q = 0), it can't touch point n, which is point shared between 2nd and 3rd 4 rings (i.e. points_i_touch[0]), if o touches the 2st ring (i.e. q = 1), it can't touch point l, which is point shared between 1st and 3rd 4 rings (i.e. points_i_touch[1]), and if o touches the 3rd ring (i.e. q = 2), it can't touch point j, which is point shared between 1st and 2nd 4 rings (i.e. points_i_touch[2])
																	violated = 1; //if o touches the inappropriate point, it is a physical violation
																	//printf("\nViolated 1!\n");
																	break;
																}
																//printf("\na[%d][%d] = %d, a[%d][%d] = %d, a[%d][%d] = %d\n",opp_shared_point[0]+1,opp_shared_point[1]+1,a[opp_shared_point[0]][opp_shared_point[1]],opp_shared_point[0]+1,opp_shared_point[2]+1,a[opp_shared_point[0]][opp_shared_point[2]],opp_shared_point[1]+1,opp_shared_point[2]+1,a[opp_shared_point[1]][opp_shared_point[2]]);
																if ((a[opp_shared_point[0]][opp_shared_point[1]] == 1) || (a[opp_shared_point[0]][opp_shared_point[2]] == 1) || (a[opp_shared_point[1]][opp_shared_point[2]] == 1)) { //if any of the points that i touches touch eachother, it is a physical violation - I've written this out separately for clarity...
																	violated = 1;
																	//printf("\nViolated 2!\n");
																	break;
																}
															}
														}
													}
													if (violated == 1)
														break;
												}
												if (violated == 1)
													break;
											}
											//printf("\nopposite to k points:");
											//for (q = 0; q < 3; ++q)
												//printf("%d,",opp_shared_point[q]+1);
											for (q = 0; q < 3; ++q) {//look for alternate point i
												points_i_touch2[q] = opp_shared_point[q];
												//printf("\nCheck setting equal: opp_shared_point[%d] = %d, points_i_touch2[%d] = %d\n",q,opp_shared_point[q]+1,q,points_i_touch2[q]+1);
											}
											points_i_touch2[3] = shared_point;
											//printf("\nPoints i2 needs to touch:");
											//for (q = 0; q < 4; ++q)
												//printf("%d,",points_i_touch2[q]+1);
											//for (q = 0; q < 4; ++q)
												//printf("\na[%d][%d] = %d\n",p+1,points_i_touch2[q]+1,a[p][points_i_touch2[q]]);
											i_exist2 = check_connections(a, points_i_touch2, not_touched_points, not_touched_idx, 4, 4, p); //this is for the 2nd connected 4 ring check
											//printf("\ni_exist2 = %d\n",i_exist2);
											if (i_exist2 == 1) { //look for the alternate point o
												//printf("\nAlternate i = %d\n",p+1);
												o1_touch_points[0] = points_i_touch[1];
												o1_touch_points[1] = points_i_touch[2];
												o2_touch_points[0] = points_i_touch[0];
												o2_touch_points[1] = points_i_touch[1];
												o3_touch_points[0] = points_i_touch[0];
												o3_touch_points[1] = points_i_touch[2];
												found_o1 = 0; //initialize
												found_o2 = 0;
												found_o3 = 0;
												//printf("points o1, o2, and o3 need to touch respectively:");
												//for (z = 0; z < 2; ++z)
													//printf("\n%d,%d\n%d,%d\n%d,%d\n",o1_touch_points[0]+1,o1_touch_points[2]+1,o2_touch_points[0]+1,o2_touch_points[1]+1,o3_touch_points[0]+1,o3_touch_points[1]+1);
												for (z = 0; z < max_col; ++z) {
													already_encountered = 0; //re-initialize
													for (s = 0; s < 4; ++s) {
														if ((z == pot_closed_4rings[s][i]) || (z == pot_closed_4rings[s][j]) || (z == pot_closed_4rings[s][k]) || (z == p)) { //only check if this point isn't one that's already been encountered
															already_encountered = 1;
														}
													}
													if (already_encountered == 0) { //only check if this point hasn't already been encountered
														o1_exist = check_connections(a, o1_touch_points, not_touched_points, not_touched_idx, 2, 2, z);
														o2_exist = check_connections(a, o2_touch_points, not_touched_points, not_touched_idx, 2, 2, z);
														o3_exist = check_connections(a, o3_touch_points, not_touched_points, not_touched_idx, 2, 2, z);
														if (o1_exist == 1) {
															//printf("\no1 = %d\n",z+1);
															o1_point = z;
															found_o1 = 1;
															if ((a[z][points_i_touch[0]] == 1) || (a[z][points_i_touch2[1]] == 1) || (a[z][points_i_touch2[2]] == 1)) { //if o1 touches n, p, or q it is a violation (where n is the point touching the shared point k, and which is shared between the 2nd and 3rd 4 rings, p is the point opposite the shared point k in the 2nd 4 ring, and q is the point opposite the shared point k in the 3rd 4 ring)
																//printf("\nViolated o1\n");
																violated = 1;
																break;
															}
														}
														if (o2_exist == 1) {
															//printf("\no2 = %d\n",z+1);
															o2_point = z;
															found_o2 = 1;
															if ((a[z][points_i_touch[2]] == 1) || (a[z][points_i_touch2[0]] == 1) || (a[z][points_i_touch2[1]] == 1)) { //if o2 touches j, m, or p it is a violation (where j is the point touching the shared point k, and which is shared between the 1st and 2nd 4 rings, m is the point opposite the shared point k in the 1st 4 ring, and p is the point opposite the shared point k in the 2nd 4 ring)
																//printf("\nViolated o2\n");
																violated = 1;
																break;
															}
														}
														if (o3_exist == 1) {
															//printf("\no3 = %d\n",z+1);
															o3_point = z;
															found_o3 = 1;
															if ((a[z][points_i_touch[1]] == 1) || (a[z][points_i_touch2[0]] == 1) || (a[z][points_i_touch2[2]] == 1)) { //if o3 touches l, m, or q it is a violation (where l is the point touching the shared point k, and which is shared between the 1st and 3rd 4 rings, m is the point opposite the shared point k in the 1st 4 ring, and q is the point opposite the shared point k in the 3rd 4 ring)
																//printf("\nViolated o3\n");
																violated = 1;
																break;
															}
														}
													}
												} //end of loop looking for alternate points o
												if (violated == 1)
													break;
												if ((found_o1 == 1) && (found_o2 == 1)) { //if any of the o points touch eachother, that is also a physical violation
													if (a[o1_point][o2_point] == 1) {
														//printf("\nViolated because o1 touches o2\n");
														violated = 1;
														break;
													}
												}
												if ((found_o1 == 1) && (found_o3 == 1)) {
													if (a[o1_point][o3_point] == 1) {
														//printf("\nViolated because o1 touches o3\n");
														violated = 1;
														break;
													}
												}
												if ((found_o2 == 1) && (found_o3 == 1)) {
													if (a[o2_point][o3_point] == 1) {
														//printf("\nViolated because o2 touches o3\n");
														violated = 1;
														break;
													}
												}
											}
										}
									} //end of outer loop looking for points i
								}
							}
						}
						if (violated == 1)
							break;
					} //end of looking through 3rd 4 ring
				}
			}
			if (violated == 1)
				break;
		}
		if (violated == 1)
			break;
	}
	
	return violated;

}

/*
Description: See if an adjacency matrix contains the non-rigid 9 particle new seed (graph 8901 from 9.3nMinus6.all saved in the Research With Michael Brenner/Self Assembly/nauty22 folder).
*/
void new_seed_9part_graph8901(int a[][MAXLEN2], int (*closed_4rings)[MAXLEN5], char **sol, double (*num_sol)[MAXLEN2], const int num_4rings, const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, count, num_conn_4rings, q_exist, z_exist, m_exist, z_point, m_point, idx, num_new_seed, recorded;
	int shared_idx[4][2], conn_4ring[MAXLEN2][6], conn_4ringsWqzm[MAXLEN2][9], points_to_touch[6], points_to_touch1[4], points_to_touch2[4], not_touched_points[6], not_touched_idx[6];
	num_conn_4rings = 0;
	num_new_seed = 0;
	recorded = 0;
	
	/*Search through 4 rings for connections amongst them*/
	for (i = 0; i < num_4rings; ++i) {
		for (j = i+1; j < num_4rings; ++j) {
			count = 0; //initialize
			for (k = 0; k < 4; ++k) {
				for (p = 0; p < 4; ++p) {
					if (closed_4rings[k][i] == closed_4rings[p][j]) {
						shared_idx[count][0] = k;
						shared_idx[count][1] = p;
						++count;
					}
				}
			}
			
			if (count == 2) { //if the count = 2, we have a connected 4 ring
				conn_4ring[num_conn_4rings][1] = closed_4rings[shared_idx[0][0]][i]; //entries 1 and 4 will contain the shared indices
				conn_4ring[num_conn_4rings][4] = closed_4rings[shared_idx[1][0]][i];
				/*The connected 4 ring will be documented in the order top 3 particles (entries 0-2), bottom 3 particles (entries 3-5), where entries 0,1,5,4 are one 4 ring, entries 1,2,3,4 are the other 4 ring, and points 1 and 4 are shared between the two 4 rings.*/
				for (k = 0; k < 4; ++k) {
					if ((a[closed_4rings[k][i]][conn_4ring[num_conn_4rings][1]] == 1) && (closed_4rings[k][i] != conn_4ring[num_conn_4rings][4]))
						conn_4ring[num_conn_4rings][0] = closed_4rings[k][i];
					if ((a[closed_4rings[k][j]][conn_4ring[num_conn_4rings][1]] == 1) && (closed_4rings[k][j] != conn_4ring[num_conn_4rings][4]))
						conn_4ring[num_conn_4rings][2] = closed_4rings[k][j]; //entry 2 connects to entry 1 and 3
					if ((a[closed_4rings[k][i]][conn_4ring[num_conn_4rings][4]] == 1) && (closed_4rings[k][i] != conn_4ring[num_conn_4rings][1]))
						conn_4ring[num_conn_4rings][5] = closed_4rings[k][i]; //entry 5 connects to entry 0 and 4
					if ((a[closed_4rings[k][j]][conn_4ring[num_conn_4rings][4]] == 1) && (closed_4rings[k][j] != conn_4ring[num_conn_4rings][1]))
						conn_4ring[num_conn_4rings][3] = closed_4rings[k][j]; //entry 3 connects to entry 2 and 4
				}
				++num_conn_4rings;
			}
		}
	}
	
	/*Search through connected 4 rings for a point, q, that touches all 6 points in the connected 4 ring, and 2 points, z and m, that touch the top and bottom 3 points in the connected 4 ring as well as point q.*/
	for (i = 0; i < num_conn_4rings; ++i) {
		for (j = 0; j < 6; ++j)
			points_to_touch[j] = conn_4ring[i][j];
		for (j = 0; j < max_col; ++j) {
			q_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 6, 6, j);
			if (q_exist == 1) { //check for points z and m
				for (k = 0; k < 3; ++k) {
					points_to_touch1[k] = conn_4ring[i][k];
					points_to_touch2[k] = conn_4ring[i][k+3];
				}
				points_to_touch1[3] = j;
				points_to_touch2[3] = j;
				z_exist = 0; //initialize
				m_exist = 0;
				for (k = 0; k < max_col; ++k) {
					if (k != j) { //don't really need this check, as k = j implies one of the points that need to be touched, is not - however, no point in doing this check if k = j, so have included it...
						if (z_exist == 0) {
							z_exist = check_connections(a, points_to_touch1, not_touched_points, not_touched_idx, 4, 4, k);
							if (z_exist == 1)
								z_point = k;
						}
						if (m_exist == 0) {
							m_exist = check_connections(a, points_to_touch2, not_touched_points, not_touched_idx, 4, 4, k);
							if (m_exist == 1)
								m_point = k;
						}
					}
				}
				if ((z_exist == 1) && (m_exist == 1)) {
					for (k = 0; k < 6; ++k)
						conn_4ringsWqzm[num_new_seed][k] = conn_4ring[i][k];
					conn_4ringsWqzm[num_new_seed][6] = j;
					conn_4ringsWqzm[num_new_seed][7] = z_point;
					conn_4ringsWqzm[num_new_seed][8] = m_point;
					/*printf("\n9 particle non-rigid seed was found here.\n");
					printf("New Seed:\n");
					for (k = 0; k < 9; ++k)
						printf("%d",conn_4ringsWqzm[num_new_seed][k]+1);*/
					++num_new_seed;
				}
			}
		}
	}
	
	/*Now have to record the appropriate solutions for this non-rigid seed.  Here, I'm only recording one possible solution option - so note that I'm doing this, and only record the solution option if a distance for that relative distance has not already been found (for example, this 9 particle non-rigid new seed is part of the 10 particle rigid new seed that has one particle stabilizing one of these adjacent 4 rings - see Self Assembly notebook #10).*/
	for (i = 0; i < num_new_seed; ++i) {
		//Assign the indices for the sqrt(2)*R distance:
		//NOTE THAT FROM THIS POINT ON, I CHANGED row_col_2_idx function to itself check that row < col - thus don't have to put this check here, or anywhere else in the code from now on (and can take the check out from all places that it currently appears).
		idx = row_col_2_idx(conn_4ringsWqzm[i][2] + 1, conn_4ringsWqzm[i][4] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(2);
			//num_sol[0][idx] = 1.41421356367146;
			sol[idx] = "sqrt(2)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		idx = row_col_2_idx(conn_4ringsWqzm[i][0] + 1, conn_4ringsWqzm[i][4] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(2);
			//num_sol[0][idx] = 1.41421356367146;
			sol[idx] = "sqrt(2)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		idx = row_col_2_idx(conn_4ringsWqzm[i][5] + 1, conn_4ringsWqzm[i][1] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(2);
			//num_sol[0][idx] = 1.41421356367146;
			sol[idx] = "sqrt(2)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		idx = row_col_2_idx(conn_4ringsWqzm[i][3] + 1, conn_4ringsWqzm[i][1] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(2);
			//num_sol[0][idx] = 1.41421356367146;
			sol[idx] = "sqrt(2)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		//Now record the sqrt(3)*R distances:
		for (j = 0; j < 3; ++j) {
			idx = row_col_2_idx(conn_4ringsWqzm[i][j] + 1, conn_4ringsWqzm[i][8] + 1, max_col);
			if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
				num_sol[0][idx] = sqrt(3);
				//num_sol[0][idx] = 1.73205080756888;
				sol[idx] = "sqrt(3)*R (feeding in only 1 possible solution of non-rigid seed)";
				recorded = 1;
			}
			idx = row_col_2_idx(conn_4ringsWqzm[i][j+3] + 1, conn_4ringsWqzm[i][7] + 1, max_col);
			if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
				num_sol[0][idx] = sqrt(3);
				//num_sol[0][idx] = 1.73205080756888;
				sol[idx] = "sqrt(3)*R (feeding in only 1 possible solution of non-rigid seed)";
				recorded = 1;
			}
		}
		//Now record the sqrt(11/3)*R distances:
		idx = row_col_2_idx(conn_4ringsWqzm[i][2] + 1, conn_4ringsWqzm[i][5] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(11.0/3.0);
			//num_sol[0][idx] = 1.91485421615195;
			sol[idx] = "sqrt(11/3)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		idx = row_col_2_idx(conn_4ringsWqzm[i][0] + 1, conn_4ringsWqzm[i][3] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = sqrt(11.0/3.0);
			//num_sol[0][idx] = 1.91485421615195;
			sol[idx] = "sqrt(11/3)*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
		//Record the 2*R distance:
		idx = row_col_2_idx(conn_4ringsWqzm[i][7] + 1, conn_4ringsWqzm[i][8] + 1, max_col);
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded by another rule before this (like the octahedral rule, for example)
			num_sol[0][idx] = 2.0;
			sol[idx] = "2*R (feeding in only 1 possible solution of non-rigid seed)";
			recorded = 1;
		}
	}
	
	if (recorded == 1)
		printf("\nThe non-rigid 9 particle new seed was encountered and 1 possible solution for the non-rigid seed was entered in.  (If this graph has been found to be unphysical by the iterative packing rule, it should be checked by hand to see whether it is truly unphysical - or whether another one of the possible solutions will satisfy it).\n");
	
}

/*
Description: See if an adjacency matrix contains the non-rigid 9 particle new seed (graph 13380 from 9.3nMinus6.all saved in the Research With Michael Brenner/Self Assembly/nauty22 folder).
*/
void new_seed_9part_graph13380(int a[][MAXLEN2], int (*closed_4rings)[MAXLEN5], char **sol, double (*num_sol)[MAXLEN2], const int num_4rings, const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, z, q, q_exist, count, count2, count3, count4, new_seed_count, p_or_z_idx, q_count;
	int points_to_touch[4], not_touched_points[4], not_touched_idx[4], rings_w_points[5][MAXLEN5], points_in_common[4], points_not_in_common[4], new_seed[MAXLEN2][9], idx[6];
	count = 0;
	new_seed_count = 0;
	
	/*See if there exist three 4 rings that each have 1 point touching all 4 points in the 4 ring (these points are called m, n, q)*/
	for (i = 0; i < num_4rings; ++i) {
		for (j = 0; j < 4; ++j)
			points_to_touch[j] = closed_4rings[j][i];
		q_count = 0; //initialize
		for (j = 0; j < max_col; ++j) {
			q_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 4, 4, j);
			if (q_exist == 1) {
				++q_count;
				for (k = 0; k < 4; ++k)
					rings_w_points[k][count] = closed_4rings[k][i];
				rings_w_points[4][count] = j;
				++count;
				if (q_count > 1)
					--count; //if more than one point touches the 4 points in the ring, decrease the count, and don't record this ring or the points - as it can't be a part of this seed (we decrease the count, as if q_count >= 2, then at least one point has been recorded, and each time this is encountered, the count should be decreased so as to record over this point.
			}
		}
	}
	
	/*Now look to see whether these 4 rings each share 2 points in common with the other two 4 rings*/
	if (count >= 3) { //if there are >= 3 closed 4 rings, each with a point touching all points in the 4 ring, then we may have this new seed - so continue checking
		for (i = 0; i < count; ++i) {
			for (j = i+1; j < count; ++j) {
				count2 = 0; //initialize
				count3 = 0;
				for (z = 0; z < 4; ++z) {
					for (k = 0; k < 4; ++k) {
						if (rings_w_points[z][i] == rings_w_points[k][j]) { //record the points between the two 4 rings that are shared
							points_in_common[count2] = rings_w_points[z][i];
							++count2;
						}
					}
				}
				if (count2 == 2) {
					count3 = 0;
					for (k = 0; k < 4; ++k) { //document the points that are not shared, as the 3rd 4 ring will have to touch all of these points
						if ((rings_w_points[k][i] != points_in_common[0]) && (rings_w_points[k][i] != points_in_common[1])) { //record the points in the first 4 ring that aren't shared
							points_not_in_common[count3] = rings_w_points[k][i];
							++count3;
						}
					}
					for (k = 0; k < 4; ++k) { //have 2 different for loops, so that we record points in the same 4 loop together
						if ((rings_w_points[k][j] != points_in_common[0]) && (rings_w_points[k][j] != points_in_common[1]) && (a[rings_w_points[k][j]][points_not_in_common[count3-1]] == 1)) { //record the points in the second 4 ring that aren't shared - 1st record the point that touches the last point recorded
							points_not_in_common[count3] = rings_w_points[k][j];
							++count3;
							break;
						}
					}
					for (k = 0; k < 4; ++k) { //have 2 different for loops, so that we record points in the same 4 loop together
						if ((rings_w_points[k][j] != points_in_common[0]) && (rings_w_points[k][j] != points_in_common[1]) && (a[rings_w_points[k][j]][points_not_in_common[count3-1]] == 1)) { //record the points in the second 4 ring that aren't shared - now record the point that touches this one
							points_not_in_common[count3] = rings_w_points[k][j];
							++count3;
							break;
						}
					}
					if (count3 == 4) {
						for (k = j+1; k < count; ++k) {
							count4 = 0;
							for (z = 0; z < 4; ++z) {
								for (q = 0; q < 4; ++q) {
									if (rings_w_points[z][k] == points_not_in_common[q])
										++count4;
								}
							}
							if (count4 == 4) { //then we have the new seed, record the points in the structure
								for (z = 0; z < 4; ++z)
									new_seed[new_seed_count][z] = points_not_in_common[z];
								new_seed[new_seed_count][4] = rings_w_points[4][k];
								for (z = 0; z < 2; ++z) { //for p, z points, first record the point (p or z) that touches the 2 points recorded in array positions 1 and 2 of the new seed array (and of the 1st 4 ring)
									if ((a[new_seed[new_seed_count][1]][points_in_common[z]] == 1) && (a[new_seed[new_seed_count][2]][points_in_common[z]] == 1)) {
										new_seed[new_seed_count][5] = points_in_common[z];
										if (z == 0)
											p_or_z_idx = 1;
										else if (z == 1)
											p_or_z_idx = 0;
										break;
									}
								}
								new_seed[new_seed_count][6] = points_in_common[p_or_z_idx]; //then record the p or z point that touches the 0th and 3rd positions of the new seed array (and necessarily of the 1st 4 ring)
								new_seed[new_seed_count][7] = rings_w_points[4][i];
								new_seed[new_seed_count][8] = rings_w_points[4][j];
								++new_seed_count;
							}
						}
					}
				}
			}
		}
	}
	
	/*Record the relative distance solutions*/
	for (i = 0; i < new_seed_count; ++i) {
		//printf("\n9 Particle Graph 13380 New Seed Found\n");
		idx[0] = row_col_2_idx(new_seed[i][7] + 1, new_seed[i][4] + 1, max_col); //record r_mq = r_mn = r_qn distances
		idx[1] = row_col_2_idx(new_seed[i][7] + 1, new_seed[i][8] + 1, max_col);
		idx[2] = row_col_2_idx(new_seed[i][4] + 1, new_seed[i][8] + 1, max_col);
		for (j = 0; j < 3; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.724744871392;
				sol[idx[j]] = "1.724744871392*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[i][7] + 1, new_seed[i][2] + 1, max_col); //record r_mk = r_mj = r_qz = r_qp = r_nl = r_ni distances
		idx[1] = row_col_2_idx(new_seed[i][7] + 1, new_seed[i][3] + 1, max_col);
		idx[2] = row_col_2_idx(new_seed[i][4] + 1, new_seed[i][6] + 1, max_col);
		idx[3] = row_col_2_idx(new_seed[i][4] + 1, new_seed[i][5] + 1, max_col);
		idx[4] = row_col_2_idx(new_seed[i][8] + 1, new_seed[i][1] + 1, max_col);
		idx[5] = row_col_2_idx(new_seed[i][8] + 1, new_seed[i][0] + 1, max_col);
		for (j = 0; j < 6; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.650680123886;
				sol[idx[j]] = "1.650680123886*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[i][1] + 1, new_seed[i][6] + 1, max_col); //record r_lz = r_lj = r_zk = r_ki = r_ip = r_pj distances
		idx[1] = row_col_2_idx(new_seed[i][1] + 1, new_seed[i][3] + 1, max_col);
		idx[2] = row_col_2_idx(new_seed[i][6] + 1, new_seed[i][2] + 1, max_col);
		idx[3] = row_col_2_idx(new_seed[i][2] + 1, new_seed[i][0] + 1, max_col);
		idx[4] = row_col_2_idx(new_seed[i][0] + 1, new_seed[i][5] + 1, max_col);
		idx[5] = row_col_2_idx(new_seed[i][5] + 1, new_seed[i][3] + 1, max_col);
		for (j = 0; j < 6; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.414213562373;
				sol[idx[j]] = "1.414213562373*R (found by Newton Iterations) (this = sqrt(2)*R)";
			}
		}
	}

}

/*
Description: See if an adjacency matrix contains the non-rigid 9 particle new seed (graph 8926 from 9.3nMinus6.all saved in the Research With Michael Brenner/Self Assembly/nauty22 folder).
*/
void new_seed_9part_graph8926(int a[][MAXLEN2], int int_mat[][MAXLEN2], int (*closed_5rings)[MAXLEN5], char **sol, double (*num_sol)[MAXLEN2], const int num_5rings, const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q, q_exist, p_exist, adjacent_part_idx1, adjacent_part_idx2, polytet_count, count, new_seed_count, next_part1, next_part2, idx, minus1called;
	int points_to_touch[5], not_touched_points[5], not_touched_idx[5], int_circ_indices[2], polytet[MAXLEN2][4], new_seed[MAXLEN2][11];
	polytet_count = 0;
	new_seed_count = 0;
	
	/*Find all instances of this new seed*/
	for (i = 0; i < num_5rings; ++i) {
		for (j = 0; j < 5; ++j)
			points_to_touch[j] = closed_5rings[j][i];
		for (j = 0; j < max_col; ++j) { //see if there exists a point, q, that touches 5/5 points in the closed 5 ring
			q_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 5, 5, j);
			if (q_exist == 1) {
				for (k = 0; k < max_col; ++k) { //look for a point, p, that touches 4/5 points of the same closed 5 ring
					if (k != j) { //here this check is necessary, as j necessarily satisfies this k condition also
						p_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 5, 4, k);
						if (p_exist == 1) {
							if (not_touched_idx[0] == 4) {
								adjacent_part_idx1 = 0;
								adjacent_part_idx2 = 3;
							}
							else if (not_touched_idx[0] == 0) {
								adjacent_part_idx1 = 1;
								adjacent_part_idx2 = 4;
							}
							else {
								adjacent_part_idx1 = not_touched_idx[0] + 1;
								adjacent_part_idx2 = not_touched_idx[0] - 1;
							}
							int_circ_indices[0] = row_col_2_idx(j+1, points_to_touch[adjacent_part_idx1] + 1, max_col);
							int_circ_indices[1] = row_col_2_idx(j+1, points_to_touch[adjacent_part_idx2] + 1, max_col);
							for (q = 0; q < 2; ++q) {
								count = 0; //initialize
								for (p = 0; p < max_col; ++p) { //check to see whether a 6 particle polytetrahedron is attached to this 5 ring, with the broken bond point as a point at the end of the polytetrahedron
									if (int_mat[p][int_circ_indices[q]] == 1) {
										polytet[polytet_count][count] = p;
										++count;
									}
								}
								if (count == 4) { //if the polytetrahedron exists, then document the new seed in an array of the following order: broken bond, adjacent particle 1, point that touches adj part1, point that touches this point and adj part2, adjacent particle 2 (this completes the closed 5 ring), q particle, p particle, polytetrahedron point that touches the broken bond, point that touches this point, point that touches that point, and point that touches that point and connects back to the closed 5 ring.
									//record the points in the 5 ring into the new seed array in the order
									new_seed[new_seed_count][0] = not_touched_points[0]; //the 'broken bond' the point in the 5 ring that the p particle does not touch
									new_seed[new_seed_count][7] = new_seed[new_seed_count][0]; //record the particle in the polytetrahedron that is the broken bond
									for (p = 0; p < 4; ++p) { //record the particle in the polytetrahedron that touches the broken bond
										if (a[polytet[polytet_count][p]][new_seed[new_seed_count][7]] == 1)
											new_seed[new_seed_count][8] = polytet[polytet_count][p];
									}
									for (p = 0; p < 4; ++p) { //record the particle in the polytetrahedron that touches the last particle recorded
										if ((a[polytet[polytet_count][p]][new_seed[new_seed_count][8]] == 1) && (polytet[polytet_count][p] != new_seed[new_seed_count][7]))
											new_seed[new_seed_count][9] = polytet[polytet_count][p];
									}
									for (p = 0; p < 4; ++p) { //record the particle in the polytetrahedron that touches the last particle recorded
										if ((a[polytet[polytet_count][p]][new_seed[new_seed_count][9]] == 1) && (polytet[polytet_count][p] != new_seed[new_seed_count][8]))
											new_seed[new_seed_count][10] = polytet[polytet_count][p];
									}
									if (a[closed_5rings[adjacent_part_idx1][i]][new_seed[new_seed_count][8]] == 1) {//record the points to the left and right of the broken bond (record the part of the 5 ring that connects to the polytetrahedron first)
										new_seed[new_seed_count][1] = closed_5rings[adjacent_part_idx1][i]; //the point to the right of the broken bond (that is in the 5 ring and touches the polytetrahedron)
										new_seed[new_seed_count][4] = closed_5rings[adjacent_part_idx2][i]; //the point to the left of the broken bond (that is in the 5 ring but doesn't touch the polytetrahedron in question)
										idx = adjacent_part_idx1;
									}
									else {
										new_seed[new_seed_count][1] = closed_5rings[adjacent_part_idx2][i];
										new_seed[new_seed_count][4] = closed_5rings[adjacent_part_idx1][i];
										idx = adjacent_part_idx2;
									}
									if (idx == 4)
										next_part1 = 0;
									else
										next_part1 = idx + 1;
									minus1called = 0;
									if (next_part1 == not_touched_idx[0]) {
										if (idx == 0)
											next_part1 = 4;
										else
											next_part1 = idx - 1;
										minus1called = 1;
									}
									new_seed[new_seed_count][2] = closed_5rings[next_part1][i]; //the particle that touches the previous one
									if (minus1called == 0) {
										if (next_part1 == 4)
											next_part2 = 0;
										else
											next_part2 = next_part1 + 1;
									}
									else {
										if (next_part1 == 0)
											next_part2 = 4;
										else
											next_part2 = next_part1 - 1;
									}
									new_seed[new_seed_count][3] = closed_5rings[next_part2][i]; //the particle that touches the previous one
									new_seed[new_seed_count][5] = j; //q particle that touches 5/5 points in the 5 ring
									new_seed[new_seed_count][6] = k; //p particle that touches 4/5 points in the 5 ring
									++polytet_count; //update counts
									++new_seed_count;
								}
							}
						}
					}
				}
			}
		}
	} //end of finding all instances where this new seed occurs
	
	/*Record relative distance solutions (prob need to record above 5 ring in specific order)*/
	for (i = 0; i < new_seed_count; ++i) {
		//printf("\n9 Particle Graph 8926 New Seed Found\n");
		idx = row_col_2_idx(new_seed[i][6] + 1, new_seed[i][9] + 1, max_col); //r_pk
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.67962000192221;
			sol[idx] = "1.67962000192221*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][6] + 1, new_seed[i][7] + 1, max_col); //r_pm
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.18324624211465;
			sol[idx] = "1.18324624211465*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][6] + 1, new_seed[i][8] + 1, max_col); //r_pz
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.74961891637700;
			sol[idx] = "1.74961891637700*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][6] + 1, new_seed[i][5] + 1, max_col); //r_pq
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.12722055426574;
			sol[idx] = "1.12722055426574*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][9] + 1, new_seed[i][3] + 1, max_col); //r_ki
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.69379006112633;
			sol[idx] = "1.69379006112633*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][9] + 1, new_seed[i][4] + 1, max_col); //r_kn
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.96667098490044;
			sol[idx] = "1.96667098490044*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][7] + 1, new_seed[i][3] + 1, max_col); //r_mi
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.66666666666667;
			sol[idx] = "1.66666666666667*R (found by Newton Iterations) (so this = 5/3*R)";
		}
		idx = row_col_2_idx(new_seed[i][3] + 1, new_seed[i][8] + 1, max_col); //r_iz
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.98383251254661;
			sol[idx] = "1.98383251254661*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][3] + 1, new_seed[i][1] + 1, max_col); //r_il
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.59199949296009;
			sol[idx] = "1.59199949296009*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][2] + 1, new_seed[i][4] + 1, max_col); //r_jn
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.59199949296009;
			sol[idx] = "1.59199949296009*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][8] + 1, new_seed[i][4] + 1, max_col); //r_zn
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.71671702138068;
			sol[idx] = "1.71671702138068*R (found by Newton Iterations)";
		}
		idx = row_col_2_idx(new_seed[i][4] + 1, new_seed[i][1] + 1, max_col); //r_nl
		if (num_sol[0][idx] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx] = 1.53446238558517;
			sol[idx] = "1.53446238558517*R (found by Newton Iterations)";
		}
	}

}

/*
Description: See if an adjacency matrix contains the non-rigid 9 particle new seed (graph 12811 from 9.3nMinus6.all saved in the Research With Michael Brenner/Self Assembly/nauty22 folder).
*/
void new_seed_9part_graph12811(int a[][MAXLEN2], int (*closed_5rings)[MAXLEN5], char **sol, double (*num_sol)[MAXLEN2], const int num_5rings, const int max_col) {

	/*Declare and initialize all local variables*/
	int i, j, k, p, z_exist, q_exist, m_exist, n_exist, adjacent_part_idx1, adjacent_part_idx2, next_part1, next_part2, found, m_point, n_point, new_seed_count;
	int points_to_touch[5], points_to_touch1[3], points_to_touch2[3], not_touched_points[5], not_touched_idx[5], test_seed[9], new_seed[9][MAXLEN5], idx[4]; 
	new_seed_count = 0;
	
	/*Find all instances of this new seed*/
	for (i = 0; i < num_5rings; ++i) {
		for (j = 0; j < 5; ++j)
			points_to_touch[j] = closed_5rings[j][i];
		for (j = 0; j < max_col; ++j) { //see if a point, z, that touches all points in the closed 5 ring
			z_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 5, 5, j);
			if (z_exist == 1) {
				for (k = 0; k < max_col; ++k) { //see if a point, q, that touches 4/5 points in the closed 5 ring exists
					if (k != j) { //here this check is necessary, as j necessarily satisfies this k condition also
						q_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 5, 4, k);
						if (q_exist == 1) {
							if (not_touched_idx[0] == 4) {
								adjacent_part_idx1 = 0;
								adjacent_part_idx2 = 3;
							}
							else if (not_touched_idx[0] == 0) {
								adjacent_part_idx1 = 1;
								adjacent_part_idx2 = 4;
							}
							else {
								adjacent_part_idx1 = not_touched_idx[0] + 1;
								adjacent_part_idx2 = not_touched_idx[0] - 1;
							}
							points_to_touch1[0] = not_touched_points[0];
							points_to_touch1[1] = k;
							points_to_touch1[2] = points_to_touch[adjacent_part_idx1];
							points_to_touch2[0] = not_touched_points[0];
							points_to_touch2[1] = k;
							points_to_touch2[2] = points_to_touch[adjacent_part_idx2];
							found = 0; //initialize
							for (p = 0; p < max_col; ++p) { //see if points, m and n, that each touch the broken bond, particle q, and one of the particles adjacent to the broken bond, exists
								if (p != j) {
									m_exist = check_connections(a, points_to_touch1, not_touched_points, not_touched_idx, 3, 3, p);
									if (m_exist == 1) {
										m_point = p; //it doesn't matter whether I call these points m or n, as the structure is symmetric
										++found;
									}
									n_exist = check_connections(a, points_to_touch2, not_touched_points, not_touched_idx, 3, 3, p);
									if (n_exist == 1) {
										n_point = p; //it doesn't matter whether I call these points m or n, as the structure is symmetric
										++found;
									}
								}
							}
							if ((found == 2) && (m_point != n_point)) { //in this case, this 9 particle new seed exists, and we should record the points here - because this strucutre is symmetric, I should check for uniqueness
								//record the closed 5 ring in the new seed in the order: broken bond, adjacent particle to right, particle that touches that, particle that touches this and adj part to left of broken bond, adjacent particle to left of broken bond.
								test_seed[0] = not_touched_points[0]; //the 'broken bond' the point in the 5 ring that the p particle does not touch
								test_seed[1] = closed_5rings[adjacent_part_idx1][i]; //the point to the right of the broken bond (i.e. the broken bond index plus one)
								if (adjacent_part_idx1 == 4)
									next_part1 = 0;
								else
									next_part1 = adjacent_part_idx1 + 1;
								test_seed[2] = closed_5rings[next_part1][i]; //the particle that touches the previous one
								if (next_part1 == 4)
									next_part2 = 0;
								else
									next_part2 = next_part1 + 1;
								test_seed[3] = closed_5rings[next_part2][i]; //the particle that touches the previous one
								test_seed[4] = closed_5rings[adjacent_part_idx2][i]; //the point to the left of the broken bond (i.e. the broken bond index minus one)
								test_seed[5] = j; //now record particle z, that touches all points in the 5 ring
								test_seed[6] = k; //now record particle q, that touches 4/5 points in the 5 ring
								test_seed[7] = m_point; //now record particle m
								test_seed[8] = n_point; //now record particle n
								unique_ring_check(new_seed, test_seed, 9, &new_seed_count); //record only unique incidences of this new seed
							}
						}
					}
				}
			}
		}
	} //end of finding all instances of the new seed
	
	/*Record the relative distances*/
	for (i = 0; i < new_seed_count; ++i) {
		//printf("\n9 Particle Graph 12811 New Seed Found\n");
		idx[0] = row_col_2_idx(new_seed[1][i] + 1, new_seed[3][i] + 1, max_col); //record the r_jl = r_ln = r_ik = r_km distances
		idx[1] = row_col_2_idx(new_seed[1][i] + 1, new_seed[8][i] + 1, max_col);
		idx[2] = row_col_2_idx(new_seed[2][i] + 1, new_seed[4][i] + 1, max_col);
		idx[3] = row_col_2_idx(new_seed[4][i] + 1, new_seed[7][i] + 1, max_col);
		for (j = 0; j < 4; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.607068616748;
				sol[idx[j]] = "1.607068616748*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[7][i] + 1, new_seed[2][i] + 1, max_col); //record the r_jn = r_im distances
		idx[1] = row_col_2_idx(new_seed[8][i] + 1, new_seed[3][i] + 1, max_col);
		for (j = 0; j < 2; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.698091312875;
				sol[idx[j]] = "1.698091312875*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[7][i] + 1, new_seed[3][i] + 1, max_col); //record the r_jm = r_in distances
		idx[1] = row_col_2_idx(new_seed[8][i] + 1, new_seed[2][i] + 1, max_col);
		for (j = 0; j < 2; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.970663367209;
				sol[idx[j]] = "1.970663367209*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[2][i] + 1, new_seed[0][i] + 1, max_col); //record the r_jp = r_nz = r_ip = r_mz distances
		idx[1] = row_col_2_idx(new_seed[7][i] + 1, new_seed[5][i] + 1, max_col);
		idx[2] = row_col_2_idx(new_seed[3][i] + 1, new_seed[0][i] + 1, max_col);
		idx[3] = row_col_2_idx(new_seed[8][i] + 1, new_seed[5][i] + 1, max_col);
		for (j = 0; j < 4; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.642586774839;
				sol[idx[j]] = "1.642586774839*R (found by Newton Iterations)";
			}
		}
		idx[0] = row_col_2_idx(new_seed[1][i] + 1, new_seed[4][i] + 1, max_col); //record the r_lk distance
		if (num_sol[0][idx[0]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
			num_sol[0][idx[0]] = 1.582669538935;
			sol[idx[0]] = "1.582669538935*R (found by Newton Iterations)";
		}
		idx[0] = row_col_2_idx(new_seed[5][i] + 1, new_seed[6][i] + 1, max_col); //record the r_zq = r_pq distances
		idx[1] = row_col_2_idx(new_seed[0][i] + 1, new_seed[6][i] + 1, max_col);
		for (j = 0; j < 2; ++j) {
			if (num_sol[0][idx[j]] == DUMMY_VAR) { //only record the solution if a solution for this distance has not already been recorded
				num_sol[0][idx[j]] = 1.085261702856;
				sol[idx[j]] = "1.085261702856*R (found by Newton Iterations)";
			}
		}
	}

}

/*
Description: See if an adjacency matrix contains the non-rigid 9 particle new seed (graph 8974 from 9.3nMinus6.all saved in the Research With Michael Brenner/Self Assembly/nauty22 folder).
*/
int new_seed_9part_graph8974(int a[][MAXLEN2], int (*closed_6rings)[MAXLEN5], const int num_6rings, const int max_col) {

	/*Declare and initialize local variables*/
	int i, j, k, p, q_exist, i_exist, p_exist, adjacent_part_idx1, adjacent_part_idx2, violated;
	int points_to_touch[6], points_to_touch2[4], not_touched_points[6], not_touched_idx[6];
	violated = 0;
	
	/*Look for point, q, that touches all points in the 6 ring, and point, i, that touches 5/6 points in the 6 ring*/
	for (i = 0; i < num_6rings; ++i) {
		for (j = 0; j < 6; ++j)
			points_to_touch[j] = closed_6rings[j][i];
		for (j = 0; j < max_col; ++j) {
			q_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 6, 6, j);
			if (q_exist == 1) {
				for (k = 0; k < max_col; ++k) { //look for point i
					if (k != j) { //here this check is necessary, as j necessarily satisfies this k condition also
						i_exist = check_connections(a, points_to_touch, not_touched_points, not_touched_idx, 6, 5, k);
						if (i_exist == 1) {
							if (not_touched_idx[0] == 5) {
								adjacent_part_idx1 = 0;
								adjacent_part_idx2 = 4;
							}
							else if (not_touched_idx[0] == 0) {
								adjacent_part_idx1 = 1;
								adjacent_part_idx2 = 5;
							}
							else {
								adjacent_part_idx1 = not_touched_idx[0] + 1;
								adjacent_part_idx2 = not_touched_idx[0] - 1;
							}
							/*Check to see whether a point, p, that touches the broken bond, the 2 adjacent particles, and the i particle exists*/
							points_to_touch2[0] = not_touched_points[0];
							points_to_touch2[1] = closed_6rings[adjacent_part_idx1][i];
							points_to_touch2[2] = closed_6rings[adjacent_part_idx2][i];
							points_to_touch2[3] = k;
							for (p = 0; p < max_col; ++p) {
								if ((p != j) && (p != k)) {
									p_exist = check_connections(a, points_to_touch2, not_touched_points, not_touched_idx, 4, 4, p);
									if (p_exist == 1) {
										violated = 1;
										break;
									}
								}
							}
						}
					}
					if (violated == 1)
						break;
				}
			}
			if (violated == 1)
				break;
		}
		if (violated == 1)
			break;
	}
	
	return violated;

}

int pot_closed_4ring_checks(int a[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col, const int only_print_sols, int *violations) {
	
	/*Declare and initialize local variables*/
	int violated;
	int num_4rings;
	int pot_closed_4rings[4][MAXLEN5];
	violated = 0;
	
	/*Find all closed 4 rings*/
	num_4rings = find_4rings_pot_closed(a, pot_closed_4rings, max_col);
	
	violated = connected_4ring_check1(a, pot_closed_4rings, num_4rings, max_col);
	if (violated == 1) {
		//if (only_print_sols == 0)
			//printf("\nGraph is unphysical because points in a series of 4 rings touch eachother, that imply a distance < R.\n");
		++violations[10];
	}
	
	return violated;	
	
}

/*
Description: Checks an adjacency matrix for new seeds of n >= 9 particles that involve closed 4 rings.

Output: No output - but the solution array is updated.
*/
void new_seed_4ring_checks(int a[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col) {

	/*Declare and initialize local variables*/
	int num_closed4rings;
	int closed_4rings[4][MAXLEN5];
	
	/*Find all closed 4 rings*/
	num_closed4rings = find_4rings(a, closed_4rings, max_col);
	
	/*Check for a 9 particle rigid new seed (graph 13380)*/
	new_seed_9part_graph13380(a, closed_4rings, sol, num_sol, num_closed4rings, max_col);
	
	/*Check for the 9 particle non-rigid new seed (graph 8901)*/
	new_seed_9part_graph8901(a, closed_4rings, sol, num_sol, num_closed4rings, max_col);

}

/*
Description: Checks an adjacency matrix for new seeds of n >= 9 particles that involve closed 5 rings.

Output: No output - but the solution array is updated.
*/
void new_seed_5ring_checks(int a[][MAXLEN2], int int_mat[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col) {

	/*Declare and initialize local variables*/
	int num_closed5rings;
	int closed_5rings[5][MAXLEN5];
	
	/*Find all closed 5 rings*/
	num_closed5rings = find_5rings(a, closed_5rings, max_col);
	
	/*Check for a 9 particle rigid new seed (graph 8926)*/
	new_seed_9part_graph8926(a, int_mat, closed_5rings, sol, num_sol, num_closed5rings, max_col);
	
	/*Check for a 9 particle rigid new seed (graph 12811)*/
	new_seed_9part_graph12811(a, closed_5rings, sol, num_sol, num_closed5rings, max_col);

}

/*
Description: Checks an adjacency matrix for new seeds of n >= 9 particles that involve closed 6 rings.

Output: 1 if the adjacency matrix contains an unphysical new seed - 0 if the adjacency matrix does not contain an unphysical new seed.  Also, the solution array is updated, if necessary.
*/
int new_seed_6ring_checks(int a[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], const int max_col, const int only_print_sols, int *violations) {

	/*Declare and initialize local variables*/
	int num_closed6rings, violated;
	int closed_6rings[6][MAXLEN5];
	violated = 0;
	
	/*Find all closed 6 rings*/
	num_closed6rings = find_6rings(a, closed_6rings, max_col);
	
	/*Check for the 9 particle unphysical new seed (graph 8974)*/
	violated = new_seed_9part_graph8974(a, closed_6rings, num_closed6rings, max_col);
	if (violated == 1) {
		if (only_print_sols == 0)
			printf("\nGraph is unphysical because it contains a 9 particle unphysical new seed (graph 8974), and thus implies a distance < R.\n");
		++violations[9];
	}
		
	
	return violated;

}

/*
Description: Calculate a new relative distance r_ij, where the points j and i can be decomposed into 2 points associated with 2 polygons, where one polygon is added to another, and they share a common base, kpq.

Output: Returns 1 if a violation has occured, and 0 if a violation has not occured.  The violation that can occur from this check is a relative distance being < R.
*/
int calc_new_rel_dists(int a[][MAXLEN2], int int_mat[][MAXLEN2], char **sol, double (*num_sol)[MAXLEN2], int *mult_num_sol_count, int *valid_sols_count, const int max_col, const int num_int_circs, const int only_print_sols, int *violations) {

	/*Declare and initialize local variables*/
	int i, j, k, q, l, m, p, z, stop, num_solved_dists, num_common_points, point1, point2, idx1, idx2, idx3, idx_check, found_base, consistent, temp, num_cons_sols, num_new_sols, temp_count, unknown_var, violated;
	int common_points[MAXLEN2], base[3];
	double mult_sol_num_temp[MAXLEN5][MAXLEN2], temp_sol[num_int_circs];
	
	stop = 0;
	violated = 0;
	temp = 0; //if temp = 0, we should use mult_sol_num; if temp = 1, we should use mult_sol_num_temp
	num_new_sols = 0;
	temp_count = 0;
	unknown_var = 0;
	for (i = 0; i < num_int_circs; ++i)
		temp_sol[i] = num_sol[0][i]; //can always set it to the 1st column, as we only want to know what distances are known in the initial check (to see if we have enough info to calculate the new distance)
	
	/*Calculate unknown relative distances, if they can be decomposed into one polygon added to another*/
	while (stop == 0) {
		num_solved_dists = 0; //initialize
		violated = 0; //re-initialize
		for (i = 0; i < num_int_circs; ++i) {
			if (temp_sol[i] == DUMMY_VAR) { //if the distance is unknown, see if we have enough information to calculate it
				unknown_var = 1;
				num_common_points = 0; //initialize
				point1 = idx2row(i,max_col); //1 based (for both row and col here)
				point2 = idx2col(i,point1,max_col);
				for (j = 1; j <= max_col; ++j) { //we let point1 and point2 = rows and then we search through all cols not equal to point1 or point2
					if (j != point1 && j != point2) {
						if (point1 < j)
							idx1 = row_col_2_idx(point1,j,max_col);
						else
							idx1 = row_col_2_idx(j,point1,max_col);
						if (point2 < j)
							idx2 = row_col_2_idx(point2,j,max_col);
						else
							idx2 = row_col_2_idx(j,point2,max_col);
						if (temp_sol[idx1] != DUMMY_VAR && temp_sol[idx2] != DUMMY_VAR) { //if the distances between point1 and point2 and this new point are known, then proceed
							common_points[num_common_points] = j; //1 based
							++num_common_points;
						}
					}
				} //end of searching through all cols to find which points have distances which are known to both point1 and point2
				if (num_common_points >= 3) { //only bother to try and solve the relative distance if we have at least 3 points whose distances are known to both point1 and point2 (as polygons must share a common base)
					//printf("\n\nSOLVING FOR RELATIVE DISTANCE BETWEEN POINTS %d AND %d.\n\n", point1, point2);
					for (j = 0; j < num_common_points; ++j) { //loop through all possible combinations of 3 points whose distances are known to both point1 and point2
						for (k = j+1; k < num_common_points; ++k) {
							for (q = k+1; q < num_common_points; ++q) {
								found_base = 0; //initialize
								idx1 = row_col_2_idx(common_points[j],common_points[k],max_col); //for common base 123, all relative distances amongst the base must be known, thus here we find the index of r12
								idx2 = row_col_2_idx(common_points[j],common_points[q],max_col); //index associated with r13
								idx3 = row_col_2_idx(common_points[k],common_points[q],max_col); //index associated with r23
								if (temp_sol[idx1] != DUMMY_VAR && temp_sol[idx2] != DUMMY_VAR && temp_sol[idx3] != DUMMY_VAR) { //if all of these distances are known, then we have enough information to calculate rij
									for (z = 0; z < *mult_num_sol_count; ++z) { //search through all possible partial solution vectors
										//printf("mult_num_sol_count = %d, z = %d\n",*mult_num_sol_count, z);
										//for (p = 0; p < num_int_circs; ++p)
											//printf("temp_sol = %f\n",temp_sol[p]);
										num_cons_sols = 0; //initialize
										//printf("\nr_%d_%d, k,p,q = %d,%d,%d\n",point1,point2,common_points[j],common_points[k],common_points[q]);
										if (temp == 0) {
											num_cons_sols = calc_rel_dist_num(mult_sol_num_temp, num_sol, z, i, &num_new_sols, max_col, num_int_circs, point1, point2, common_points[j], common_points[k], common_points[q]);
											//printf("temp = %d, write to temp = mn_sol1\n",temp);
											//for (p = 0; p < num_int_circs; ++p)
												//printf("mn_sol = %f\n",mult_sol_num_temp[z][p]);
										}
										else {
											num_cons_sols = calc_rel_dist_num(num_sol, mult_sol_num_temp, z, i, &num_new_sols, max_col, num_int_circs, point1, point2, common_points[j], common_points[k], common_points[q]);
											//printf("temp = %d, write to num_sol = mn_sol2\n",temp);
											//for (p = 0; p < num_int_circs; ++p)
												//printf("mn_sol = %f\n",num_sol[z][p]);
										}
										//printf("num consistent sols = %d\n",num_cons_sols);
										//printf("\ntemp_sol[%d] = %f; found_base = %d",i,temp_sol[i],found_base);
										if (num_cons_sols != 0) {
											found_base = 1;
											//printf("\ntemp_solved_dists = %d\n",temp_solved_dists);
										} 
									} //end of looping through all potential solution columns
									if (found_base == 1) { //perhaps combine this with one of the later found_base statments, if are going to eliminate violated from them
										++num_solved_dists;
										*mult_num_sol_count = num_new_sols; //if it has enough info to calc a rel dist, but can't find a consistent solution, then there is no solution - but found_base = 1 implies at least one solution was found
										num_new_sols = 0; //re-initialize
										/*for (z = 0; z < *mult_num_sol_count; ++z) {
											printf("\n");
											for (p = 0; p < num_int_circs; ++p)
												printf("mn_sol1 = %f\n",mult_sol_num_temp[z][p]);
										}
										for (z = 0; z < *mult_num_sol_count; ++z) {
											printf("\n");
											for (p = 0; p < num_int_circs; ++p)
												printf("mn_sol2 = %f\n",num_sol[z][p]);
										}*/
										if (temp == 0) { //if we've solved a distance, then we'll have to document the solutions in the new solution array, and we should now use that one to feed in solutions, so switch temp values to indicate this -- also, refill temp_sol vector with first column of the most recently documented solution array
											for (z = 0; z < num_int_circs; ++z)
												temp_sol[z] = mult_sol_num_temp[0][z];
											temp = 1;
											/*for (z = 0; z < *mult_num_sol_count; ++z) {
												printf("\n");
												for (p = 0; p < num_int_circs; ++p)
													printf("mn_sol = %f\n",mult_sol_num_temp[z][p]);
											}*/
										}
										else {
											for (z = 0; z < num_int_circs; ++z)
												temp_sol[z] = num_sol[0][z];
											temp = 0;
											/*for (z = 0; z < *mult_num_sol_count; ++z) {
												printf("\n");
												for (p = 0; p < num_int_circs; ++p)
													printf("mn_sol = %f\n",num_sol[z][p]);
											}*/
										}
									}
									else { //if there was enough information to calculate the new relative distance, and no consistent relative distance for any of the potential solutions was found, then there is no solution
										//printf("\nThis violated was called!\n");
										violated = 1;
									}
								}
								if (found_base == 1 || violated == 1)
									break;
							}
							if (found_base == 1 || violated == 1)
								break;
						}
						if (found_base == 1 || violated == 1)
							break;
					} //end of looping through all possible combinations of 3 points
				}
			}
			if (violated == 1)
				break;
		} //end of looping through all relative distances
		if (num_solved_dists == 0 || violated == 1) //if we haven't solved any new distances, stop searching as either all distances have been solved, or we don't have enough information to solve any new distances
			stop = 1;
	} //end of calculating new relative distances
	
	
	//printf("temp2 = %d",temp);
	if (temp != 0) { //set the solution storage to the num_sol array, regardless of where it ended it up
		//printf("\nWrite to mn_sol1 = temp (temp != 0)\n");
		for (i = 0; i < *mult_num_sol_count; ++i) {
			printf("\n");
			for (j = 0; j < num_int_circs; ++j) {
				num_sol[i][j] = mult_sol_num_temp[i][j];
				//printf("num_sol = %f\n",num_sol[i][j]);
			}
		}
	}
	//else
		//printf("\nWrite to mn_sol2 (temp = 0)\n");
	
	//printf("\nmult_num_sol_count = %d\n",*mult_num_sol_count);
	
	//printf("\nNum valid solutions1 = %d\n",*valid_sols_count);
	
	/*printf("\n%d Solutions so far...before consistency check:\n", *mult_num_sol_count);
	for (i = 0; i < *mult_num_sol_count; ++i) {
		printf("\nSol %d:\n",i+1);
		for (j = 0; j < num_int_circs; ++j) {
			printf("num_sol[%d][%d] = %f\n",idx2row(j,max_col),idx2col(j,idx2row(j,max_col),max_col),num_sol[i][j]);
		}
	}*/
	
	//printf("\nBefore 1st global check, violated = %d, unknown_var = %d.\n",violated,unknown_var);
	/*Check that the global solution found is consistent - HAVE TO DO THIS FOR ALL POTENTIAL SOLUTIONS FOUND - SO WRITE THE EXTRA FOR LOOP FOR THIS*/
	if ((violated == 0) && (unknown_var == 1)) { //the unknown_var = 1 part only checks for global consistency if there was an unknown relative distance to look for
		//printf("\nFirst Global Check was Called...\n");
		for (z = 0; z < *mult_num_sol_count; ++z) { //loop through all potential solutions
			consistent = 1; //initialize
			for (k = 0; k < num_int_circs; ++k)
				temp_sol[k] = num_sol[z][k]; //write the temporary solution vector to test
			for (i = 1; i <= max_col; ++i) { //must enter in 1 based row and column indices
				for (j = i+1; j <= max_col; ++j) {
					idx_check = row_col_2_idx(i,j,max_col);
					if (temp_sol[idx_check] != DUMMY_VAR) {
						//printf("\nTHE NEW CHECK!!\n");
						consistent = consistency_check(temp_sol, max_col, i, j, temp_sol[idx_check]);
					}
					if (consistent == 0) {
						//printf("\nSol %d was found globally inconsistent by consistency check\n",z+1);
						violated = 1;
						break;
					}
				}
				if (violated == 1)
					break;
			}
			//printf("\nsol_col = %d, violated = %d\n",z,violated);
			if (violated == 0) { //if the solution is consistent, document it
				for (i = 0; i < num_int_circs; ++i)
					num_sol[*valid_sols_count][i] = num_sol[z][i];
				++*valid_sols_count;
			}
			violated = 0; //re-initialize
		} //end of looping through all potential solutions
	}
	
	//printf("\nNum valid solutions2 = %d\n",*valid_sols_count);
	
	/*If relative distances have been calculated, check for global consistency.*/
	if ((found_base == 1) && (violated == 0) && (*valid_sols_count > 0) && (unknown_var == 1)) { //the unknown_var = 1 part only checks for global consistency if there was an unknown relative distance to look for
		for (i = 0; i < *valid_sols_count; ++i) {
			for (j = 0; j < num_int_circs; ++j)
				temp_sol[j] = num_sol[i][j];
			violated = global_consistency_check(temp_sol, max_col, only_print_sols, violations);
			//printf("\nvalid_col_count = %d, violated = %d\n",i,violated);
			if (violated == 0) { //if the solution is consistent, document it
				for (k = 0; k < num_int_circs; ++k)
					num_sol[temp_count][k] = num_sol[i][k];
				++temp_count;
			}
			violated = 0; //re-initialize
		} //end of looping through all potential solutions
		*valid_sols_count = temp_count;
	}
	
	if (unknown_var == 0) //if all relative distances were already found, then document that there exists one solution
		*valid_sols_count = 1;
	
	if ((found_base == 1) && (*valid_sols_count == 0))
		violated = 1;
		
	if (violated == 1) {
		if (only_print_sols == 0)
			printf("\nGraph is unphysical because it has no solution or it implies that a distance is < R.\n");
		++violations[8];
	}
	
	//printf("Num valid solutions3 = %d\n",*valid_sols_count);
	
	return violated;
	
}

/*
Description: If new relative distances have been found - then check for global consistency (I may also write in a check later to check for global consistency of matrices before or if new distances have not been calculated).

Output: Returns 1 if a violation has occured, and 0 if a violation has not occured.
*/
int global_consistency_check(double sol_num[], const int max_col, const int only_print_sols, int *violations) {

	/*Declare and initialize local variables*/
	int i, j, k, q, p, z, s, l, idx_ij, idx_ik, idx_ip, idx_iq, idx_jk, idx_jp, idx_jq, idx_pk, idx_pq, idx_kq, zero_consistency, zero_check, nums_A1, nums_A2, nums_A3, violated, round_num1, round_num2, round_num3, round_num4, normal_check, num_points_on_line, num_dists_geq4, num_dists_geq3, num_dists_geq2;
	int points_on_line[4];
	double a1, b1, c1, a2, b2, c2, a3, b3, c3, num_A1, num_A2, num_A3, arg_A1, arg_A2, arg_A3, A1, A2, A3, PI;
	double A1_array[4], A2_array[4], A3_array[4];
	violated = 0;
	round_num1 = 5; //original runs have this set to 5 - this is precision to which dihedral angle = 0 or Pi is checked
	round_num2 = 5; //original runs have this set to 5 - this is precision to which A1 = sum or difference of A2 and A3 is checked
	round_num3 = 5; //original runs have this set to 5 - this is precision to which both numerator and denominator = 0, for the zero canceling case, is checked
	round_num4 = 5; //if don't use a rounded number for the in line check, then graph 11633 for n = 9 is erroneously ruled out
	//PI = 3.141592653589793238462643;
	//Pi to 1024 decimal places:
	PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327887;
	
	for (i = 1; i <= max_col; ++i) {
		for (j = i+1; j <= max_col; ++j) {
			for (k = j+1; k <= max_col; ++k) {
				for (p = k+1; p <= max_col; ++p) {
					for (q = p+1; q <= max_col; ++q) {
						//printf("\ni = %d, j = %d, k = %d, p = %d, q = %d\n",i,j,k,p,q);
						idx_ij = row_col_2_idx(i,j,max_col);
						idx_ik = row_col_2_idx(i,k,max_col);
						idx_ip = row_col_2_idx(i,p,max_col);
						idx_iq = row_col_2_idx(i,q,max_col);
						idx_jk = row_col_2_idx(j,k,max_col);
						idx_jp = row_col_2_idx(j,p,max_col);
						idx_jq = row_col_2_idx(j,q,max_col);
						idx_pk = row_col_2_idx(k,p,max_col);
						idx_pq = row_col_2_idx(p,q,max_col);
						idx_kq = row_col_2_idx(k,q,max_col);
						a1 = law_of_cosines_num(sol_num[idx_ip],sol_num[idx_jp],sol_num[idx_ij],3);
						b1 = law_of_cosines_num(sol_num[idx_jp],sol_num[idx_jk],sol_num[idx_pk],2);
						c1 = law_of_cosines_num(sol_num[idx_ip],sol_num[idx_ik],sol_num[idx_pk],2);
						a2 = law_of_cosines_num(sol_num[idx_ip],sol_num[idx_iq],sol_num[idx_pq],2);
						b2 = law_of_cosines_num(sol_num[idx_kq],sol_num[idx_pk],sol_num[idx_pq],1);
						c2 = law_of_cosines_num(sol_num[idx_ip],sol_num[idx_ik],sol_num[idx_pk],2);
						a3 = law_of_cosines_num(sol_num[idx_jp],sol_num[idx_jq],sol_num[idx_pq],2);
						b3 = law_of_cosines_num(sol_num[idx_jk],sol_num[idx_jp],sol_num[idx_pk],1);
						c3 = law_of_cosines_num(sol_num[idx_kq],sol_num[idx_pk],sol_num[idx_pq],1);
						//printf("\na1 = %f, b1 = %f, c1 = %f, arg_A1 = %f\n",a1,b1,c1,arg_A1);
						//printf("\na2 = %f, b2 = %f, c2 = %f, arg_A2 = %f\n",a2,b2,c2,arg_A2);
						//printf("\na3 = %f, b3 = %f, c3 = %f, arg_A3 = %f\n",a3,b3,c3,arg_A3);
						//printf("\nr_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f\n",i,j,sol_num[idx_ij],i,k,sol_num[idx_ik],i,p,sol_num[idx_ip],i,q,sol_num[idx_iq],j,k,sol_num[idx_jk],j,p,sol_num[idx_jp],j,q,sol_num[idx_jq],p,k,sol_num[idx_pk],p,q,sol_num[idx_pq],k,q,sol_num[idx_kq]);
						//printf("\n1!");
						normal_check = 0; //initialize
						if ((round(a1*pow(10,round_num1)) >= 1) && (round(b1*pow(10,round_num1)) >= 1) && (round(c1*pow(10,round_num1)) >= 1) && (round(a2*pow(10,round_num1)) >= 1) && (round(b2*pow(10,round_num1)) >= 1) && (round(c2*pow(10,round_num1)) >= 1) && (round(a3*pow(10,round_num1)) >= 1) && (round(b3*pow(10,round_num1)) >= 1) && (round(c3*pow(10,round_num1)) >= 1)) { //check that none of the angles = 0 to round_num1 decimal places, and only perform dihedral angle check if this is satisfied
							if ((fabs(round(a1*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(b1*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(c1*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(a2*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(b2*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(c2*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(a3*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(b3*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1) && (fabs(round(c3*pow(10,round_num1)) - round(PI*pow(10,round_num1))) >= 1)) { //check that none of the angles = pi to round_num1 decimal places, and only perform dihedral angle check if this is satisfied
								normal_check = 1;
								//printf("2\n");
								num_A1 = cos(a1) - cos(b1)*cos(c1);
								num_A2 = cos(a2) - cos(b2)*cos(c2);
								num_A3 = cos(a3) - cos(b3)*cos(c3);
								arg_A1 = num_A1/(sin(c1)*sin(b1));
								arg_A2 = num_A2/(sin(c2)*sin(b2));
								arg_A3 = num_A3/(sin(c3)*sin(b3));
								arg_A1 = check_round(arg_A1,num_A1);
								arg_A2 = check_round(arg_A2,num_A2);
								arg_A3 = check_round(arg_A3,num_A3);
								A1 = acos(arg_A1);
								A2 = acos(arg_A2);
								A3 = acos(arg_A3);
								//printf("\na1 = %f, b1 = %f, c1 = %f, arg_A1 = %f\n",a1,b1,c1,arg_A1);
								//printf("\na2 = %f, b2 = %f, c2 = %f, arg_A2 = %f\n",a2,b2,c2,arg_A2);
								//printf("\na3 = %f, b3 = %f, c3 = %f, arg_A3 = %f\n",a3,b3,c3,arg_A3);
								//printf("\nAHHHH:\n A1 = %f, A2 = %f, A3 = %f, A2 + A3 = %f, A3 - A2 = %f, A2 - A3 = %f, 2pi - (A1 + A2) = %f, 2pi - (A3 - A2) = %f, 2pi - (A2 - A3) = %f\n",A1,A2,A3, A2+A3, A3-A2, A2-A3, 2*PI-(A2+A3),2*PI-(A3-A2),2*PI-(A2-A3));
								//printf("\n%f = %f\n",round(A1*pow(10,4)),round((A2 + A3)*pow(10,4)));
								//printf("\nr_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f, r_%d%d = %f\n",i,j,sol_num[idx_ij],i,k,sol_num[idx_ik],i,p,sol_num[idx_ip],i,q,sol_num[idx_iq],j,k,sol_num[idx_jk],j,p,sol_num[idx_jp],j,q,sol_num[idx_jq],p,k,sol_num[idx_pk],p,q,sol_num[idx_pq],k,q,sol_num[idx_kq]);
								//higher precision
								//printf("\nr_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f, r_%d%d = %6.12f\n",i,j,sol_num[idx_ij],i,k,sol_num[idx_ik],i,p,sol_num[idx_ip],i,q,sol_num[idx_iq],j,k,sol_num[idx_jk],j,p,sol_num[idx_jp],j,q,sol_num[idx_jq],p,k,sol_num[idx_pk],p,q,sol_num[idx_pq],k,q,sol_num[idx_kq]);
								//if ((round(A1*pow(10,4)) != round((A2 + A3)*pow(10,4))) && (round(A1*pow(10,4)) != round((A2 - A3)*pow(10,4))) && (round(A1*pow(10,4)) != round((A3 - A2)*pow(10,4))) && (round(A1*pow(10,4)) != round((2*PI-(A2 + A3))*pow(10,4))) && (round(A1*pow(10,4)) != round((2*PI-(A2 - A3))*pow(10,4))) && (round(A1*pow(10,4)) != round((2*PI-(A3 - A2))*pow(10,4)))) {//if A1 is not the sum or difference of the other 2 dihedral angles, then there is an inconsistency
								//printf("\nI'm here and I didn't go through the check because dihedral angle is not well defined...\n");
								if ((fabs(round(A1*pow(10,round_num2)) - round((A2 + A3)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((A2 - A3)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((A3 - A2)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A2 + A3))*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A2 - A3))*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A3 - A2))*pow(10,round_num2))) > 1)) {//if A1 is not the sum or difference of the other 2 dihedral angles, then there is an inconsistency (this is done to a precision of round_num2)
									//printf("\nNOPE, I DID!!!\n");
									zero_consistency = 1; //initialize
									zero_check = 0; //initialize
									A1_array[0] = A1; //initialize
									nums_A1 = 1;
									A2_array[0] = A2;
									nums_A2 = 1;
									A3_array[0] = A3;
									nums_A3 = 1;
									if (((fabs(round(sin(b1)*pow(10,round_num3))) - 0 <= 1) || (fabs(round(sin(c1)*pow(10,round_num3))) - 0 <= 1)) && (fabs(round((cos(a1) - cos(b1)*cos(c1))*pow(10,round_num3))) - 0 <= 1)) { //if terms in denominator and numerator = 0, check for cancelling of zeros (check this for each tetrahedron)
										zero_check = 1;
										//printf("\nA1 recalc was called!\n");
										//A1 = recalc_dihed_ang_0cancelling(const double r_bc, const double r_ab, const double r_ac, const double r_b, const double r_a, const double r_c);
										//A1 = recalc_dihed_ang_0cancelling(sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_ip],sol_num[idx_jk],sol_num[idx_ij],sol_num[idx_ik]); //recalculate the dihedral angles with the cancelled out zeros...
										recalc_dihed_ang_0cancelling(A1_array,&nums_A1,sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_ip],sol_num[idx_jk],sol_num[idx_ij],sol_num[idx_ik]); //recalculate the dihedral angles with the cancelled out zeros...
									}
									if (((fabs(round(sin(b2)*pow(10,round_num3))) - 0 <= 1) || (fabs(round(sin(c2)*pow(10,round_num3))) - 0 <= 1)) && (fabs(round((cos(a2) - cos(b2)*cos(c2))*pow(10,round_num3))) - 0 <= 1)) {
										zero_check = 1;
										//printf("\nA2 recalc was called!\n");
										//A2 = recalc_dihed_ang_0cancelling(sol_num[idx_pk],sol_num[idx_pq],sol_num[idx_ip],sol_num[idx_kq],sol_num[idx_iq],sol_num[idx_ik]);
										recalc_dihed_ang_0cancelling(A2_array,&nums_A2,sol_num[idx_pk],sol_num[idx_pq],sol_num[idx_ip],sol_num[idx_kq],sol_num[idx_iq],sol_num[idx_ik]);
									}
									if (((fabs(round(sin(b3)*pow(10,round_num3))) - 0 <= 1) || (fabs(round(sin(c3)*pow(10,round_num3))) - 0 <= 1)) && (fabs(round((cos(a3) - cos(b3)*cos(c3))*pow(10,round_num3))) - 0 <= 1)) {
										zero_check = 1;
										//printf("\nA3 recalc was called!\n");
										//A3 = recalc_dihed_ang_0cancelling(sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_pq],sol_num[idx_jk],sol_num[idx_jq],sol_num[idx_kq]);
										recalc_dihed_ang_0cancelling(A3_array,&nums_A3,sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_pq],sol_num[idx_jk],sol_num[idx_jq],sol_num[idx_kq]);
										//printf("\nDists fed in to calc A3: r_pk = %f, r_jp = %f, r_pq = %f, r_jk = %f, r_jq = %f, r_kq = %f.\n",sol_num[idx_pk],sol_num[idx_jp],sol_num[idx_pq],sol_num[idx_jk],sol_num[idx_jq],sol_num[idx_kq]);
									}
									if (zero_check == 1) { //recheck global consistency after zero cancelling
										//printf("\nnums_A1 = %d, nums_A2 = %d, nums_A3 = %d\n",nums_A1,nums_A2,nums_A3);
										for (z = 0; z < nums_A1; ++z) {
											A1 = A1_array[z];
											//printf("\nA1 = %f from A1 value #%d, arg_A1 = %f\n",A1,z+1,cos(A1));
											for (s = 0; s < nums_A2; ++s) {
												A2 = A2_array[s];
												//printf("\nA2 = %f from A2 value #%d, arg_A2 = %f\n",A2,s+1,cos(A2));
												for (l = 0; l < nums_A3; ++l) {
													A3 = A3_array[l];
													//printf("\nA3 = %f from A3 value #%d, arg_A3 = %f\n",A3,l+1,cos(A3));
													zero_consistency = 1; //re-initialize before checking
													//printf("\nFrom Recalc After Zero Cancelling:\n A1 = %f, A2 = %f, A3 = %f, A2 + A3 = %f, A3 - A2 = %f, A2 - A3 = %f, 2pi - (A1 + A2) = %f, 2pi - (A3 - A2) = %f, 2pi - (A2 - A3) = %f\n",A1,A2,A3, A2+A3, A3-A2, A2-A3, 2*PI-(A2+A3),2*PI-(A3-A2),2*PI-(A2-A3));
													if ((fabs(round(A1*pow(10,round_num2)) - round((A2 + A3)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((A2 - A3)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((A3 - A2)*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A2 + A3))*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A2 - A3))*pow(10,round_num2))) > 1) && (fabs(round(A1*pow(10,round_num2)) - round((2*PI-(A3 - A2))*pow(10,round_num2))) > 1)) {//if A1 is not the sum or difference of the other 2 dihedral angles, then there is an inconsistency
														zero_consistency = 0;
													}
													else { //if we find a set of consistent dihedral angles, then the solution is consistent, and break this search
														break;
													}
												}
												if (zero_consistency == 1) //break if found a consistent set of dihedral angles
													break;
											}
											if (zero_consistency == 1) //break if found a consistent set of dihedral angles
												break;
										}
									}
									else { //if no zero cancelling occurs, then the solution is inconsistent
										zero_consistency = 0;
									}
									if (zero_consistency == 0) { //if it's still inconsistent after zero cancelling...
										//printf("\nGlobal Inconsistency!\n");
										violated = 1;
										break;
									}
								}
							} //end of consistency check performed only if none of the angles (ai,bi,ci) = pi, i.e. performed only if the dihedral angle is well defined
						} //end of consistency check performed only if none of the angles (ai,bi,ci) = 0, i.e. performed only if the dihedral angle is well defined
						//normal_check = 1; //set it to this if don't want to do the below check
						if (normal_check == 0) { //if the dihedral angle was not well defined because 3 or more points lie in a line, then perform the alternative consistency check, for points in a line.  Note that I started to write this for p points in a line, but that the below is only appropriately establishing when 3 points lie in a line, not when >3 points lie in a line, and thus if want to include the >3 points in a line check, then have to first write the appropriate check that establishes > 3 points lie in a line.  Write now it establishes various sets of 3 points totaling either 3, 4, or 5 points lie in a line - but it does not establish whether the angle(s) between the various combinations of 3/4 or 3/5 points are also 0 or pi - for example can have 2 overlapping sets of 3 points that total 5 unique points that are hinged around 1 point - i.e. 2 unique points in each set, and one overlapping point that forms the hinge for which the angle is not 0 or pi, and thus there exists 2 seperate sets of 3 points that lie on 2 seperate lines, and thus define a plane, not a line.  So, the 4 and 5 point checks are commented out below and should remain commented out until the appropriate check is inserted.  Also note that this check does not seem to be necessary for up to n = 10, the 3 point check seems to be sufficient.
							for (z = 0; z < 5; ++z) //initialize the points_on_line vector
								points_on_line[z] = 0;
							points_on_line[3] = 1; //this is point p and is involved in every angle, and so is always = 1
							if ((round(a1*pow(10,round_num1)) == 0) || (fabs(round((a1-PI)*pow(10,round_num1))) == 0)) {
								points_on_line[0] = 1; //point i
								points_on_line[1] = 1; //point j
							}
							if ((round(b1*pow(10,round_num1)) == 0) || (fabs(round((b1-PI)*pow(10,round_num1))) == 0) || (round(b3*pow(10,round_num1)) == 0) || (fabs(round((b3-PI)*pow(10,round_num1))) == 0)) { //b1 and b3 are the same angle, so either = 0 or PI gives the same result
								points_on_line[2] = 1; //point k
								points_on_line[1] = 1; //point j
							}
							if ((round(c1*pow(10,round_num1)) == 0) || (fabs(round((c1-PI)*pow(10,round_num1))) == 0) || (round(c2*pow(10,round_num1)) == 0) || (fabs(round((c2-PI)*pow(10,round_num1))) == 0)) { //c1 and c2 are the same angle, so either = 0 or PI gives the same result
								points_on_line[0] = 1; //point i
								points_on_line[2] = 1; //point k
							}
							if ((round(a2*pow(10,round_num1)) == 0) || (fabs(round((a2-PI)*pow(10,round_num1))) == 0)) {
								points_on_line[0] = 1; //point i
								points_on_line[4] = 1; //point q
							}
							if ((round(b2*pow(10,round_num1)) == 0) || (fabs(round((b2-PI)*pow(10,round_num1))) == 0) || (round(c3*pow(10,round_num1)) == 0) || (fabs(round((c3-PI)*pow(10,round_num1))) == 0)) {
								points_on_line[2] = 1; //point k
								points_on_line[4] = 1; //point q
							}
							if ((round(a3*pow(10,round_num1)) == 0) || (fabs(round((a3-PI)*pow(10,round_num1))) == 0)) {
								points_on_line[1] = 1; //point j
								points_on_line[4] = 1; //point q
							}
							num_points_on_line = 0; //initialize
							for (z = 0; z < 5; ++z) {
								if (points_on_line[z] == 1)
									++num_points_on_line;
							}
							/*printf("\nNum points on line = %d:",num_points_on_line);
							for (z = 0; z < 5; ++z)
							    printf("%d",points_on_line[z]);
							printf("\n");*/
							/*if (num_points_on_line == 5) {
								num_dists_geq4 = 0; //initialize
								num_dists_geq3 = 0;
								num_dists_geq2 = 0;
								if ((round((sol_num[idx_ij] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ik] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_iq] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jk] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jq] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 4)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_kq] - 4)*pow(10,round_num4)) >= 0)) // check that 1 distance is >= 4
									++num_dists_geq4;
								if (num_dists_geq4 < 1)
									violated = 1;
								printf("\nnum_dists_geq4 = %d\n",num_dists_geq4);
								if (violated == 0) { //only bother to continue to check if hasn't already been violated
									if (round((sol_num[idx_ij] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_ik] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_ip] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_iq] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_jk] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_jp] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_jq] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_pk] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_pq] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (round((sol_num[idx_kq] - 3)*pow(10,round_num4)) >= 0)
										++num_dists_geq3;
									if (num_dists_geq3 < 2)
										violated = 1;
									printf("\nnum_dists_geq3 = %d\n",num_dists_geq3);
								}
								if (violated == 0) { //only bother to continue to check if hasn't already been violated
									if (round((sol_num[idx_ij] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_ik] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_iq] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_jk] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_jq] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (round((sol_num[idx_kq] - 2)*pow(10,round_num4)) >= 0)
										++num_dists_geq2;
									if (num_dists_geq2 < 3)
										violated = 1;
									printf("\nnum_dists_geq2 = %d\n",num_dists_geq2);
								}
							}
							else if (num_points_on_line == 4) {
								num_dists_geq3 = 0; //initialize
								num_dists_geq2 = 0;
								if (points_on_line[0] == 0) { //point i is not on the line
									if ((round((sol_num[idx_jk] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_kq] - 3)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 3
										++num_dists_geq3;
									if (num_dists_geq3 < 1)
										violated = 1;
									if (violated == 0) { //only bother to continue checking if haven't already violated consistency
										if (round((sol_num[idx_jk] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_kq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (num_dists_geq2 < 2)
											violated = 1;
									}
								}
								if (points_on_line[1] == 0) { //point j is not on the line
									if ((round((sol_num[idx_ik] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_iq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_kq] - 3)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 3
										++num_dists_geq3;
									if (num_dists_geq3 < 1)
										violated = 1;
									if (violated == 0) { //only bother to continue checking if haven't already violated consistency
										if (round((sol_num[idx_ik] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_iq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_kq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (num_dists_geq2 < 2)
											violated = 1;
									}
								}
								if (points_on_line[2] == 0) { //point k is not on the line
									if ((round((sol_num[idx_ij] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_iq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jq] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 3)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 3
										++num_dists_geq3;
									if (num_dists_geq3 < 1)
										violated = 1;
									if (violated == 0) { //only bother to continue checking if haven't already violated consistency
										if (round((sol_num[idx_ij] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_iq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (num_dists_geq2 < 2)
											violated = 1;
									}
								}
								if (points_on_line[4] == 0) { //point q is not on the line (no need to check point p, as if we're in this situation it is always on the line)
									if ((round((sol_num[idx_ij] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ik] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jk] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 3)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 3)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 3
										++num_dists_geq3;
									if (num_dists_geq3 < 1)
										violated = 1;
									if (violated == 0) { //only bother to continue checking if haven't already violated consistency
										if (round((sol_num[idx_ij] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_ik] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jk] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)
											++num_dists_geq2;
										if (num_dists_geq2 < 2)
											violated = 1;
									}
								}
							}
							else if (num_points_on_line == 3) {*/
							if (num_points_on_line == 3) {
								num_dists_geq2 = 0; //initialize
								if ((points_on_line[0] == 0) && (points_on_line[1] == 0)) { //points i and j are not on the line
									if ((round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_kq] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2 to precision of round_num4
										++num_dists_geq2;
								}
								if ((points_on_line[0] == 0) && (points_on_line[2] == 0)) { //points i and k are not on the line
									if ((round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jq] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2
										++num_dists_geq2;
								}
								if ((points_on_line[0] == 0) && (points_on_line[4] == 0)) { //points i and q are not on the line
									if ((round((sol_num[idx_jk] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2
										++num_dists_geq2;
								}
								if ((points_on_line[1] == 0) && (points_on_line[2] == 0)) { //points j and k are not on the line
									if ((round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_iq] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pq] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2
										++num_dists_geq2;
								}
								if ((points_on_line[1] == 0) && (points_on_line[4] == 0)) { //points j and q are not on the line
									if ((round((sol_num[idx_ik] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_pk] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2
										++num_dists_geq2;
								}
								if ((points_on_line[2] == 0) && (points_on_line[4] == 0)) { //points k and q are not on the line
									if ((round((sol_num[idx_ij] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_ip] - 2)*pow(10,round_num4)) >= 0) || (round((sol_num[idx_jp] - 2)*pow(10,round_num4)) >= 0)) //check that at least one of the relevant distances is >= 2
										++num_dists_geq2;
								}
								if (num_dists_geq2 < 1)
									violated = 1;
							}
							/*if (violated == 0)
								printf("\nConsistency check in the line was performed and nothing was violated\n");
							else if (violated == 1)
								printf("\nConsistency check in the line was performed and WAS violated\n");*/
						} //end of check for 3 or more points lying on a line
					}
					if (violated == 1)
						break;
				}
				if (violated == 1)
					break;
			}
			if (violated == 1)
				break;
		}
		if (violated == 1)
			break;
	}

	return violated;
	
}