//=====================================================================================================================================
// Assing initial values for the verticies
int id_tmp = 0;
double Angle_start = 3.14159265/2.0 - Nx/2.0*Theta_step;
for (k=0; k<2; k++) {
	for (j=0; j<Ny; j++) {
		for (i=0; i<Nx; i++) {
			id_tmp = i + Nx*j + Nx*Ny*k;
			
			
			if (Sheet_condition == 1) {
				if (k==0) {
					Vx[id_tmp] = r_out*cos(i*Theta_step + Angle_start);
					Vy[id_tmp] = j*sqrt(pow(r_out*(cos(2*Theta_step) - cos(Theta_step)), 2) + pow(r_out*(sin(2*Theta_step) - sin(Theta_step)), 2));
					Vz[id_tmp] = r_out*sin(i*Theta_step+ Angle_start);
				}
				else {
					Vx[id_tmp] = r_in*cos(i*Theta_step + Angle_start);
					Vy[id_tmp] = j*sqrt(pow(r_out*(cos(2*Theta_step) - cos(Theta_step)), 2) + pow(r_out*(sin(2*Theta_step) - sin(Theta_step)), 2));;
					Vz[id_tmp] = r_in*sin(i*Theta_step + Angle_start);
				}
			}
			else if (Sheet_condition == 2){
				
				if (fixed_first_last_vert_xyz == 1) {
					if (i==0) {
						Vx[id_tmp] = i*X_spacing - 3*X_spacing;
					}
					else if (i==(Nx-1)){
						Vx[id_tmp] = i*X_spacing + 3*X_spacing;
					}
					else {
						Vx[id_tmp] = i*X_spacing;
					}
					
				}
				else {
					Vx[id_tmp] = i*X_spacing;
				}
				
				Vy[id_tmp] = j*Y_spacing;
				
				if (k==0) {
					Vz[id_tmp] = Z_spacing;
				}
				else {
					Vz[id_tmp] = 0;
				}
			}
			
			// Now assing the index of the cells that are affected by this vertex
			if ((i>0) && (i < (Nx-1))) {
				if ((j>0) && (j<(Ny-1))) {
					V_cells[id_tmp][0] = Cell_Id(i, j, Nx);
					V_cells[id_tmp][1] = Cell_Id(i-1, j, Nx);
					V_cells[id_tmp][2] = Cell_Id(i, j-1, Nx);
					V_cells[id_tmp][3] = Cell_Id(i-1, j-1, Nx);
				}
				else if (j==0){
					V_cells[id_tmp][0] = Cell_Id(i, j, Nx);
					V_cells[id_tmp][1] = Cell_Id(i-1, j, Nx);
					V_cells[id_tmp][2] = -1;
					V_cells[id_tmp][3] = -1;
				}
				else if (j==(Ny-1)){
					V_cells[id_tmp][0] = -1;
					V_cells[id_tmp][1] = -1;
					V_cells[id_tmp][2] = Cell_Id(i, j-1, Nx);
					V_cells[id_tmp][3] = Cell_Id(i-1, j-1, Nx);
				}
			}
			else if (i==0) {
				if ((j>0) && (j<(Ny-1))) {
					V_cells[id_tmp][0] = Cell_Id(i, j, Nx);
					V_cells[id_tmp][1] = -1;
					V_cells[id_tmp][2] = Cell_Id(i, j-1, Nx);
					V_cells[id_tmp][3] = -1;
				}
				else if (j==0){
					V_cells[id_tmp][0] = Cell_Id(i, j, Nx);
					V_cells[id_tmp][1] = -1;
					V_cells[id_tmp][2] = -1;
					V_cells[id_tmp][3] = -1;
				}
				else if (j==(Ny-1)){
					V_cells[id_tmp][0] = -1;
					V_cells[id_tmp][1] = -1;
					V_cells[id_tmp][2] = Cell_Id(i, j-1, Nx);
					V_cells[id_tmp][3] = -1;
				}
			}
			else if (i==(Nx-1)){
				if ((j>0) && (j<(Ny-1))) {
					V_cells[id_tmp][0] = -1;
					V_cells[id_tmp][1] = Cell_Id(i-1, j, Nx);
					V_cells[id_tmp][2] = -1;
					V_cells[id_tmp][3] = Cell_Id(i-1, j-1, Nx);
				}
				else if (j==0){
					V_cells[id_tmp][0] = -1;
					V_cells[id_tmp][1] = Cell_Id(i-1, j, Nx);
					V_cells[id_tmp][2] = -1;
					V_cells[id_tmp][3] = -1;
				}
				else if (j==(Ny-1)){
					V_cells[id_tmp][0] = -1;
					V_cells[id_tmp][1] = -1;
					V_cells[id_tmp][2] = -1;
					V_cells[id_tmp][3] = Cell_Id(i-1, j-1, Nx);
				}
			}
			
			
			// Define the movement properties of the verticies
			// Fix the first and the last y layer to lie within the same y plane
			if ((fixed_y == 1) && ((j==0) || (j==(Ny-1)))) {
				V_not_move_y[id_tmp] = 0;
			}
			else {
				V_not_move_y[id_tmp] = 1;
			}
			
			// Fix the first and the last cells to not move along z
			if ((fixed_first_last_cell_z == 1) && ((i==0) || (i==1) || (i==(Nx-1)) || (i==(Nx-2)))) {
				V_not_move_z[id_tmp] = 0;
			}
			else {
				V_not_move_z[id_tmp] = 1;
			}
			
			
			
			if ((fixed_first_last_vert_xyz == 1) && ((i==0) || (i == (Nx-1)))) {
				V_not_move_x[id_tmp] = 0;
				V_not_move_y[id_tmp] = 0;
				V_not_move_z[id_tmp] = 0;
			}
			else {
				V_not_move_x[id_tmp] = 1;
			}
		}
	}
}

