
#include <math.h>
#include <string.h>

class Embryo{
	/*
	Class object to implement the model
	By Oleg Polyakov | Princeton University | 2008-2013

	*/

	double **Vx;
	double **Vy;

	double **Gx;
	double **Gy;

	double **La;
	double **Ka;
	double **Ka_con;
	double **Kl;
	
	double **A0;
	double **A0_con;
	double **L0;

	double **Atmp;
	double **Ltmp;

	double *Kv;
	double *Vtmp;
	double *V0;

	double Ky;
	double Y0;

	double K_angle;

	double *R;

	int N;
	int N_layers;

	double Alpha;
	double Sigma;

	double mu;

	double Energy;
	bool Energy_initialized;

	char * filename;
	int loop_count;

public:
	Embryo(int N_n, int N_l){
		N = N_n;
		N_layers = N_l;

		Alpha = 0;
		Sigma = 0;
		mu = .0025;
		Energy = 0;
		K_angle = 0.0;
		loop_count = 0;

		Energy_initialized = false;

		Vx = new double*[N_layers];
		Vy = new double*[N_layers];

		Gx = new double*[N_layers];
		Gy = new double*[N_layers];

		La = new double*[N_layers];
		Ka = new double*[N_layers];
		Ka_con = new double*[N_layers];
		A0 = new double*[N_layers];
		A0_con = new double*[N_layers];
		Atmp = new double*[N_layers];
		Kl = new double*[N_layers-1];
		L0 = new double*[N_layers-1];
		Ltmp = new double*[N_layers-1];

		for (int i=0;i<N_layers;++i){
			Vx[i] = new double[N];
			Vy[i] = new double[N];
			Gx[i] = new double[N];
			Gy[i] = new double[N];

			La[i] = new double[N];
			Ka[i] = new double[N];
			Ka_con[i] = new double[N];
			A0[i] = new double[N];
			A0_con[i] = new double[N];
			Atmp[i] = new double[N];
			if (i == (N_layers-1)) continue;
			Kl[i] = new double[N];
			L0[i] = new double[N];
			Ltmp[i] = new double[N];
		}

		Kv = new double[N];
		Vtmp = new double[N];
		V0 = new double[N];

		R = new double[N_layers];
		for (int i=0;i<N_layers;++i)
			R[i] = 14.0 - i*(14.0 - 8.9)/(N_layers - 1); 
		

		double theta, pi = 3.141592653589;
		// Initialize the position of the verticies
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				theta = -pi/2.0 + j*2.0*pi/N;
				Vx[i][j] = R[i]*cos(theta);
				Vy[i][j] = R[i]*sin(theta);
			}
		}

		// Initialize the initial geometrical parameters
		Y0 = 0.0;Ky = 0.0;
		for (int i=0;i<N;i++){
			V0[i] = 0.0; Kv[i] = 1.0;
		}

		double tmpx = 0, tmpy = 0; int jN = 0;
		for (int j=0;j<N;j++){
			jN = (j+1)%N;
			// Top layer
			tmpx = Vx[0][jN] - Vx[0][j];tmpy = Vy[0][jN] - Vy[0][j];
			V0[j] = (Vx[0][jN] + Vx[0][j])*tmpy - (Vy[0][jN] + Vy[0][j])*tmpx;
			A0[0][j] = sqrt(tmpx*tmpx + tmpy*tmpy);A0_con[0][j] = A0[0][j];

			// Bottom layer
			tmpx = Vx[N_layers-1][j] - Vx[N_layers-1][jN];tmpy = Vy[N_layers-1][j] - Vy[N_layers-1][jN];
			V0[j] += (Vx[N_layers-1][jN] + Vx[N_layers-1][j])*tmpy - (Vy[N_layers-1][jN] + Vy[N_layers-1][j])*tmpx;
			A0[N_layers-1][j] = sqrt(tmpx*tmpx + tmpy*tmpy);A0_con[N_layers-1][j] = A0[N_layers-1][j];

			Y0 -= (Vx[N_layers-1][jN] + Vx[N_layers-1][j])*tmpy/4.0 - (Vy[N_layers-1][jN] + Vy[N_layers-1][j])*tmpx/4.0;
	
			// Middle layer
			for (int k=0;k < N_layers-1;k++){
				tmpx = Vx[k+1][jN] - Vx[k+1][j];tmpy = Vy[k+1][jN] - Vy[k+1][j];
				A0[k+1][j] = sqrt(tmpx*tmpx + tmpy*tmpy);A0_con[k+1][j] = A0[k+1][j];

				tmpx = Vx[k+1][jN] - Vx[k][jN];tmpy = Vy[k+1][jN] - Vy[k][jN];
			
				V0[j] -= (Vx[k][jN] + Vx[k+1][jN])*tmpy - (Vy[k][jN] + Vy[k+1][jN])*tmpx;
				V0[jN] += (Vx[k][jN] + Vx[k+1][jN])*tmpy - (Vy[k][jN] + Vy[k+1][jN])*tmpx;

				L0[k][jN] = sqrt(pow(tmpx, 2) + pow(tmpy, 2));
			}
			V0[j] /= 4.0;
		}

		
		for (int j=0;j<N;j++){
			for (int k=0;k<N_layers;k++){
				Atmp[k][j] = A0[k][j];
				Ka[k][j] = 0.0;
				Ka_con[k][j] = 0.0;
				La[k][j] = 0.0;
				if (k == (N_layers-1)) continue;
					Kl[k][j] = 0.0;
					Ltmp[k][j] = L0[k][j];
			}
		}

	};

	void set_kangle(double Kangle_in){
		K_angle = Kangle_in;
	};

	void set_ky(double Ky_in){
		Ky = Ky_in;
	};

	void set_kv(double Kv_in){
		for (int i=0;i<N;i++)
			Kv[i] = Kv_in;
	};

	void set_vo(double vo_in){
		for (int i=0;i<N;i++)
			V0[i] = vo_in;
	};

	void set_ka(double Ka_in){
		for (int j=0;j<N;j++)
			for (int k=0;k<N_layers;k++)
				Ka[k][j] = Ka_in;
	};

	void set_ka(double Ka_in, int layer_id){
		for (int j=0;j<N;j++)
			Ka[layer_id][j] = Ka_in;
	};

	void set_kl(double Kl_in){
		for (int j=0;j<N;j++)
			for (int k=0;k<N_layers-1;k++)
				Kl[k][j] = Kl_in;
	};

	void set_kl(double Kl_in, int layer_id){
		for (int j=0;j<N;j++)
			Kl[layer_id][j] = Kl_in;
	};

	void set_lo(double lo_in){
		for (int j=0;j<N;j++)
			for (int k=0;k<N_layers-1;k++)
				L0[k][j] = lo_in;
	};

	void set_lo(double n_ecto_Lo, double K_ecto_Lo, double Lo_factor){
		for (int j=0;j<N;j++)
			for (int k=0;k<N_layers-1;k++)
				L0[k][j] *= ((1-Lo_factor)/(1+pow((j-N/2.0)/K_ecto_Lo, n_ecto_Lo)) + Lo_factor);
	};

	void set_ka_con(double Alpha_in, double Sigma_in){
		double theta, pi = 3.141592653589;
		Alpha = Alpha_in; Sigma = Sigma_in;
		for (int i=0;i<N;i++){
			for (int k=0;k<N_layers-1;k++){
				theta = (-pi*R[0] + i*2*pi*R[0]/(N-1))/Sigma;
				Ka_con[k][i] = Alpha*exp(-theta*theta/2.0);
			}
		}
	};

	void set_ka_con(double Alpha_in, double Sigma_in, int layer_id){
		double theta, pi = 3.141592653589;
		Alpha = Alpha_in; Sigma = Sigma_in;
		for (int i=0;i<N;i++){
			theta = (-pi*R[0] + i*2*pi*R[0]/(N-1))/Sigma;
			Ka_con[layer_id][i] = Alpha*exp(-theta*theta/2.0);	
		}
	};

	void set_a0_con(double K_meso_Ao_drop, double n_meso_Ao_drop){
		for (int i=0;i<N;i++)
			for (int k=0;k<N_layers;k++)
				A0_con[k][i] *= (1 - 1/(1+pow((i+1-(N+1)/2.0)/K_meso_Ao_drop, n_meso_Ao_drop)));	
	};

	void set_a0_con(double K_meso_Ao_drop, double n_meso_Ao_drop, int layer_id){
		for (int i=0;i<N;i++)
			A0_con[layer_id][i] *= (1 - 1/(1+pow((i+1-(N+1)/2.0)/K_meso_Ao_drop, n_meso_Ao_drop)));	
	};

	void set_la(double Alpha_in, double Sigma_in){
		double theta, pi = 3.141592653589;
		//Alpha = Alpha_in; Sigma = Sigma_in;
		for (int i=0;i<N;i++){
			for (int k=0;k<N_layers-1;k++){
				theta = (-pi*R[0] + i*2*pi*R[0]/(N-1))/Sigma_in;
				La[k][i] = Alpha_in*exp(-theta*theta/2.0);
			}
		}
	};

	void set_la(double Alpha_in, double Sigma_in, int layer_id){
		double theta, pi = 3.141592653589;
		//Alpha = Alpha_in; Sigma = Sigma_in;
		for (int i=0;i<N;i++){
			theta = (-pi*R[0] + i*2*pi*R[0]/(N-1))/Sigma_in;
			La[layer_id][i] = Alpha_in*exp(-theta*theta/2.0);	
		}
	};

	bool Evolve(int N_steps){
		double Y_tmp = 0, E_total_old = Energy, E_total;
		double tmp_dax, tmp_day, tmp_a, tmp_dlx, tmp_dly, tmp_l, tmp, tmp2, R_tmp, A_tmp, B_tmp, C_tmp;
		double m_a_f = 0.06;
		double E_o = 5000.0*Kv[0], K = 500.0;
		mu = .0015/Kv[0];
		int jN, jN2;

		for (int n_st=0;n_st < N_steps;++n_st){
			
			for (int j=0;j<N;j++){
				Vtmp[j] = 0.0;
				for (int k=0;k<N_layers;k++){
					Gx[k][j] = 0.0;Gy[k][j] = 0.0;
				}
			}
			Y_tmp = 0.0; E_total = 0.0;

			// Loop over all the verticies
			for (int j=0;j<N;j++){
				jN = (j+1)%N;
				// Top layer
				tmp_dax = Vx[0][jN] - Vx[0][j];tmp_day = Vy[0][jN] - Vy[0][j];
				tmp_a = sqrt(tmp_dax*tmp_dax + tmp_day*tmp_day);
			
				// Energy for apical side
				if (tmp_a < m_a_f*A0[0][j]){
					tmp = (2.0*Ka[0][j]*(tmp_a - A0[0][j]) + 2.0*Ka_con[0][j]*(tmp_a - A0_con[0][j]) + La[0][j]*tmp_a/(m_a_f*A0[0][j]))/tmp_a;
					E_total += Ka[0][j]*(tmp_a - A0[0][j])*(tmp_a - A0[0][j]) + Ka_con[0][j]*(tmp_a - A0_con[0][j])*(tmp_a - A0_con[0][j]) + La[0][j]*tmp_a*tmp_a/(m_a_f*A0[0][j]);
				}
				else{
					tmp = (2.0*Ka[0][j]*(tmp_a - A0[0][j]) + 2.0*Ka_con[0][j]*(tmp_a - A0_con[0][j]) + La[0][j])/tmp_a;
					E_total += Ka[0][j]*(tmp_a - A0[0][j])*(tmp_a - A0[0][j]) + Ka_con[0][j]*(tmp_a - A0_con[0][j])*(tmp_a - A0_con[0][j]) + La[0][j]*tmp_a;
				}
			
				Gx[0][jN] -= tmp*tmp_dax; Gy[0][jN] -= tmp*tmp_day;
				Gx[0][j] += tmp*tmp_dax; Gy[0][j] += tmp*tmp_day;
				// Sum up volume for the apical side
				Vtmp[j] += (Vx[0][jN] + Vx[0][j])*tmp_day/4.0 - (Vy[0][jN] + Vy[0][j])*tmp_dax/4.0;
		
				for (int k=0;k<N_layers-1;k++){
					// Bottom layer
					tmp_dax = Vx[k+1][jN] - Vx[k+1][j];tmp_day = Vy[k+1][jN] - Vy[k+1][j];
					tmp_a = sqrt(tmp_dax*tmp_dax + tmp_day*tmp_day);
					tmp = (2.0*Ka[k+1][j]*(tmp_a - A0[k+1][j]) + 2.0*Ka_con[k+1][j]*(tmp_a - A0_con[k+1][j]) + La[k+1][j])/tmp_a;
					// Energy for basal side
					E_total += Ka[k+1][j]*(tmp_a - A0[k+1][j])*(tmp_a - A0[k+1][j]) + Ka_con[k+1][j]*(tmp_a - A0_con[k+1][j])*(tmp_a - A0_con[k+1][j]) + La[k+1][j]*tmp_a;

					Gx[k+1][jN] -= tmp*tmp_dax; Gy[k+1][jN] -= tmp*tmp_day;
					Gx[k+1][j] += tmp*tmp_dax; Gy[k+1][j] += tmp*tmp_day;
					// Sum up volume for the basal side, the yolk too
					if (k == (N_layers-2)){
						tmp = (Vx[k+1][jN] + Vx[k+1][j])*tmp_day/4.0 - (Vy[k+1][jN] + Vy[k+1][j])*tmp_dax/4.0;
						Vtmp[j] -= tmp;
						Y_tmp += tmp;
					}
					// Lateral layer
					tmp_dlx = Vx[k+1][jN] - Vx[k][jN];tmp_dly = Vy[k+1][jN] - Vy[k][jN];
					tmp_l = sqrt(tmp_dlx*tmp_dlx + tmp_dly*tmp_dly);
					tmp = (Vx[k][jN] + Vx[k+1][jN])*tmp_dly/4.0 - (Vy[k][jN] + Vy[k+1][jN])*tmp_dlx/4.0;
					// Sum up volume for the right side
					Vtmp[jN] -= tmp;
					// Sum up volume for the left side
					Vtmp[j] += tmp;
					tmp = 2.0*Kl[k][j]*(tmp_l - L0[k][j])/tmp_l;
					// Energy for lateral side
					E_total += Kl[k][jN]*(tmp_l - L0[k][jN])*(tmp_l - L0[k][jN]);
			
					Gx[k][jN] += tmp*tmp_dlx; Gy[k][jN] += tmp*tmp_dly;
					Gx[k+1][jN] -= tmp*tmp_dlx; Gy[k+1][jN] -= tmp*tmp_dly;

					// Apply the membrane force
					if (j>45 && j<65){
						tmp = sqrt((Vx[k+1][j+1] - Vx[k][j+1])*(Vx[k+1][j+1] - Vx[k][j+1]) + (Vy[k+1][j+1] - Vy[k][j+1])*(Vy[k+1][j+1] - Vy[k][j+1]));
						R_tmp = (-(Vy[k+1][j+1] - Vy[k][j+1])*(Vx[k+1][j] - Vx[k][j+1]) + (Vy[k+1][j] - Vy[k][j+1])*(Vx[k+1][j+1] - Vx[k][j+1]))/tmp;
						E_total += E_o/(1.0 + exp(K*R_tmp));
						tmp2 = -E_o*K*exp(K*R_tmp)/(1.0+exp(K*R_tmp))/(1.0+exp(K*R_tmp));
						Gx[k+1][j] += tmp2*(Vy[k+1][j+1] - Vy[k][j+1])/tmp;
						Gy[k+1][j] -= tmp2*(Vx[k+1][j+1] - Vx[k][j+1])/tmp;
					}

					if (j>15 && j<35){
						tmp = sqrt((Vx[k+1][j-1] - Vx[k][j-1])*(Vx[k+1][j-1] - Vx[k][j-1]) + (Vy[k+1][j-1] - Vy[k][j-1])*(Vy[k+1][j-1] - Vy[k][j-1]));
						R_tmp = (-(Vy[k+1][j-1] - Vy[k][j-1])*(Vx[k+1][j] - Vx[k][j-1]) + (Vy[k+1][j] - Vy[k][j-1])*(Vx[k+1][j-1] - Vx[k][j-1]))/tmp;
						E_total += E_o/(1.0 + exp(K*R_tmp));
						tmp2 = E_o*K*exp(K*R_tmp)/(1.0+exp(K*R_tmp))/(1.0+exp(K*R_tmp));
						Gx[k+1][j] += tmp2*(Vy[k+1][j-1] - Vy[k][j-1])/tmp;
						Gy[k+1][j] -= tmp2*(Vx[k+1][j-1] - Vx[k][j-1])/tmp;
					}
					
				}	
			}

			// Introduce Angles
			if (N_layers > 2){
				for (int j=0;j<N;++j){
					for (int k=1;k<N_layers-1;k++){
						A_tmp = sqrt((Vx[k-1][j] - Vx[k][j])*(Vx[k-1][j] - Vx[k][j]) + (Vy[k-1][j] - Vy[k][j])*(Vy[k-1][j] - Vy[k][j]));
						B_tmp = sqrt((Vx[k+1][j] - Vx[k][j])*(Vx[k+1][j] - Vx[k][j]) + (Vy[k+1][j] - Vy[k][j])*(Vy[k+1][j] - Vy[k][j]));
						C_tmp = ((Vx[k-1][j] - Vx[k][j])*(Vy[k+1][j] - Vy[k][j]) - (Vx[k+1][j] - Vx[k][j])*(Vy[k-1][j] - Vy[k][j]));
						tmp = C_tmp/A_tmp/B_tmp;
						E_total += K_angle*tmp*tmp;

						Gx[k][j] -= 2.0*K_angle*tmp*((Vy[k-1][j] - Vy[k+1][j]) + C_tmp*((Vx[k-1][j] - Vx[k][j])/A_tmp/A_tmp + (Vx[k+1][j] - Vx[k][j])/B_tmp/B_tmp))/A_tmp/B_tmp;
						Gy[k][j] -= 2.0*K_angle*tmp*((Vx[k+1][j] - Vx[k-1][j]) + C_tmp*((Vy[k-1][j] - Vy[k][j])/A_tmp/A_tmp + (Vy[k+1][j] - Vy[k][j])/B_tmp/B_tmp))/A_tmp/B_tmp;
					}
				}
			}

			// Energy for the Yolk
			E_total += Ky*(Y_tmp - Y0)*(Y_tmp - Y0);
			for (int j=0;j<N;j++){
				jN = (j+1)%N;jN2 = (j+2)%N;
				tmp = Kv[j]*(Vtmp[j] - V0[j]);tmp2 = Kv[jN]*(Vtmp[jN] - V0[jN]);
				// Energy for volume of cell
				E_total += Kv[j]*(Vtmp[j] - V0[j])*(Vtmp[j] - V0[j]);
		
				// Gradient for volume of cell
				Gx[0][jN] -= tmp*(Vy[1][jN] - Vy[0][j]) + tmp2*(Vy[0][jN2] - Vy[1][jN]);
				Gy[0][jN] -= tmp*(Vx[0][j] - Vx[1][jN]) + tmp2*(Vx[1][jN] - Vx[0][jN2]);

				Gx[N_layers-1][jN] -= tmp*(Vy[N_layers-1][j] - Vy[N_layers-2][jN]) + tmp2*(Vy[N_layers-2][jN] - Vy[N_layers-1][jN2]);
				Gy[N_layers-1][jN] -= tmp*(Vx[N_layers-2][jN] - Vx[N_layers-1][j]) + tmp2*(Vx[N_layers-1][jN2] - Vx[N_layers-2][jN]);

				for (int k=1;k<N_layers-1;k++){
					Gx[k][jN] -= (tmp - tmp2)*(Vy[k+1][jN] - Vy[k-1][jN]);
					Gy[k][jN] -= (tmp - tmp2)*(Vx[k-1][jN] - Vx[k+1][jN]);
				}
				// Gradient for the yolk
				Gx[N_layers-1][jN] -= Ky*(Y_tmp - Y0)*(Vy[N_layers-1][jN2] - Vy[N_layers-1][j]);
				Gy[N_layers-1][jN] -= Ky*(Y_tmp - Y0)*(Vx[N_layers-1][j] - Vx[N_layers-1][jN2]);
			}

			
			if (!Energy_initialized){
				Energy_initialized = true;
				Energy = E_total;
				E_total_old = E_total;
			}
			
				

			if (E_total_old < E_total){
				Energy = E_total_old;
				return false;
			}
			else
				E_total_old = E_total;
			
			

			// Move the verticies
			for (int j=0;j<N;j++){
				for (int k=0;k<N_layers;k++){
					Vx[k][j] += Gx[k][j]*mu; Vy[k][j] += Gy[k][j]*mu;
				}
			}

			
			// Apply Boundary Conditions
			for (int j=0;j<N;++j){
				R_tmp = sqrt(Vx[0][j]*Vx[0][j] + Vy[0][j]*Vy[0][j]);
				if (R_tmp > R[0]){
					Vx[0][j] *= R[0]/R_tmp;Vy[0][j] *= R[0]/R_tmp;}
			
				if (j<=39 && Vx[0][j] < 0)
					Vx[0][j] = 0.0;
				else if (j>=40 && Vx[0][j] > 0)
					Vx[0][j] = 0.0;
			}
			
		}
		Energy = E_total;
		loop_count += N_steps;
		return true;
	};

	void create_output_file(string filename_in){
		
		int st_size = filename_in.length(); 
		filename = new char[st_size];
		for (int i=0;i<st_size;++i)
			filename[i] = filename_in[i];

		filename[st_size] = '\0';
		// Record everything
		ofstream* file = new ofstream(filename);
		
		(*file)<<"i=1;"<<endl;
		(*file)<<"Vx{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<Vx[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"Vy{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<Vy[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"Gx{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<0<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"Gy{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<0<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;
		
		(*file)<<"P.A{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<A0[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.L{i} = [";
		for (int i=0;i<N_layers-1;i++){
			for (int j=0;j<N;j++){
				(*file)<<L0[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.V{i} = [";
		for (int j=0;j<N;j++){
			(*file)<<V0[j]<<" ";}
		(*file)<<";];"<<endl;
		
		(*file)<<"P.Y(i) = ["<<Y0<<"];"<<endl;

		(*file)<<"P.KA = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<Ka[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.KA_con = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<Ka_con[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.A0_con = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<A0_con[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.LA = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				(*file)<<La[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.KL = [";
		for (int i=0;i<N_layers-1;i++){
			for (int j=0;j<N;j++){
				(*file)<<Kl[i][j]<<" ";}
			if (i != N_layers-1)
				(*file)<<";";
		}
		(*file)<<"];"<<endl;

		(*file)<<"P.KY = ["<< Ky << "];"<<endl;

		(*file)<<"E(i) = "<<Energy << ";" << endl;

		file->close();
		delete file;

	};

	void output(){

		// Record the positions and other stuff
		//cout << loop_count << " steps completed" << endl;
		
		double tmpx = 0, tmpy = 0; int jN = 0;
		for (int j=0;j<N;j++){
			jN = (j+1)%N;
			// Top layer
			tmpx = Vx[0][jN] - Vx[0][j];tmpy = Vy[0][jN] - Vy[0][j];
			Atmp[0][j] = sqrt(tmpx*tmpx + tmpy*tmpy);

			// Bottom layer
			tmpx = Vx[N_layers-1][j] - Vx[N_layers-1][jN];tmpy = Vy[N_layers-1][j] - Vy[N_layers-1][jN];
			Atmp[N_layers-1][j] = sqrt(tmpx*tmpx + tmpy*tmpy);
	
			// Middle layer
			for (int k=0;k < N_layers-1;k++){
				tmpx = Vx[k+1][jN] - Vx[k+1][j];tmpy = Vy[k+1][jN] - Vy[k+1][j];
				Atmp[k+1][j] = sqrt(tmpx*tmpx + tmpy*tmpy);

				tmpx = Vx[k+1][jN] - Vx[k][jN];tmpy = Vy[k+1][jN] - Vy[k][jN];
				Ltmp[k][jN] = sqrt(pow(tmpx, 2) + pow(tmpy, 2));
			}
		}

		ofstream file;
		file.open(filename, ios::app);
		
		file<<"i=i+1;"<<endl;
		file<<"Vx{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				file<<Vx[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;

		file<<"Vy{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				file<<Vy[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;

		file<<"Gx{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				file<<Gx[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;

		file<<"Gy{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				file<<Gy[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;
		
		file<<"P.A{i} = [";
		for (int i=0;i<N_layers;i++){
			for (int j=0;j<N;j++){
				file<<Atmp[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;

		file<<"P.L{i} = [";
		for (int i=0;i<N_layers-1;i++){
			for (int j=0;j<N;j++){
				file<<Ltmp[i][j]<<" ";}
			if (i!= N_layers-1)
				file<<";";
		}
		file<<"];"<<endl;

		file<<"P.V{i} = [";
		for (int j=0;j<N;j++){
			file<<Vtmp[j]<<" ";}
		file<<";];"<<endl;

		file<<"P.Y(i) = ["<<Y0<<"];"<<endl;
		file<<"E(i) = "<<Energy << ";" << endl;
		
		file.close();

	};

	void draw(){

		double scale = 10;
		double outer_rad = R[0]*scale;

		int x_off = 150, y_off = 150;

		BCX_Circle(hConWnd, x_off, y_off, outer_rad*1.1, 1, Red);
		//BCX_Circle(hConWnd, x_off, y_off, outer_rad, Red);

		int jN;

		for (int k=0;k<N_layers;k++){
			for (int j=0;j<N;j++){
				jN = (j+1)%N;
				BCX_Line(hConWnd, Vx[k][j]*scale + x_off, -Vy[k][j]*scale + y_off, Vx[k][jN]*scale + x_off, -Vy[k][jN]*scale + y_off, Red);
			}
			if (k==N_layers-1) continue;
			for (int j=0;j<N;++j){
				BCX_Line(hConWnd, Vx[k][j]*scale + x_off, -Vy[k][j]*scale + y_off, Vx[k+1][j]*scale + x_off, -Vy[k+1][j]*scale + y_off, Red);
			}
		}


	};

};

