for (t=0;t<T_end; t++){
	
	// Calculate initial yolk volumes
	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];
	}
	
	//cout << "Yolk Volume is " << Yolk_V << endl;
	
	// Store the previous and the next vectors
	for (i=0; i<Nv; i++) {
		Vx_nxt[i] = Vx[i];
		Vy_nxt[i] = Vy[i];
		Vz_nxt[i] = Vz[i];
		
		Vx_prv[i] = Vx[i];
		Vy_prv[i] = Vy[i];
		Vz_prv[i] = Vz[i];
	}
	
	
	
	// Loop over all verticies and calculate the gradient
	#include "Gradient_Calc.h"
	
	//*
	// Scale the gradient so that the maximum step would be equal to dx
	V_grad_mag = 0;
	for (i=0; i<Nv; i++) {// Calculate the maximum value of the step
		V_grad_tmp = sqrt(V_grad[i][0]*V_grad[i][0] + V_grad[i][1]*V_grad[i][1] + V_grad[i][2]*V_grad[i][2]);
		if (V_grad_tmp > V_grad_mag) {
			V_grad_mag = V_grad_tmp;
		}
	}
	
	for (i=0; i<Nv; i++) { // Rescale the gradient
		V_grad[i][0] *= dx/V_grad_mag;
		V_grad[i][1] *= dx/V_grad_mag;
		V_grad[i][2] *= dx/V_grad_mag;
	}
	
	// Calculate the energy
	E_old = 0;
	// Energy from cells
	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];
		
		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];
		E_old += E_from_cell_V2(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_old += 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_old += K_yolk*pow(Yolk_V - Yolk_O, 2);
	}
	
	
	
	// Now keep moving the verticies until the energy is no longer decreasing
	
	while (1){
		//cout << "E_old = " << E_old << "  E_new = " << E_new << "  t = " << t <<  endl;
		
		//E_old = E_new;
		// Record the old positions 
		for (i=0; i<Nv; i++) {
			Vx_prv[i] = Vx[i];
			Vy_prv[i] = Vy[i];
			Vz_prv[i] = Vz[i];
		}
		
		// Move the verticies
		// Move Apical
		for (i=0; i<(Nv/2); i++) {
			Vx[i] += V_grad[i][0]*V_not_move_x[i];//*V_not_move[i];
			Vy[i] += V_grad[i][1]*V_not_move_y[i];
			Vz[i] += V_grad[i][2]*V_not_move_z[i];//*V_not_move[i];
		}
		
		// Move Basal
		for (i=Nv/2; i<Nv; i++) {
			Vx[i] += V_grad[i][0]*V_not_move_x[i]*Basal_friction;//*V_not_move[i];
			Vy[i] += V_grad[i][1]*V_not_move_y[i]*Basal_friction;
			Vz[i] += V_grad[i][2]*V_not_move_z[i]*Basal_friction;//*V_not_move[i];
		}
		
		
		// Apply boundary conditions 
		if (fixed_top == 1) {
			if (Sheet_condition == 2) {
				for (i=0; i<Nv; i++) {
					if (Vz[i] > Z_spacing) {
						Vz[i] = Z_spacing;
					}
				}
			}
			else if (Sheet_condition == 1){
				for (i=0; i<Nv; i++) {
					R_tmp = sqrt(Vx[i]*Vx[i] + Vy[i]*Vy[i] + Vz[i]*Vz[i]);
					if (R_tmp > r_out) {
						Vx[i] = Vx[i]*r_out/R_tmp;
						Vy[i] = Vy[i]*r_out/R_tmp;
						Vz[i] = Vz[i]*r_out/R_tmp;
					}
				}
			}
		}
		
		// Calculate the new energy
		E_new = 0;
		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]];
			}
			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];
			
			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];
			E_new += E_from_cell_V2(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_new += 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_new += K_yolk*pow(Yolk_V - Yolk_O, 2);
		}
		
		if (E_new > E_old) {
			break;
		}
		else {
			E_old = E_new;
		}
		
		
		//cout << "E_old = " << E_old << "  E_new = " << E_new << "  t = " << t <<  endl << endl;
	}
	
	
	
	for (i=0; i<Nv; i++) {
		Vx[i] = Vx_prv[i];
		Vy[i] = Vy_prv[i];
		Vz[i] = Vz_prv[i];
	}
	
	
	
	// 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];
		
		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];
		
		E_total += E_from_cell_V2(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) << endl;
		
		if (t % T_record == 0){
			
			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));
			}
			
			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);
		}
		
		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;
		
		
		}