/*
 cout << "V_not_move_x = [";
 for (i=0; i<Nv; i++) {
 cout << V_not_move_x[i] << "  ";
 }
 cout << "]" << endl;
 
 cout << "V_not_move_y = [";
 for (i=0; i<Nv; i++) {
 cout << V_not_move_y[i] << "  ";
 }
 cout << "]" << endl;
 
 cout << "V_not_move_z = [";
 for (i=0; i<Nv; i++) {
 cout << V_not_move_z[i] << "  ";
 }
 cout << "]" << endl;
 */

//=====================================================================================================================================
//=====================================================================================================================================
// Assing values for the yolk sheet

for (j=0; j<Ny; j++) {
	for (i=0; i<Nx; i++) {
		id_tmp = i + Nx*j;
		
		
		if (Sheet_condition == 1) {
			
			Yolk_x[id_tmp] = 0;
			Yolk_y[id_tmp] = j*Y_spacing;
			Yolk_z[id_tmp] = 0;
		}
		else if (Sheet_condition == 2)
		{
			if (fixed_first_last_vert_xyz == 1) {
				if (i==0) {
					Yolk_x[id_tmp] = i*X_spacing - 3*X_spacing;
				}
				else if (i==(Nx-1)){
					Yolk_x[id_tmp] = i*X_spacing + 3*X_spacing;
				}
				else {
					Yolk_x[id_tmp] = i*X_spacing;
				}
				
			}
			else {
				Yolk_x[id_tmp] = i*X_spacing;
			}
			
			
			Yolk_y[id_tmp] = j*Y_spacing;
			Yolk_z[id_tmp] = -4*Z_spacing;
		}
		
	}
}


