// This is for the parallel processing

// Declare start and end flags
int flag1_start = 0, flag2_start = 0, flag3_start = 0, flag4_start = 0;
int flag1_end = 0, flag2_end = 0, flag3_end = 0, flag4_end = 0;

double V_grad_max1 = 0, V_grad_max2 = 0, V_grad_max3 = 0, V_grad_max4 = 0, V_grad_max_all = 0;

double E1 = 0, E2 = 0, E3 = 0, E4 = 0;
double VG1 = 0, VG2 = 0, VG3 = 0, VG4 = 0;
double E_grad_dir = 0;

// Declare iostream files for flags and intermediate variables
ifstream file_flag_read; ofstream file_flag_write;
ofstream file_par_write; ifstream file_par_read;

// Set all start flags to 0
#include "write_start_flags.h"

// Set all end flags to 0
#include "write_end_flags.h"


//******************************************************************************************************
// Here starts the parallel loop

int Change_A_flag = 1;
int Exit_flag = 0;
file_par_write.open("Exit_flag.m", ios::out);
file_par_write << Exit_flag << endl;
file_par_write.close();

const double A_change_amount = .01;
double Current_con_factor = 1.0;

int t_count = 0;




t = 0; t_count = 0;
while (1) {	
	
	
	if (Exit_flag == 1 && Change_A_flag == 1) {
		file_par_write.open("Exit_flag.m", ios::out);
		file_par_write << Exit_flag << endl;
		file_par_write.close();
		break;
	}
	
	// Isotropically decrease A0
	if (Change_A_flag == 1){
		if ((Current_con_factor - A_change_amount) >= A_con_factor){
			Current_con_factor -= A_change_amount;
		}
		else {
			Exit_flag = 1;
		}
		
		cout << endl <<  "The Apical cells are now constricting to " << Current_con_factor << " of their original lengths" << endl;
		Change_A_flag = 0;
		file_par_write.open("con_factor.m", ios::out);
		file_par_write << Current_con_factor << endl;
		file_par_write.close();
	}
		
		//******************************************************************************************************
		// Start slave proceduce 1: Export the corresponding intermediate verticies to other slaves
		// update Current_con_factor
		cout << "Starting Slave Proceduces 1 for t = " << t << endl;
		flag1_end = 0; flag2_end = 0; flag3_end = 0; flag4_end = 0;
		#include "write_end_flags.h"
		flag1_start = 1; flag2_start = 1; flag3_start = 1; flag4_start = 1;
		#include "write_start_flags.h"
		cout << "Slave Procedure Started" << endl; 
		cout << "Waiting for procedure to end" << endl;
		#include "wait_for_end_flags.h"
		
		//******************************************************************************************************
		// Start slave proceduce 2: Import the corresponding intermediate verticies from other slaves
		// Calculate the gradient and export the corresponding imtermediate gradient to other slaves
		cout << "Starting Slave Proceduces 2 for t = " << t << endl;
		flag1_end = 0; flag2_end = 0; flag3_end = 0; flag4_end = 0;
		#include "write_end_flags.h"
		flag1_start = 1; flag2_start = 1; flag3_start = 1; flag4_start = 1;
		#include "write_start_flags.h"
		cout << "Slave Procedure Started" << endl; 
		cout << "Waiting for procedure to end" << endl;
		#include "wait_for_end_flags.h"
		
		//******************************************************************************************************
		// Start slave proceduce 3: Import the corresponding intermediate gradient from other slaves
		// Export the maximum mag to the master program
		// Master program determines the overall max and exports it back to the slaves
		cout << "Starting Slave Proceduces 3 for t = " << t << endl;
		flag1_end = 0; flag2_end = 0; flag3_end = 0; flag4_end = 0;
		#include "write_end_flags.h"
		flag1_start = 1; flag2_start = 1; flag3_start = 1; flag4_start = 1;
		#include "write_start_flags.h"
		cout << "Slave Procedure Started" << endl; 
		cout << "Waiting for procedure to end" << endl;
		#include "wait_for_end_flags.h"
		#include "determine_grad_max.h"
		
		//******************************************************************************************************
		// Start slave proceduce 4: Rescale the gradient
		// Calculate E_old for each slave and export it back to the master program
		cout << "Starting Slave Proceduces 4 for t = " << t << endl;
		flag1_end = 0; flag2_end = 0; flag3_end = 0; flag4_end = 0;
		#include "write_end_flags.h"
		flag1_start = 1; flag2_start = 1; flag3_start = 1; flag4_start = 1;
		#include "write_start_flags.h"
		cout << "Slave Procedure Started" << endl; 
		cout << "Waiting for procedure to end" << endl;
		#include "wait_for_end_flags.h"
		cout << "Importing E_old" << endl; 
		
		E_old = 0;
		file_par_read.open("Energy_1.m", ios::in); file_par_read >> E1; file_par_read.close();
		file_par_read.open("Energy_2.m", ios::in); file_par_read >> E2; file_par_read.close();
		file_par_read.open("Energy_3.m", ios::in); file_par_read >> E3; file_par_read.close();
		file_par_read.open("Energy_4.m", ios::in); file_par_read >> E4; file_par_read.close();
		
		E_old = E1 + E2 + E3 + E4;
		//******************************************************************************************************
		// Start slave proceduce 5: Loop until (E_new - E_old <= E_grad_dir*c)
		// Calculate E_old for each slave and export it back to the master program
		cout << "Starting Final Slave Proceduce 5: the time loop for t = " << t << endl;
		while (1) {	
			flag1_end = 0; flag2_end = 0; flag3_end = 0; flag4_end = 0;
			#include "write_end_flags.h"
			flag1_start = 1; flag2_start = 1; flag3_start = 1; flag4_start = 1;
			#include "write_start_flags.h"
			
			#include "wait_for_end_flags.h"
			
			E_new = 0; E_grad_dir = 0;
			file_par_read.open("Energy_1.m", ios::in); file_par_read >> E1; file_par_read >> VG1; file_par_read.close();
			file_par_read.open("Energy_2.m", ios::in); file_par_read >> E2; file_par_read >> VG2; file_par_read.close();
			file_par_read.open("Energy_3.m", ios::in); file_par_read >> E3; file_par_read >> VG3; file_par_read.close();
			file_par_read.open("Energy_4.m", ios::in); file_par_read >> E4; file_par_read >> VG4; file_par_read.close();
			
			E_new = E1 + E2 + E3 + E4;
			E_grad_dir = VG1 + VG2 + VG3 + VG4;	
			cout << "E_old = " << E_old << " E_new = " << E_new << endl;
			if (E_new - E_old <= E_grad_dir) break;
			
		}
		
		E_total = E_new;
		//******************************************************************************************************
		// If a sufficient energy minimum is reached then set Change_A_flag to 1
		if ((fabs(E_total - E_total_old)/E_total < Percent_stop) || t_count == T_count_max) {
			Change_A_flag = 1;
			cout << "Energy Minimum Reached after = " << t_count << " interations" << endl;
			t_count = 0;
		}
		else {
			t_count++;
		}
		
		E_total_old = E_total;
		t++;
	
}		
		
