int i=0, j=0, k=0, t=0, count = 1;


const int Nv = 2*Nx*Ny; // This is the total number of verticies
const int Nc = (Nx-1)*(Ny-1);// This is the number of cells

// Define initial values of the values

double A_0;
double V_0;
double B_0;
double L_0;

// Define Par so you could easily pass it into the function
double* K_par = new double[10];
double* O_par = new double[8];

double* other_par = new double[7];

// Declare arrays of verticies position
double* Vx = new double[Nv];
double* Vy = new double[Nv];
double* Vz = new double[Nv];

// For parallel processing
int* V_slave_id = new int[Nv]; // Which slave will work on which V
int* C_slave_id = new int[Nc]; // Which slave will work on which V

// Declare arrays of yolk positions
double* Yolk_x = new double[Nv/2];
double* Yolk_y = new double[Nv/2];
double* Yolk_z = new double[Nv/2];

int* V_not_move_x = new int[Nv];
int* V_not_move_y = new int[Nv];
int* V_not_move_z = new int[Nv];

double* Vx_nxt = new double[Nv];
double* Vy_nxt = new double[Nv];
double* Vz_nxt = new double[Nv];

double* Vx_prv = new double[Nv];
double* Vy_prv = new double[Nv];
double* Vz_prv = new double[Nv];

double E_nxt_x = 0;
double E_prv_x = 0;
double E_nxt_y = 0;
double E_prv_y = 0;
double E_nxt_z = 0;
double E_prv_z = 0;

double E_total = 0;
double E_total_old = 0;

double E_new = 0;
double E_old = 0;
int Cell_sel = 0;


double V_grad_mag = 0;
double V_grad_tmp = 0;

long * idum = new long;
*(idum) = 100;

double** V_grad = new double*[Nv];
for (i=0; i<Nv; i++) {
	V_grad[i] = new double[3];
}


int** V_cells = new int*[Nv]; // The 4 cells that are affected by each vertex
int** Plane_test_p1 = new int*[Nv]; // The index of the vertices that define the plane which bound each point
int** Plane_test_p2 = new int*[Nv];
int** Plane_test_p3 = new int*[Nv];

for (i=0; i<Nv; i++) {
	V_cells[i] = new int[4];
	Plane_test_p1[i] = new int[24];
	Plane_test_p2[i] = new int[24];
	Plane_test_p3[i] = new int[24];
}

double* Point1 = new double[3];
double* Point2 = new double[3];
double* Point3 = new double[3];
double* Point0 = new double[3];

// Declare arrays of cell parameters
int** C_apver = new int*[Nc]; // indices of the apical verticies
int** C_bsver = new int*[Nc]; // indicies of the basal verticies
int** C_ykver = new int*[Nc];
int** C_side_factor = new int*[Nc];
for (i=0; i<Nc; i++) {
	C_apver[i] = new int[4];
	C_bsver[i] = new int[4];
	C_ykver[i] = new int[4];
	C_side_factor[i] = new int[4];
}

double* A_K_cell = new double[Nc]; // A_K for each cell
double* A_0_cell = new double[Nc]; // A0 for each cell
double* V_K_cell = new double[Nc]; // V_K for each cell
double* V_0_cell = new double[Nc]; // V0 for each cell
double* B_K_cell = new double[Nc]; // V_K for each cell
double* B_0_cell = new double[Nc]; // V0 for each cell
double* L_K_cell = new double[Nc]; // L_K for each cell
double* L_0_cell = new double[Nc]; // L0 for each cell
double* D_K_cell = new double[Nc]; // D_K for each cell
double* D_0_cell = new double[Nc]; // D0 for each cell

// Declare arrays that hold the parameters of the cell
double* Cell_Apical = new double[Nc];
double* Cell_Basal = new double[Nc];
double* Cell_Volume = new double[Nc];
double* Cell_Lengths_1 = new double[Nc];
double* Cell_Lengths_2 = new double[Nc];
double* Cell_Lengths_3 = new double[Nc];
double* Cell_Lengths_4 = new double[Nc];

double* Cell_Diag_1 = new double[Nc];
double* Cell_Diag_2 = new double[Nc];
double* Cell_Diag_3 = new double[Nc];
double* Cell_Diag_4 = new double[Nc];

double* Yolk_volumes = new double[Nc];
double Yolk_V = 0;

double* E_Apical = new double[Nc]; // A_K for each cell
double* E_Basal = new double[Nc]; // A0 for each cell
double* E_Sides = new double[Nc]; // V_K for each cell
double* E_Volume = new double[Nc]; // V0 for each cell
double* E_Planes = new double[Nv]; // V_K for each cell

// Declare variables for computation
// Apical verticies and basal verticies

double** AP = new double*[4];
double** BS = new double*[4];

for (i=0; i<4; i++) {
	AP[i] = new double[3];
	BS[i] = new double[3];
}

double R_tmp = 0;
double Dis_tmp = 0;
double Yolk_V_tmp = 0;

// Stuff for the potential
double R_max = .1;
const int N_plane = 10000;
double dx_plane = R_max/N_plane;

double par_a = 1000;
double par_b = 1000;

double* Plane_potential = new double[N_plane];

for (i=0;i<N_plane;i++){
	Plane_potential[i] = par_a/(1 + exp(-(-dx_plane*i + R_max/2.0)*par_b));
}