//=====================================================================================================================================
//=====================================================================================================================================
// Assing Id to the cells
for (j=0; j<(Ny-1); j++) {
	for (i=0; i<(Nx-1); i++) {
		id_tmp = Cell_Id(i, j, Nx);
		C_apver[id_tmp][0] = Ver_Id(i, j, 0, Nx, Ny);
		C_apver[id_tmp][1] = Ver_Id(i+1, j, 0, Nx, Ny);
		C_apver[id_tmp][2] = Ver_Id(i, j+1, 0, Nx, Ny);
		C_apver[id_tmp][3] = Ver_Id(i+1, j+1, 0, Nx, Ny);
		
		C_bsver[id_tmp][0] = Ver_Id(i, j, 1, Nx, Ny);
		C_bsver[id_tmp][1] = Ver_Id(i+1, j, 1, Nx, Ny);
		C_bsver[id_tmp][2] = Ver_Id(i, j+1, 1, Nx, Ny);
		C_bsver[id_tmp][3] = Ver_Id(i+1, j+1, 1, Nx, Ny);
		
		C_ykver[id_tmp][0] = Ver_Id(i, j, 0, Nx, Ny);
		C_ykver[id_tmp][1] = Ver_Id(i+1, j, 0, Nx, Ny);
		C_ykver[id_tmp][2] = Ver_Id(i, j+1, 0, Nx, Ny);
		C_ykver[id_tmp][3] = Ver_Id(i+1, j+1, 0, Nx, Ny);
	}
}

//====================================================================================================================================
// Assing the index to the points that determine the contraint