return 0;
		
		
		
		
		
		
/*	
return 0;

	// Calculate total energy
	E_total = 0;
	for (i=0; i<Nc; i++) {
		// Get the apical verticies
		for (k=0; k<4; k++) {
			AP[k][0] = Vx[C_apver[i][k]]; AP[k][1] = Vy[C_apver[i][k]]; AP[k][2] = Vz[C_apver[i][k]];
			BS[k][0] = Vx[C_bsver[i][k]]; BS[k][1] = Vy[C_bsver[i][k]]; BS[k][2] = Vz[C_bsver[i][k]];
		}
		K_par[0] = A_K_cell[i]; K_par[1] = B_K_cell[i]; K_par[2] = V_K_cell[i];
		K_par[3] = L_K_cell[i]*C_side_factor[i][0]; K_par[4] = L_K_cell[i]*C_side_factor[i][1];
		K_par[5] = L_K_cell[i]*C_side_factor[i][2]; K_par[6] = L_K_cell[i]*C_side_factor[i][3];
		K_par[7] = D_K_cell[i];
		
		O_par[0] = A_0_cell[i]; O_par[1] = B_0_cell[i]; O_par[2] = V_0_cell[i];
		O_par[3] = L_0_cell[i]; O_par[4] = L_0_cell[i]; O_par[5] = L_0_cell[i];
		O_par[6] = L_0_cell[i]; O_par[7] = D_0_cell[i];
		E_total += E_from_cell_V2_diag(AP, BS, K_par, O_par);
	}
	
	// Energy from Plane constraint
	if (plane_constraint_present == 1){
		for (i=0; i<Nv; i++) {
			for (j=0; j<24; j++) {
				
				if ((Plane_test_p1[i][j] < 0) || (Plane_test_p2[i][j] < 0) || (Plane_test_p3[i][j] < 0))
					continue;
				
				Point1[0] = Vx[Plane_test_p1[i][j]];Point1[1] = Vy[Plane_test_p1[i][j]];Point1[2] = Vz[Plane_test_p1[i][j]];
				Point2[0] = Vx[Plane_test_p2[i][j]];Point2[1] = Vy[Plane_test_p2[i][j]];Point2[2] = Vz[Plane_test_p2[i][j]];
				Point3[0] = Vx[Plane_test_p3[i][j]];Point3[1] = Vy[Plane_test_p3[i][j]];Point3[2] = Vz[Plane_test_p3[i][j]];
				
				// Pertubation for x
				Point0[0] = Vx[i];Point0[1] = Vy[i];Point0[2] = Vz[i];				
				Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
				if (Dis_tmp < R_max) {
					E_total += Plane_potential[int(Dis_tmp/dx_plane)];
				}
			}
		}
	}
	
	// Energy from Yolk
	if (yolk_present == 1){
		Yolk_V = 0;
		for (i=0; i<Nc; i++) {
			
			for (k=0; k<4; k++) {
				BS[k][0] = Vx[C_bsver[i][k]]; BS[k][1] = Vy[C_bsver[i][k]]; BS[k][2] = Vz[C_bsver[i][k]];
				AP[k][0] = Yolk_x[C_ykver[i][k]]; AP[k][1] = Yolk_y[C_ykver[i][k]]; AP[k][2] = Yolk_z[C_ykver[i][k]];
			}
			
			Yolk_volumes[i] = Volume_Calc(AP, BS);
			Yolk_V += Yolk_volumes[i];
		}
		E_total += K_yolk*pow(Yolk_V - Yolk_O, 2);
	}
	
	
	if (t % T_display == 0)
		cout << "t = " << t << " E = " << E_total<< " E_diff = " << fabs(E_total - E_total_old) << " dx average = " << Average_Displacement << " count = " << t_count << endl;
		
	if ((fabs(E_total - E_total_old)/E_total < Percent_stop) || t_count == T_count_max) {
		Change_A_flag = 1;
		cout << "Energy Minimum Reached after = " << t_count << " interations" << endl;
		t_count = 0;
	}
	else {
		t_count++;
	}

		
		if (Change_A_flag == 1){
			
			for (i=0; i<Nc; i++) {
				
				for (k=0; k<4; k++) {
					AP[k][0] = Vx[C_apver[i][k]]; AP[k][1] = Vy[C_apver[i][k]]; AP[k][2] = Vz[C_apver[i][k]];
					BS[k][0] = Vx[C_bsver[i][k]]; BS[k][1] = Vy[C_bsver[i][k]]; BS[k][2] = Vz[C_bsver[i][k]];
				}
				
				Cell_Apical[i] = Area_Calc(AP);//A_0 = A_from_cell(AP, BS);
				Cell_Basal[i] = Area_Calc(BS);//B_0 = B_from_cell(AP, BS);
				Cell_Volume[i] = Volume_Calc(AP, BS);//V_0 = V_from_cell(AP, BS);
				Cell_Lengths_1[i] = sqrt(pow(AP[0][0] - BS[0][0], 2) + pow(AP[0][1] - BS[0][1], 2) + pow(AP[0][2] - BS[0][2], 2));
				Cell_Lengths_2[i] = sqrt(pow(AP[1][0] - BS[1][0], 2) + pow(AP[1][1] - BS[1][1], 2) + pow(AP[1][2] - BS[1][2], 2));
				Cell_Lengths_3[i] = sqrt(pow(AP[2][0] - BS[2][0], 2) + pow(AP[2][1] - BS[2][1], 2) + pow(AP[2][2] - BS[2][2], 2));
				Cell_Lengths_4[i] = sqrt(pow(AP[3][0] - BS[3][0], 2) + pow(AP[3][1] - BS[3][1], 2) + pow(AP[3][2] - BS[3][2], 2));
				
				E_Apical[i] = K_par[0]*pow(Cell_Apical[i] - A_0_cell[i], 2);
				E_Basal[i] = K_par[1]*pow(Cell_Basal[i] - B_0_cell[i], 2);
				E_Sides[i] = K_par[3]*pow(Cell_Lengths_1[i] - L_0_cell[i], 2);
				E_Sides[i] += K_par[4]*pow(Cell_Lengths_2[i] - L_0_cell[i], 2);
				E_Sides[i] += K_par[5]*pow(Cell_Lengths_3[i] - L_0_cell[i], 2);
				E_Sides[i] += K_par[6]*pow(Cell_Lengths_4[i] - L_0_cell[i], 2);
				E_Volume[i] = K_par[2]*pow(Cell_Volume[i] - V_0_cell[i], 2);				
			}
			
			Record_Vert_data(Vx, Vy, Vz, Nv, count, p);
			Record_Cell_data(Cell_Apical, Cell_Basal, Cell_Volume, Cell_Lengths_1, Cell_Lengths_2, Cell_Lengths_3, Cell_Lengths_4, Yolk_V_tmp, Nc, count, p);
			Record_Energy_data(E_Apical, E_Basal, E_Sides, E_Volume, E_Planes, Nc, Nv, count++, p);
			
		
		
		}
		
		if (E_total_old < E_total && t>0) 
		cout << "WARNING: Energy is increasing." << endl;
		
		if (Percent_stop_flag == 1 && fabs(E_total - E_total_old)/E_total < Percent_stop) {
			cout << "Energy change is less than .1%" << endl;
			break;
		}
		
		E_total_old = E_total;
		t++;
		
		
		}

*/