for (k=0; k<2; k++) {
	for (j=0; j<Ny; j++) {
		for (i=0; i<Nx; i++) {
			id_tmp = Ver_Id(i, j, k, Nx, Ny);
			
			// Left side
			Plane_test_p1[id_tmp][0] = Ver_Id(i-1, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][0] = Ver_Id(i-1, j, 0, Nx, Ny);
			Plane_test_p3[id_tmp][0] = Ver_Id(i-1, j-1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][1] = Ver_Id(i-1, j-1, 1, Nx, Ny); 
			Plane_test_p2[id_tmp][1] = Ver_Id(i-1, j, 1, Nx, Ny);
			Plane_test_p3[id_tmp][1] = Ver_Id(i-1, j, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][2] = Ver_Id(i-1, j+1, 1, Nx, Ny); 
			Plane_test_p2[id_tmp][2] = Ver_Id(i-1, j, 1, Nx, Ny);
			Plane_test_p3[id_tmp][2] = Ver_Id(i-1, j, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][3] = Ver_Id(i-1, j, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][3] = Ver_Id(i-1, j+1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][3] = Ver_Id(i-1, j+1, 0, Nx, Ny);
			
			// Right side
			Plane_test_p1[id_tmp][4] = Ver_Id(i+1, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][4] = Ver_Id(i+1, j, 0, Nx, Ny);
			Plane_test_p3[id_tmp][4] = Ver_Id(i+1, j-1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][5] = Ver_Id(i+1, j-1, 1, Nx, Ny); 
			Plane_test_p2[id_tmp][5] = Ver_Id(i+1, j, 1, Nx, Ny);
			Plane_test_p3[id_tmp][5] = Ver_Id(i+1, j, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][6] = Ver_Id(i+1, j+1, 1, Nx, Ny); 
			Plane_test_p2[id_tmp][6] = Ver_Id(i+1, j, 1, Nx, Ny);
			Plane_test_p3[id_tmp][6] = Ver_Id(i+1, j, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][7] = Ver_Id(i+1, j, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][7] = Ver_Id(i+1, j+1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][7] = Ver_Id(i+1, j+1, 0, Nx, Ny);
			
			// Forward side
			Plane_test_p1[id_tmp][8] = Ver_Id(i-1, j+1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][8] = Ver_Id(i-1, j+1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][8] = Ver_Id(i, j+1, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][9] = Ver_Id(i, j+1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][9] = Ver_Id(i, j+1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][9] = Ver_Id(i-1, j+1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][10] = Ver_Id(i, j+1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][10] = Ver_Id(i, j+1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][10] = Ver_Id(i+1, j+1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][11] = Ver_Id(i, j+1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][11] = Ver_Id(i+1, j+1, 0, Nx, Ny);
			Plane_test_p3[id_tmp][11] = Ver_Id(i+1, j+1, 1, Nx, Ny);
			
			// Back side
			Plane_test_p1[id_tmp][12] = Ver_Id(i-1, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][12] = Ver_Id(i-1, j-1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][12] = Ver_Id(i, j-1, 0, Nx, Ny);
			
			Plane_test_p1[id_tmp][13] = Ver_Id(i, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][13] = Ver_Id(i, j-1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][13] = Ver_Id(i-1, j-1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][14] = Ver_Id(i, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][14] = Ver_Id(i, j-1, 1, Nx, Ny);
			Plane_test_p3[id_tmp][14] = Ver_Id(i+1, j-1, 1, Nx, Ny);
			
			Plane_test_p1[id_tmp][15] = Ver_Id(i, j-1, 0, Nx, Ny); 
			Plane_test_p2[id_tmp][15] = Ver_Id(i+1, j-1, 0, Nx, Ny);
			Plane_test_p3[id_tmp][15] = Ver_Id(i+1, j-1, 1, Nx, Ny);
			
			
			if (k==0) {
			// Bottom side
				Plane_test_p1[id_tmp][16] = Ver_Id(i-1, j-1, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][16] = Ver_Id(i, j-1, 1, Nx, Ny);
				Plane_test_p3[id_tmp][16] = Ver_Id(i-1, j, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][17] = Ver_Id(i-1, j, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][17] = Ver_Id(i, j-1, 1, Nx, Ny);
				Plane_test_p3[id_tmp][17] = Ver_Id(i, j, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][18] = Ver_Id(i, j+1, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][18] = Ver_Id(i-1, j, 1, Nx, Ny);
				Plane_test_p3[id_tmp][18] = Ver_Id(i, j, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][19] = Ver_Id(i-1, j, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][19] = Ver_Id(i-1, j+1, 1, Nx, Ny);
				Plane_test_p3[id_tmp][19] = Ver_Id(i, j+1, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][20] = Ver_Id(i+1, j-1, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][20] = Ver_Id(i+1, j, 1, Nx, Ny);
				Plane_test_p3[id_tmp][20] = Ver_Id(i, j-1, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][21] = Ver_Id(i, j, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][21] = Ver_Id(i, j-1, 1, Nx, Ny);
				Plane_test_p3[id_tmp][21] = Ver_Id(i+1, j, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][22] = Ver_Id(i+1, j+1, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][22] = Ver_Id(i+1, j, 1, Nx, Ny);
				Plane_test_p3[id_tmp][22] = Ver_Id(i, j, 1, Nx, Ny);
				
				Plane_test_p1[id_tmp][23] = Ver_Id(i+1, j+1, 1, Nx, Ny); 
				Plane_test_p2[id_tmp][23] = Ver_Id(i, j, 1, Nx, Ny);
				Plane_test_p3[id_tmp][23] = Ver_Id(i, j+1, 1, Nx, Ny);
			}
			else {
				// Top side
				Plane_test_p1[id_tmp][16] = Ver_Id(i-1, j-1, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][16] = Ver_Id(i, j-1, 0, Nx, Ny);
				Plane_test_p3[id_tmp][16] = Ver_Id(i-1, j, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][17] = Ver_Id(i-1, j, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][17] = Ver_Id(i, j-1, 0, Nx, Ny);
				Plane_test_p3[id_tmp][17] = Ver_Id(i, j, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][18] = Ver_Id(i, j+1, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][18] = Ver_Id(i-1, j, 0, Nx, Ny);
				Plane_test_p3[id_tmp][18] = Ver_Id(i, j, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][19] = Ver_Id(i-1, j, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][19] = Ver_Id(i-1, j+1, 0, Nx, Ny);
				Plane_test_p3[id_tmp][19] = Ver_Id(i, j+1, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][20] = Ver_Id(i+1, j-1, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][20] = Ver_Id(i+1, j, 0, Nx, Ny);
				Plane_test_p3[id_tmp][20] = Ver_Id(i, j-1, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][21] = Ver_Id(i, j, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][21] = Ver_Id(i, j-1, 0, Nx, Ny);
				Plane_test_p3[id_tmp][21] = Ver_Id(i+1, j, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][22] = Ver_Id(i+1, j+1, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][22] = Ver_Id(i+1, j, 0, Nx, Ny);
				Plane_test_p3[id_tmp][22] = Ver_Id(i, j, 0, Nx, Ny);
				
				Plane_test_p1[id_tmp][23] = Ver_Id(i+1, j+1, 0, Nx, Ny); 
				Plane_test_p2[id_tmp][23] = Ver_Id(i, j, 0, Nx, Ny);
				Plane_test_p3[id_tmp][23] = Ver_Id(i, j+1, 0, Nx, Ny);
				
			}


			
		}
	}
}


//====================================================================================================================================
// Assing the right inital parameters to each cell
// Calculate A_0, B_0, V_0, L_0 for a cell

for (i=0; i<Nc; i++) {
	Cell_sel = i;
	
	for (k=0; k<4; k++) {
		AP[k][0] = Vx[C_apver[Cell_sel][k]]; AP[k][1] = Vy[C_apver[Cell_sel][k]]; AP[k][2] = Vz[C_apver[Cell_sel][k]];
		BS[k][0] = Vx[C_bsver[Cell_sel][k]]; BS[k][1] = Vy[C_bsver[Cell_sel][k]]; BS[k][2] = Vz[C_bsver[Cell_sel][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));
	
	Cell_Diag_1[i] = sqrt(pow(AP[0][0] - BS[2][0], 2) + pow(AP[0][1] - BS[2][1], 2) + pow(AP[0][2] - BS[2][2], 2));
	Cell_Diag_2[i] = sqrt(pow(AP[1][0] - BS[3][0], 2) + pow(AP[1][1] - BS[3][1], 2) + pow(AP[1][2] - BS[3][2], 2));
	Cell_Diag_3[i] = sqrt(pow(AP[2][0] - BS[0][0], 2) + pow(AP[2][1] - BS[0][1], 2) + pow(AP[2][2] - BS[0][2], 2));
	Cell_Diag_4[i] = sqrt(pow(AP[3][0] - BS[1][0], 2) + pow(AP[3][1] - BS[1][1], 2) + pow(AP[3][2] - BS[1][2], 2));
	
	//cout << "For cell " << i << " A = " << Cell_Apical[i] << " B_0 = " << Cell_Basal[i] << " V_0 = " << Cell_Volume[i] << " L_0 = " << Cell_Lengths_1[i] << endl;
	}
	
	
	// Declare the initial yolk volume
	double Yolk_O = 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_O += Volume_Calc(AP, BS);
	}
	
	// Display and compare to the theoretical values
	cout << endl << "A_0 = " << Cell_Apical[0] <<endl;
	cout << "B_0 = " << Cell_Basal[0] <<endl;
	cout << "V_0 = " << Cell_Volume[0] <<endl;
	cout << "L_0 = " << Cell_Lengths_1[0] <<endl;
	cout << "Y_0 = " << Yolk_O << endl << endl;
	
	cout << "K_apical = " << K_apical << endl;
	cout << "K_basal = " << K_basal << endl;
	cout << "K_volume = " << K_volume << endl;
	cout << "K_length = " << K_length << endl;
	cout << "K_yolk = " << K_yolk << endl << endl << endl;
	
	
	
	for (j=0; j<(Ny-1); j++) {
		for (i=0; i<(Nx-1); i++) {
			id_tmp = Cell_Id(i, j, Nx);
			// Assing the parameters
			// Assing the weights for each sides to compensate for ones that are counted twice
			if (i==0) {
				if (j==0) {
					C_side_factor[id_tmp][0] = 4; C_side_factor[id_tmp][1] = 2; C_side_factor[id_tmp][2] = 2; C_side_factor[id_tmp][3] = 1;
				}
				else if (j==(Ny-2)){
					C_side_factor[id_tmp][0] = 2; C_side_factor[id_tmp][1] = 1; C_side_factor[id_tmp][2] = 4; C_side_factor[id_tmp][3] = 2;
				}
				else {
					C_side_factor[id_tmp][0] = 2; C_side_factor[id_tmp][1] = 1; C_side_factor[id_tmp][2] = 2; C_side_factor[id_tmp][3] = 1;
				}
			}
			else if (i==(Nx-2)) {
				if (j==0) {
					C_side_factor[id_tmp][0] = 2; C_side_factor[id_tmp][1] = 4; C_side_factor[id_tmp][2] = 1; C_side_factor[id_tmp][3] = 2;
				}
				else if (j==(Ny-2)){
					C_side_factor[id_tmp][0] = 1; C_side_factor[id_tmp][1] = 2; C_side_factor[id_tmp][2] = 2; C_side_factor[id_tmp][3] = 4;
				}
				else {
					C_side_factor[id_tmp][0] = 1; C_side_factor[id_tmp][1] = 2; C_side_factor[id_tmp][2] = 1; C_side_factor[id_tmp][3] = 2;
				}
			}
			else if ((j==0) && (i>0) && (i<(Nx-2))){
				C_side_factor[id_tmp][0] = 2; C_side_factor[id_tmp][1] = 2; C_side_factor[id_tmp][2] = 1; C_side_factor[id_tmp][3] = 1;
			}
			else if ((j==(Ny-2)) && (i>0) && (i<(Nx-2))){
				C_side_factor[id_tmp][0] = 1; C_side_factor[id_tmp][1] = 1; C_side_factor[id_tmp][2] = 2; C_side_factor[id_tmp][3] = 2;
			}
			else {
				C_side_factor[id_tmp][0] = 1; C_side_factor[id_tmp][1] = 1; C_side_factor[id_tmp][2] = 1; C_side_factor[id_tmp][3] = 1;
			}
			
			// Assing which cells will constrict
			if (gaussian_KA_flag == 1){
				if ((i>= Constrict_start) && (i<= Constrict_end) && (j >= y_start) && (j <= y_stop)) {
					A_0_cell[id_tmp] = A_con_factor*Cell_Apical[id_tmp];
					//A_K_cell[id_tmp] = K_apical + KA_con_factor*K_apical*exp(-pow((i-(Nx-1)/2.0)/KA_sigma, 2));
				}
				else {
					A_0_cell[id_tmp] = Cell_Apical[id_tmp];
					//A_K_cell[id_tmp] = K_apical;
				}
				
				/*
				if ((j > y_start) && (j < y_stop)) {
					A_K_cell[id_tmp] = K_apical + KA_con_factor*K_apical*exp(-pow((i + 0.5 - (Nx-1)/2.0)/KA_sigma, 2));
				}
				else if ((j == y_start) || (j == y_stop)) {
					A_K_cell[id_tmp] = K_apical + 0.5*KA_con_factor*K_apical*exp(-pow((i + 0.5 - (Nx-1)/2.0)/KA_sigma, 2));
				}
				else {
					A_K_cell[id_tmp] = K_apical;
				}
				 */
				A_K_cell[id_tmp] = K_apical + KA_con_factor*K_apical*exp(-pow((i + 0.5 - (Nx-1)/2.0)/KA_sigma, 2));
				

				//A_K_cell[id_tmp] = K_apical + KA_con_factor*K_apical/2.0*(1.0/(1.0 + exp(-(i-Constrict_start)*KA_sigma)) + 1.0/(1.0 + exp((i-Constrict_end)*KA_sigma)));
			}
			else {
				if ((i>= Constrict_start) && (i<= Constrict_end) && (j >= y_start) && (j <= y_stop)) {
					A_0_cell[id_tmp] = A_con_factor*Cell_Apical[id_tmp];
					A_K_cell[id_tmp] = KA_con_factor*K_apical;
				}
				else {
					A_0_cell[id_tmp] = Cell_Apical[id_tmp];
					A_K_cell[id_tmp] = K_apical;
				}
			}
			
			// Assing the values of A_K, B_K, V_K, L_K to each cell
			
			if (neighbor_apical_active == 1) {
				if ((i >= (Constrict_start - neighbor_cell_N) && i < Constrict_start) || (i <= (Constrict_end + neighbor_cell_N) && i > Constrict_end) && (j >= y_start) && (j <= y_stop)) {
					A_0_cell[id_tmp] = A_con_factor*Cell_Apical[id_tmp];
				}
			}
			
			B_K_cell[id_tmp] = K_basal;
			B_0_cell[id_tmp] = Cell_Basal[id_tmp];
			
			if (neighbor_basal_strong == 1) {
				if ((i >= (Constrict_start - neighbor_cell_N) && i < Constrict_start) || (i <= (Constrict_end + neighbor_cell_N) && i > Constrict_end) && (j >= y_start) && (j <= y_stop)) {
					B_K_cell[id_tmp] = KA_con_factor*K_apical;
				}
			}
			
			V_K_cell[id_tmp] = K_volume;
			V_0_cell[id_tmp] = Cell_Volume[id_tmp];
			
			L_K_cell[id_tmp] = K_length;
			L_0_cell[id_tmp] = Cell_Lengths_1[id_tmp];
			
			D_K_cell[id_tmp] = K_diagonal;
			D_0_cell[id_tmp] = Cell_Diag_1[id_tmp];
			
			//C_side_factor[id_tmp][0] = 1; C_side_factor[id_tmp][1] = 1; C_side_factor[id_tmp][2] = 1; C_side_factor[id_tmp][3] = 1;
			
		}
	}
	
	if (Ny == 2){
		A_K_cell[0] /= 100; A_K_cell[Nc-1] /= 100;
		B_K_cell[0] /= 100; B_K_cell[Nc-1] /= 100;
		V_K_cell[0] /= 100; V_K_cell[Nc-1] /= 100;
		L_K_cell[0] /= 100; L_K_cell[Nc-1] /= 100;
	}
	
	K_par[0] = A_K_cell[1];
	K_par[1] = B_K_cell[1];
	K_par[2] = V_K_cell[1];
	K_par[3] = L_K_cell[1];
	K_par[4] = L_K_cell[1];
	K_par[5] = L_K_cell[1];
	K_par[6] = L_K_cell[1];
	K_par[7] = D_K_cell[1];
	
	O_par[0] = A_0_cell[1];
	O_par[1] = B_0_cell[1];
	O_par[2] = V_0_cell[1];
	O_par[3] = L_0_cell[1];
	O_par[4] = L_0_cell[1];
	O_par[5] = L_0_cell[1];
	O_par[6] = L_0_cell[1];
	O_par[7] = D_0_cell[1];
	
	other_par[0] = y_start;
	other_par[1] = y_stop;
	other_par[2] = Constrict_start;
	other_par[3] = Constrict_end;
	other_par[4] = A_0_cell[1]*A_con_factor;
	other_par[5] = K_yolk;
	other_par[6] = Yolk_O;
	
	for (i=0; i<Nc; i++) {
		cout << "For cell " << i << " A_0 = " << A_0_cell[i] << " B_0 = " << B_0_cell[i] << " V_0 = " << V_0_cell[i] << " L_0 = " << L_0_cell[i];
		//cout << "For cell " << i << " A_K = " << A_K_cell[i] << " B_K = " << B_K_cell[i] << " V_K = " << V_K_cell[i] << " L_K = " << L_K_cell[i] << endl << end;
		cout << " A_K = " << A_K_cell[i] << " B_K = " << B_K_cell[i] << " V_K = " << V_K_cell[i] << " L_K = " << L_K_cell[i] << endl;
	}
	
	
	
	