diff --git a/src/EoS.h b/src/EoS.h
index 2284bb9..cbf143f 100644
--- a/src/EoS.h
+++ b/src/EoS.h
@@ -1,15 +1,15 @@
 #pragma once
 
-namespace EoS{
+namespace EoS {
 	//////////////////
 	//abstract class
 	//////////////////
-	template <typename type> class EoS_t{
-		public:
-		EoS_t(){
+	template <typename type> class EoS_t {
+	  public:
+		EoS_t() {
 			return ;
 		}
-		/* virtual */ ~EoS_t(){
+		/* virtual */ ~EoS_t() {
 			return ;
 		}
 		virtual type Pressure  (const type dens, const type eng) const = 0;
@@ -18,75 +18,75 @@ namespace EoS{
 	//////////////////
 	//EoSs
 	//////////////////
-	template <typename type> class IdealGas : public EoS_t<type>{
+	template <typename type> class IdealGas : public EoS_t<type> {
 		const type hcr;//heat capacity ratio;
-		public:
-		IdealGas(const type _hcr) : hcr(_hcr){
+	  public:
+		IdealGas(const type _hcr) : hcr(_hcr) {
 		}
-		inline type Pressure(const type dens, const type eng) const{
+		inline type Pressure(const type dens, const type eng) const {
 			return (hcr - 1.0) * dens * eng;
 		}
-		inline type SoundSpeed(const type dens, const type eng) const{
+		inline type SoundSpeed(const type dens, const type eng) const {
 			return sqrt(hcr * (hcr - 1.0) * eng);
 		}
-		inline type HeatCapacityRatio() const{
+		inline type HeatCapacityRatio() const {
 			return hcr;
 		}
 	};
 
-    //Tillotson equation "is applicable to the prediction of the shock and release of materials undergoing hypervelocity impacts."
-	template <typename type> class Tillotson : public EoS_t<type>{
+	//Tillotson equation "is applicable to the prediction of the shock and release of materials undergoing hypervelocity impacts."
+	template <typename type> class Tillotson : public EoS_t<type> {
 		type rho0, a, b, A, B, u0, alpha, beta, uiv, ucv;
-		inline type P_co(const type dens, const type eng) const{
+		inline type P_co(const type dens, const type eng) const {
 			const type eta = dens / rho0;
 			const type mu  = eta - 1.0;
 			return (a + b / (eng / u0 / eta / eta + 1.0)) * dens * eng + A * mu + B * mu * mu;
 		}
-		inline type P_ex(const type dens, const type eng) const{
+		inline type P_ex(const type dens, const type eng) const {
 			const type eta = dens / rho0;
 			const type mu  = eta - 1.0;
 			return a * dens * eng + (b * dens * eng / (eng / u0 / eta / eta + 1.0) + A * mu * exp(- alpha * (1.0 / eta - 1.0))) * exp(- beta * (1.0 / eta - 1.0) * (1.0 / eta - 1.0));
 		}
-		inline type dPdrho(const type rho, const type u) const{
+		inline type dPdrho(const type rho, const type u) const {
 			const type drho = 0.0001;
 			return (Pressure(rho + drho, u) - Pressure(rho - drho, u)) / (2.0 * drho);
 		}
-		inline type dPdu(const type rho, const type u) const{
+		inline type dPdu(const type rho, const type u) const {
 			const type du = 0.0001;
 			return (Pressure(rho, u + du) - Pressure(rho, u - du)) / (2.0 * du);
 		}
-		public:
-		Tillotson(const type a_rho0, const type a_u0, const type a_uiv, const type a_ucv, const type a_A, const type a_B, const type a_a, const type a_b, const type a_alpha, const type a_beta){
-            //in MKS unit...
-            //From Brundage 2013
-            //"Implementation of Tillotson Equation of State for Hypervelocity Impact of Metals, Geologic Materials, and Liquids"
-            rho0  = a_rho0; // Density                          [kg/m^3]
-            u0    = a_u0;   // Initial Energy                   [J/kg]
-            uiv   = a_uiv;  // Energy at incipient vaporization [J/kg]
-            ucv   = a_ucv;  // Energy at complete vaporization  [J/kg]
-            A     = a_A;    // Bulk modulus                     [Pa]
-            B     = a_B;    // Tillotson parameter              [Pa]
-            a     = a_a;    // Tillotson parameter              [dimension-less]
-            b     = a_b;    // Tillotson parameter
-            alpha = a_alpha;// Tillotson parameter
-            beta  = a_beta; // Tillotson parameter
+	  public:
+		Tillotson(const type a_rho0, const type a_u0, const type a_uiv, const type a_ucv, const type a_A, const type a_B, const type a_a, const type a_b, const type a_alpha, const type a_beta) {
+			//in MKS unit...
+			//From Brundage 2013
+			//"Implementation of Tillotson Equation of State for Hypervelocity Impact of Metals, Geologic Materials, and Liquids"
+			rho0  = a_rho0; // Density                          [kg/m^3]
+			u0    = a_u0;   // Initial Energy                   [J/kg]
+			uiv   = a_uiv;  // Energy at incipient vaporization [J/kg]
+			ucv   = a_ucv;  // Energy at complete vaporization  [J/kg]
+			A     = a_A;    // Bulk modulus                     [Pa]
+			B     = a_B;    // Tillotson parameter              [Pa]
+			a     = a_a;    // Tillotson parameter              [dimension-less]
+			b     = a_b;    // Tillotson parameter
+			alpha = a_alpha;// Tillotson parameter
+			beta  = a_beta; // Tillotson parameter
 		}
-		inline type Pressure(const type dens, const type eng) const{
+		inline type Pressure(const type dens, const type eng) const {
 			const type p_min = 1.0e+7;
-			if(dens >= rho0 || eng < uiv){
+			if(dens >= rho0 || eng < uiv) {
 				return std::max(P_co(dens, eng), p_min);
-			}else if(dens < rho0 && eng > ucv){
+			} else if(dens < rho0 && eng > ucv) {
 				return std::max(P_ex(dens, eng), p_min);
-			}else{
+			} else {
 				return std::max(((eng - uiv) * P_ex(dens, eng) + (ucv - eng) * P_co(dens, eng)) / (ucv - uiv), p_min);
 			}
 		}
-		inline type SoundSpeed(const type dens, const type eng) const{
+		inline type SoundSpeed(const type dens, const type eng) const {
 			return sqrt(std::max(Pressure(dens, eng) / (dens * dens) * dPdu(dens, eng) + dPdrho(dens, eng), 0.0) + 1.0e-16);
 		}
 	};
 
-	template <typename type> class ANEOS : public EoS_t<type>{
+	template <typename type> class ANEOS : public EoS_t<type> {
 		/**
 		 * A member variable that contains all the EoS data from ANEOS. The lines
 		 * represent density, the columns energy. The data values for each entry
@@ -111,7 +111,7 @@ namespace EoS{
 		 * The function returns the table index of the next larger density and energy.
 		 */
 		std::pair<unsigned int, unsigned int>
-		get_table_index(const type density, const type energy) const{
+		get_table_index(const type density, const type energy) const {
 			auto line = std::lower_bound(densities.begin(),densities.end(),density);
 			const unsigned int line_index = std::distance(densities.begin(),line);
 			auto column = std::lower_bound(energies.begin(),energies.end(),energy);
@@ -125,45 +125,42 @@ namespace EoS{
 
 		const type
 		get_interpolated_value (const type density,
-				const type energy,
-				const unsigned int property_index) const{
+		                        const type energy,
+		                        const unsigned int property_index) const {
 			const std::pair<unsigned int, unsigned int> data_index = get_table_index(density, energy);
 
 			// If we are not at the boundaries of the data table
 			if ((data_index.first < densities.size()-1) &&
-					(data_index.second < energies.size()-1))
-			{
+			        (data_index.second < energies.size()-1)) {
 				// compute the interpolation weights of this density and energy
 				const type xi = (density - densities[data_index.first-1]) /
-						(densities[data_index.first] - densities[data_index.first-1]);
+				                (densities[data_index.first] - densities[data_index.first-1]);
 
 				const type eta = (energy - energies[data_index.second-1]) /
-						(energies[data_index.second] - energies[data_index.second-1]);
+				                 (energies[data_index.second] - energies[data_index.second-1]);
 
 				// use these coordinates for a bilinear interpolation
 				return  (1-xi)*(1-eta) * eos_data[data_index.first-1][data_index.second-1][property_index] +
-						xi    *(1-eta) * eos_data[data_index.first]  [data_index.second-1][property_index] +
-						(1-xi)*eta     * eos_data[data_index.first-1][data_index.second]  [property_index] +
-						xi    *eta     * eos_data[data_index.first]  [data_index.second]  [property_index];
-			}
-			else
-			{
+				        xi    *(1-eta) * eos_data[data_index.first]  [data_index.second-1][property_index] +
+				        (1-xi)*eta     * eos_data[data_index.first-1][data_index.second]  [property_index] +
+				        xi    *eta     * eos_data[data_index.first]  [data_index.second]  [property_index];
+			} else {
 				// Return the boundary value
 				return eos_data[data_index.first][data_index.second][property_index];
 			}
 			return eos_data[data_index.first][data_index.second][property_index];
 		}
 
-		public:
+	  public:
 		/**
 		 * Construct the EoS by reading its data from a file. The file format looks
 		 * as follows:
 		 * TODO fill in format description, check that my implementation is correct
 		 */
-		ANEOS(const std::string &filename){
+		ANEOS(const std::string &filename) {
 			std::ifstream input;
 			input.open(filename, std::ios::in);
-			if(input.fail() == true){
+			if(input.fail() == true) {
 				std::cout << "Cannot open file..." << filename << std::endl;
 				return;
 			}
@@ -176,9 +173,9 @@ namespace EoS{
 			unsigned int column_index = 0;
 			std::string line;
 
-			while(std::getline(input, line)){
+			while(std::getline(input, line)) {
 				// read the table size, which has to be of the format # n_rows n_cols
-				if (n_lines == 0 && n_expected_lines == 0 && n_expected_columns == 0){
+				if (n_lines == 0 && n_expected_lines == 0 && n_expected_columns == 0) {
 					std::string tmp;
 					std::istringstream stream(line);
 
@@ -196,11 +193,11 @@ namespace EoS{
 					continue;
 
 				// if we have reached the expected number of columns add a new row
-				if(column_index == n_expected_columns){
+				if(column_index == n_expected_columns) {
 					++n_lines;
 					column_index = 0;
 					eos_data.push_back(std::vector<std::array<type,6> >());
-					}
+				}
 
 				const unsigned int line_index = n_lines - 1;
 				eos_data[line_index].push_back(std::array<type,6>());
@@ -210,8 +207,7 @@ namespace EoS{
 				unsigned int field_index = 0;
 				type tmp;
 
-				while(stream >> tmp)
-				{
+				while(stream >> tmp) {
 					eos_data[line_index][column_index][field_index] = tmp;
 
 					if (column_index == 0 && field_index == 0)
@@ -229,20 +225,20 @@ namespace EoS{
 			// test_data();
 		}
 
-		inline type Pressure(const type dens, const type eng) const{
+		inline type Pressure(const type dens, const type eng) const {
 			return get_interpolated_value(dens,eng,3);
 		}
 
-		inline type SoundSpeed(const type dens, const type eng) const{
+		inline type SoundSpeed(const type dens, const type eng) const {
 			return get_interpolated_value(dens,eng,4);
 		}
 
-		void test_data() const{
+		void test_data() const {
 			std::ofstream output;
 			output.open("test_output.txt", std::ios::out);
 			output.setf(std::ios::scientific);
 			output.precision(8);
-			if(output.fail() == true){
+			if(output.fail() == true) {
 				std::cout << "Cannot open file test_output.txt" << std::endl;
 				return;
 			}
@@ -252,11 +248,9 @@ namespace EoS{
 			unsigned int n_expected_columns = energies.size();
 
 			for (unsigned int line = 0; line < n_expected_lines; ++line)
-				for (unsigned int column = 0; column < n_expected_columns; ++column)
-				{
+				for (unsigned int column = 0; column < n_expected_columns; ++column) {
 					output << ' ';
-					for (unsigned int field = 0; field < n_expected_fields; ++field)
-					{
+					for (unsigned int field = 0; field < n_expected_fields; ++field) {
 						output << eos_data[line][column][field] << "  ";
 					}
 
@@ -270,6 +264,6 @@ namespace EoS{
 static const EoS::IdealGas<PS::F64>  Monoatomic(5./3.);
 static const EoS::IdealGas<PS::F64>  Diatomic  (1.4);
 static const EoS::Tillotson<PS::F64> Granite   (2680.0, 16.0e+6, 3.5e+6, 18.00e+6,  18.0e+9,  18.0e+9, 0.5, 1.3, 5.0, 5.0);
-static const EoS::Tillotson<PS::F64> Iron      (7800.0,  9.5e+6, 2.4e+6 , 8.67e+6, 128.0e+9, 105.0e+9, 0.5, 1.5, 5.0, 5.0);
+static const EoS::Tillotson<PS::F64> Iron      (7800.0,  9.5e+6, 2.4e+6, 8.67e+6, 128.0e+9, 105.0e+9, 0.5, 1.5, 5.0, 5.0);
 static const EoS::ANEOS<PS::F64> AGranite      ("eos/granite.rho_u.txt");
 
diff --git a/src/GI.h b/src/GI.h
index 921f355..90aa18a 100644
--- a/src/GI.h
+++ b/src/GI.h
@@ -5,12 +5,12 @@
 #ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 #error
 #endif
-template <class Ptcl> class GI : public Problem<Ptcl>{
-    public:
+template <class Ptcl> class GI : public Problem<Ptcl> {
+  public:
 	static double end_time;
-    static double damping;
+	static double damping;
 	static void setupIC(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo, PS::DomainInfo& dinfo,
-                        const std::string &input_file){
+	                    const std::string &input_file) {
 		const double Corr = .98;//Correction Term
 		/////////
 		//place ptcls
@@ -21,18 +21,18 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		/////////
 
 		// Use parameters from input file, or defaults if none provided
-        ParameterFile parameter_file(input_file);
-        std::cout << "Reading parameters from " << input_file << std::endl;
+		ParameterFile parameter_file(input_file);
+		std::cout << "Reading parameters from " << input_file << std::endl;
 		PS::F64 UnitMass = parameter_file.getValueOf("UnitMass", 6.0e+24);
 		PS::F64 UnitRadi = parameter_file.getValueOf("UnitRadi", 6400e+3);
 		PS::F64 coreFracRadi = parameter_file.getValueOf("coreFracRadi", 3500.0e+3 / 6400.0e+3);
 		PS::F64 coreFracMass = parameter_file.getValueOf("coreFracMass", 0.3);
 		PS::F64 imptarMassRatio = parameter_file.getValueOf("imptarMassRatio", 0.1);
-        int mode = parameter_file.getValueOf("mode", 2 );
-        PS::F64 impVel = parameter_file.getValueOf("impVel",0.);
-        end_time = parameter_file.getValueOf("end_time",1.0e+4);
-        damping = parameter_file.getValueOf("damping",1.);
-        
+		int mode = parameter_file.getValueOf("mode", 2 );
+		PS::F64 impVel = parameter_file.getValueOf("impVel",0.);
+		end_time = parameter_file.getValueOf("end_time",1.0e+4);
+		damping = parameter_file.getValueOf("damping",1.);
+
 		const PS::F64 Expand = 1.1;
 		const PS::F64 tarMass = UnitMass;
 		const PS::F64 tarRadi = UnitRadi;
@@ -53,9 +53,9 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		///////////////////
 		//target
 		int tarNmntl = 0;
-		for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-			for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-				for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
+		for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+			for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+				for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
 					const PS::F64 r = sqrt(x*x + y*y + z*z) * UnitRadi;
 					if(r >= tarRadi || r <= tarCoreRadi) continue;
 					++ tarNmntl;
@@ -64,11 +64,11 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		}
 		int tarNcore;
 		double tarCoreShrinkFactor = 1.0;
-		while(tarCoreShrinkFactor *= 0.99){
+		while(tarCoreShrinkFactor *= 0.99) {
 			tarNcore = 0;
-			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
 						const PS::F64 r = tarCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
 						if(r >= Corr * tarCoreRadi) continue;
 						++ tarNcore;
@@ -79,9 +79,9 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		}
 		//imp
 		int impNmntl = 0;
-		for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-			for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-				for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
+		for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+			for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+				for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
 					const PS::F64 r = Expand * sqrt(x*x + y*y + z*z) * UnitRadi;
 					if(r >= impRadi || r <= impCoreRadi) continue;
 					++ impNmntl;
@@ -90,11 +90,11 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		}
 		double impCoreShrinkFactor = 1.0;
 		int impNcore;
-		while(impCoreShrinkFactor *= 0.99){
+		while(impCoreShrinkFactor *= 0.99) {
 			impNcore = 0;
-			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
 						const PS::F64 r = Expand * impCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
 						if(r >= Corr * impCoreRadi) continue;
 						++ impNcore;
@@ -125,7 +125,7 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		std::cout << "    core density     : " << impCoreMass / (4.0 * math::pi / 3.0 * impCoreRadi * impCoreRadi * impCoreRadi * Corr * Corr * Corr) << std::endl;
 		std::cout << "    mantle density   : " << (impMass - impCoreMass) / (4.0 * math::pi / 3.0 * (impRadi * impRadi * impRadi - impCoreRadi * impCoreRadi * impCoreRadi)) << std::endl;
 		std::cout << "    mean density     : " << impMass / (4.0 * math::pi / 3.0 * impRadi * impRadi * impRadi) << std::endl;
-        std::cout << "    velocity         : " << impVel << std::endl;
+		std::cout << "    velocity         : " << impVel << std::endl;
 		std::cout << "Total:" << Nptcl << std::endl;
 		std::cout << "Tar-to-Imp mass ratio: " << (double)(impNmntl) / (double)(tarNmntl) << std::endl;
 		const int NptclIn1Node = Nptcl / PS::Comm::getNumberOfProc();
@@ -133,153 +133,153 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		//Real put
 		///////////////////
 		PS::S32 id = 0;
-        
-        switch (mode){
-            case 1:
-                std::cout << "creating target from tar.dat" << std::endl;
-                FILE * tarFile;
-                tarFile = fopen("input/tar.dat","r");
-                FileHeader tarheader;
-                int nptcltar;
-                nptcltar = tarheader.readAscii(tarFile);
-                std::cout << "num tar ptcl: " << nptcltar << std::endl;
-                for(int i=0; i<nptcltar; i++){
-                    Ptcl ith;
-                    ith.readAscii(tarFile);
-                    if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
-                }
-                for(PS::U32 i = 0 ; i < tar.size() ; ++ i){
-                    tar[i].mass /= (PS::F64)(Nptcl);
-                }
-                for(PS::U32 i = 0 ; i < tar.size() ; ++ i){
-                    ptcl.push_back(tar[i]);
-                }
-                
-                std::cout << "creating impactor from imp.dat" << std::endl;
-                FILE * impFile;
-                impFile = fopen("input/imp.dat","r");
-                FileHeader impheader;
-                int nptclimp;
-                nptclimp = impheader.readAscii(impFile);
-                std::cout << "num imp ptcl: " << nptclimp << std::endl;
-                for(int i=0; i<nptclimp; i++){
-                    Ptcl ith;
-                    ith.readAscii(impFile);
-                    ith.vel.x += (-1) * impVel;
-                    if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
-                }
-                for(PS::U32 i = 0 ; i < imp.size() ; ++ i){
-                    imp[i].mass /= (PS::F64)(Nptcl);
-                }
-                for(PS::U32 i = 0 ; i < imp.size() ; ++ i){
-                    ptcl.push_back(imp[i]);
-                }
-                break;
-            case 2:
-                //Put Tar.
-                std::cout << "creating target" << std::endl;
-                for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-                    for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-                        for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
-                            const PS::F64 r = sqrt(x*x + y*y + z*z) * UnitRadi;
-                            if(r >= tarRadi || r <= tarCoreRadi) continue;
-                            Ptcl ith;
-                            ith.pos.x = UnitRadi * x;
-                            ith.pos.y = UnitRadi * y;
-                            ith.pos.z = UnitRadi * z;
-                            ith.dens = (tarMass - tarCoreMass) / (4.0 / 3.0 * math::pi * (tarRadi * tarRadi * tarRadi - tarCoreRadi * tarCoreRadi * tarCoreRadi));
-                            ith.mass = tarMass + impMass;
-                            ith.eng  = 0.1 * Grav * tarMass / tarRadi;
-                            ith.id   = id++;
-                            // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&AGranite);
-                            ith.tag = 0;
-                            if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
-                        }
-                    }
-                }
-                for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-                    for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-                        for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
-                            const PS::F64 r = tarCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
-                            if(r >= Corr * tarCoreRadi) continue;
-                            Ptcl ith;
-                            ith.pos.x = tarCoreShrinkFactor * UnitRadi * x;
-                            ith.pos.y = tarCoreShrinkFactor * UnitRadi * y;
-                            ith.pos.z = tarCoreShrinkFactor * UnitRadi * z;
-                            ith.dens = tarCoreMass / (4.0 / 3.0 * math::pi * tarCoreRadi * tarCoreRadi * tarCoreRadi * Corr * Corr * Corr);
-                            ith.mass = tarMass + impMass;
-                            ith.eng  = 0.1 * Grav * tarMass / tarRadi;
-                            ith.id   = id++;
-                            // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&Iron);
-                            ith.tag = 1;
-                            if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
-                        }
-                    }
-                }
-                for(PS::U32 i = 0 ; i < tar.size() ; ++ i){
-                    tar[i].mass /= (PS::F64)(Nptcl);
-                }
-                for(PS::U32 i = 0 ; i < tar.size() ; ++ i){
-                    ptcl.push_back(tar[i]);
-                }
-                break;
-            case 3:
-                //imp
-                std::cout << "creating impactor" << std::endl;
-                for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-                    for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-                        for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
-                            const PS::F64 r = Expand * sqrt(x*x + y*y + z*z) * UnitRadi;
-                            if(r >= impRadi || r <= impCoreRadi) continue;
-                            Ptcl ith;
-                            ith.pos.x = Expand * UnitRadi * x + offset;
-                            ith.pos.y = Expand * UnitRadi * y;
-                            ith.pos.z = Expand * UnitRadi * z;
-                            ith.dens = (impMass - impCoreMass) / (4.0 / 3.0 * math::pi * (impRadi * impRadi * impRadi - impCoreRadi * impCoreRadi * impCoreRadi));
-                            ith.mass = tarMass + impMass;
-                            ith.eng  = 0.1 * Grav * tarMass / tarRadi;
-                            ith.id   = id++;
-                            // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&AGranite);
-                            ith.tag = 2;
-                            if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
-                        }
-                    }
-                }
-                for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx){
-                    for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx){
-                        for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx){
-                            const PS::F64 r = Expand * impCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
-                            if(r >= impCoreRadi) continue;
-                            Ptcl ith;
-                            ith.pos.x = Expand * impCoreShrinkFactor * UnitRadi * x + offset;
-                            ith.pos.y = Expand * impCoreShrinkFactor * UnitRadi * y;
-                            ith.pos.z = Expand * impCoreShrinkFactor * UnitRadi * z;
-                            ith.dens = impCoreMass / (4.0 / 3.0 * math::pi * impCoreRadi * impCoreRadi * impCoreRadi * Corr * Corr * Corr);
-                            ith.mass = tarMass + impMass;
-                            ith.eng  = 0.1 * Grav * tarMass / tarRadi;
-                            ith.id   = id++;
-                            // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&Iron);
-                            ith.tag = 3;
-                            if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
-                        }
-                    }
-                }
-                for(PS::U32 i = 0 ; i < imp.size() ; ++ i){
-                    imp[i].mass /= (PS::F64)(Nptcl);
-                }
-                for(PS::U32 i = 0 ; i < imp.size() ; ++ i){
-                    ptcl.push_back(imp[i]);
-                }
-                break;
-        }
+
+		switch (mode) {
+		case 1:
+			std::cout << "creating target from tar.dat" << std::endl;
+			FILE * tarFile;
+			tarFile = fopen("input/tar.dat","r");
+			FileHeader tarheader;
+			int nptcltar;
+			nptcltar = tarheader.readAscii(tarFile);
+			std::cout << "num tar ptcl: " << nptcltar << std::endl;
+			for(int i=0; i<nptcltar; i++) {
+				Ptcl ith;
+				ith.readAscii(tarFile);
+				if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
+			}
+			for(PS::U32 i = 0 ; i < tar.size() ; ++ i) {
+				tar[i].mass /= (PS::F64)(Nptcl);
+			}
+			for(PS::U32 i = 0 ; i < tar.size() ; ++ i) {
+				ptcl.push_back(tar[i]);
+			}
+
+			std::cout << "creating impactor from imp.dat" << std::endl;
+			FILE * impFile;
+			impFile = fopen("input/imp.dat","r");
+			FileHeader impheader;
+			int nptclimp;
+			nptclimp = impheader.readAscii(impFile);
+			std::cout << "num imp ptcl: " << nptclimp << std::endl;
+			for(int i=0; i<nptclimp; i++) {
+				Ptcl ith;
+				ith.readAscii(impFile);
+				ith.vel.x += (-1) * impVel;
+				if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
+			}
+			for(PS::U32 i = 0 ; i < imp.size() ; ++ i) {
+				imp[i].mass /= (PS::F64)(Nptcl);
+			}
+			for(PS::U32 i = 0 ; i < imp.size() ; ++ i) {
+				ptcl.push_back(imp[i]);
+			}
+			break;
+		case 2:
+			//Put Tar.
+			std::cout << "creating target" << std::endl;
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
+						const PS::F64 r = sqrt(x*x + y*y + z*z) * UnitRadi;
+						if(r >= tarRadi || r <= tarCoreRadi) continue;
+						Ptcl ith;
+						ith.pos.x = UnitRadi * x;
+						ith.pos.y = UnitRadi * y;
+						ith.pos.z = UnitRadi * z;
+						ith.dens = (tarMass - tarCoreMass) / (4.0 / 3.0 * math::pi * (tarRadi * tarRadi * tarRadi - tarCoreRadi * tarCoreRadi * tarCoreRadi));
+						ith.mass = tarMass + impMass;
+						ith.eng  = 0.1 * Grav * tarMass / tarRadi;
+						ith.id   = id++;
+						// TODO: Modify this line for all particles that need new EoS
+						ith.setPressure(&AGranite);
+						ith.tag = 0;
+						if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
+					}
+				}
+			}
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
+						const PS::F64 r = tarCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
+						if(r >= Corr * tarCoreRadi) continue;
+						Ptcl ith;
+						ith.pos.x = tarCoreShrinkFactor * UnitRadi * x;
+						ith.pos.y = tarCoreShrinkFactor * UnitRadi * y;
+						ith.pos.z = tarCoreShrinkFactor * UnitRadi * z;
+						ith.dens = tarCoreMass / (4.0 / 3.0 * math::pi * tarCoreRadi * tarCoreRadi * tarCoreRadi * Corr * Corr * Corr);
+						ith.mass = tarMass + impMass;
+						ith.eng  = 0.1 * Grav * tarMass / tarRadi;
+						ith.id   = id++;
+						// TODO: Modify this line for all particles that need new EoS
+						ith.setPressure(&Iron);
+						ith.tag = 1;
+						if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
+					}
+				}
+			}
+			for(PS::U32 i = 0 ; i < tar.size() ; ++ i) {
+				tar[i].mass /= (PS::F64)(Nptcl);
+			}
+			for(PS::U32 i = 0 ; i < tar.size() ; ++ i) {
+				ptcl.push_back(tar[i]);
+			}
+			break;
+		case 3:
+			//imp
+			std::cout << "creating impactor" << std::endl;
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
+						const PS::F64 r = Expand * sqrt(x*x + y*y + z*z) * UnitRadi;
+						if(r >= impRadi || r <= impCoreRadi) continue;
+						Ptcl ith;
+						ith.pos.x = Expand * UnitRadi * x + offset;
+						ith.pos.y = Expand * UnitRadi * y;
+						ith.pos.z = Expand * UnitRadi * z;
+						ith.dens = (impMass - impCoreMass) / (4.0 / 3.0 * math::pi * (impRadi * impRadi * impRadi - impCoreRadi * impCoreRadi * impCoreRadi));
+						ith.mass = tarMass + impMass;
+						ith.eng  = 0.1 * Grav * tarMass / tarRadi;
+						ith.id   = id++;
+						// TODO: Modify this line for all particles that need new EoS
+						ith.setPressure(&AGranite);
+						ith.tag = 2;
+						if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
+					}
+				}
+			}
+			for(PS::F64 x = -1.0 ; x <= 1.0 ; x += dx) {
+				for(PS::F64 y = -1.0 ; y <= 1.0 ; y += dx) {
+					for(PS::F64 z = -1.0 ; z <= 1.0 ; z += dx) {
+						const PS::F64 r = Expand * impCoreShrinkFactor * sqrt(x*x + y*y + z*z) * UnitRadi;
+						if(r >= impCoreRadi) continue;
+						Ptcl ith;
+						ith.pos.x = Expand * impCoreShrinkFactor * UnitRadi * x + offset;
+						ith.pos.y = Expand * impCoreShrinkFactor * UnitRadi * y;
+						ith.pos.z = Expand * impCoreShrinkFactor * UnitRadi * z;
+						ith.dens = impCoreMass / (4.0 / 3.0 * math::pi * impCoreRadi * impCoreRadi * impCoreRadi * Corr * Corr * Corr);
+						ith.mass = tarMass + impMass;
+						ith.eng  = 0.1 * Grav * tarMass / tarRadi;
+						ith.id   = id++;
+						// TODO: Modify this line for all particles that need new EoS
+						ith.setPressure(&Iron);
+						ith.tag = 3;
+						if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
+					}
+				}
+			}
+			for(PS::U32 i = 0 ; i < imp.size() ; ++ i) {
+				imp[i].mass /= (PS::F64)(Nptcl);
+			}
+			for(PS::U32 i = 0 ; i < imp.size() ; ++ i) {
+				ptcl.push_back(imp[i]);
+			}
+			break;
+		}
 
 		const PS::S32 numPtclLocal = ptcl.size();
 		sph_system.setNumberOfParticleLocal(numPtclLocal);
-		for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i){
+		for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i) {
 			sph_system[i] = ptcl[i];
 		}
 		//Fin.
@@ -287,22 +287,22 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		std::cout << "setup..." << std::endl;
 	}
 
-	static void setEoS(PS::ParticleSystem<Ptcl>& sph_system){
-		for(PS::U64 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+	static void setEoS(PS::ParticleSystem<Ptcl>& sph_system) {
+		for(PS::U64 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 			// TODO: Modify the lines below for all particles that need new EoS
-			if(sph_system[i].tag % 2 == 0){
+			if(sph_system[i].tag % 2 == 0) {
 				sph_system[i].setPressure(&AGranite);
-			}else{
+			} else {
 				sph_system[i].setPressure(&Iron);
 			}
 		}
 	}
 
-	static void addExternalForce(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo){
+	static void addExternalForce(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo) {
 		if(sysinfo.time >= 5000) return;
 		std::cout << "Add Ext. Force!!!" << std::endl;
 		#pragma omp parallel for
-		for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+		for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 			sph_system[i].acc += - sph_system[i].vel * 0.05 / sph_system[i].dt;
 		}
 	}
diff --git a/src/class.h b/src/class.h
index 5af1796..bff1603 100644
--- a/src/class.h
+++ b/src/class.h
@@ -1,74 +1,74 @@
 #pragma once
 
-enum TYPE{
+enum TYPE {
 	HYDRO,
 	FREEZE,
 };
 
-struct system_t{
+struct system_t {
 	PS::F64 dt, time, output_time;
 	PS::S64 step, output_id;
-	system_t() : step(0), time(0.0), output_id(0), output_time(0.0){
+	system_t() : step(0), time(0.0), output_id(0), output_time(0.0) {
 	}
 };
 
-class FileHeader{
-public:
+class FileHeader {
+  public:
 	int Nbody;
 	double time;
-	int readAscii(FILE* fp){
+	int readAscii(FILE* fp) {
 		fscanf(fp, "%lf\n", &time);
 		fscanf(fp, "%d\n", &Nbody);
 		return Nbody;
 	}
-	void writeAscii(FILE* fp) const{
+	void writeAscii(FILE* fp) const {
 		fprintf(fp, "%e\n", time);
 		fprintf(fp, "%d\n", Nbody);
 	}
 };
 
-namespace STD{
-	namespace RESULT{
+namespace STD {
+	namespace RESULT {
 		//Density summation
-		class Dens{
-			public:
+		class Dens {
+		  public:
 			PS::F64 dens;
 			PS::F64 smth;
-			void clear(){
+			void clear() {
 				dens = smth = 0;
 			}
 		};
 		//for Balsara switch
-		class Drvt{
-			public:
+		class Drvt {
+		  public:
 			PS::F64 div_v;
 			PS::F64vec rot_v;
 			PS::F64 grad_smth;
-			void clear(){
+			void clear() {
 				div_v = 0.0;
 				grad_smth = 0.0;
 				rot_v = 0.0;
 			}
 		};
 		//Hydro force
-		class Hydro{
-			public:
+		class Hydro {
+		  public:
 			PS::F64vec acc;
 			PS::F64 eng_dot;
 			PS::F64 dt;
-			void clear(){
+			void clear() {
 				acc = 0;
 				eng_dot = 0;
 				dt = 1.0e+30;
 			}
 		};
 		//Self gravity
-		class Grav{
-			public:
+		class Grav {
+		  public:
 			PS::F64vec acc;
 			PS::F64    pot;
 			PS::F64    dt;
-			void clear(){
+			void clear() {
 				acc = 0.0;
 				pot = 0.0;
 				dt = 1.0e+30;
@@ -76,8 +76,8 @@ namespace STD{
 		};
 	}
 
-	class RealPtcl{
-		public:
+	class RealPtcl {
+	  public:
 		PS::F64 mass;
 		PS::F64vec pos, vel, acc;
 		PS::F64 dens;//DENSity
@@ -103,144 +103,144 @@ namespace STD{
 
 		TYPE type;
 		//Constructor
-		RealPtcl(){
+		RealPtcl() {
 			type = HYDRO;
 			AVa = 1.0;
 			Bal = 1.0;
 		}
 		//Copy functions
-		void copyFromForce(const RESULT::Dens& dens){
+		void copyFromForce(const RESULT::Dens& dens) {
 			this->dens = dens.dens;
 			this->smth = dens.smth;
 		}
-		void copyFromForce(const RESULT::Drvt& drvt){
+		void copyFromForce(const RESULT::Drvt& drvt) {
 			this->div_v = drvt.div_v;
 			this->rot_v = drvt.rot_v;
-			if(PARAM::FLAG_B95 == true){
+			if(PARAM::FLAG_B95 == true) {
 				this->Bal = std::abs(drvt.div_v) / (std::abs(drvt.div_v) + sqrt(drvt.rot_v * drvt.rot_v) + 1.0e-4 * this->snds / this->smth); //Balsala switch
-			}else{
+			} else {
 				this->Bal = 1.0;
 			}
 			this->grad_smth = drvt.grad_smth;
 		}
-		void copyFromForce(const RESULT::Hydro& force){
+		void copyFromForce(const RESULT::Hydro& force) {
 			this->acc     = force.acc;
 			this->eng_dot = force.eng_dot;
 			this->dt      = force.dt;
 		}
-		void copyFromForce(const RESULT::Grav& force){
+		void copyFromForce(const RESULT::Grav& force) {
 			this->grav = force.acc;
 			this->pot  = force.pot;
 			//not copy dt
 		}
 		//Give necessary values to FDPS
-		PS::F64 getCharge() const{
+		PS::F64 getCharge() const {
 			return this->mass;
 		}
-		PS::F64vec getPos() const{
+		PS::F64vec getPos() const {
 			return this->pos;
 		}
-		PS::F64 getRSearch() const{
+		PS::F64 getRSearch() const {
 			return kernel_t::supportRadius() * this->smth;
 		}
-		void setPos(const PS::F64vec& pos){
+		void setPos(const PS::F64vec& pos) {
 			this->pos = pos;
 		}
-		void writeAscii(FILE* fp) const{
-			#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
-			fprintf(fp, "%ld\t%ld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",  id,  tag,  mass,  pos.x,  pos.y,  0.0  ,  vel.x,  vel.y,  0.0  ,  dens,  eng,  pres, pot);
-			#else
+		void writeAscii(FILE* fp) const {
+#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
+			fprintf(fp, "%ld\t%ld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",  id,  tag,  mass,  pos.x,  pos.y,  0.0,  vel.x,  vel.y,  0.0,  dens,  eng,  pres, pot);
+#else
 			fprintf(fp, "%ld\t%ld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",  id,  tag,  mass,  pos.x,  pos.y,  pos.z,  vel.x,  vel.y,  vel.z,  dens,  eng,  pres, pot);
-			#endif
+#endif
 		}
-		void readAscii(FILE* fp){
-			#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
+		void readAscii(FILE* fp) {
+#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 			double dummy1, dummy2;
 			fscanf (fp, "%ld\t%ld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", &id, &tag, &mass, &pos.x, &pos.y, &dummy1, &vel.x, &vel.y, &dummy2, &dens, &eng, &pres, &pot);
-			#else
+#else
 			fscanf (fp, "%ld\t%ld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", &id, &tag, &mass, &pos.x, &pos.y, &pos.z, &vel.x, &vel.y, &vel.z, &dens, &eng, &pres, &pot);
-			#endif
+#endif
 		}
-		void setPressure(const EoS::EoS_t<PS::F64>* const _EoS){
+		void setPressure(const EoS::EoS_t<PS::F64>* const _EoS) {
 			EoS = _EoS;
 		}
-		void initialize(){
+		void initialize() {
 			smth = PARAM::SMTH * pow(mass / dens, 1.0/(PS::F64)(PARAM::Dim));
 			grad_smth = 1.0;
 		}
-		void initialKick(const PS::F64 dt_glb){
+		void initialKick(const PS::F64 dt_glb) {
 			//if(type == FREEZE) return;
 			vel_half = vel + 0.5 * dt_glb * (acc + grav);
 			eng_half = eng + 0.5 * dt_glb * eng_dot;
 		}
-		void fullDrift(const PS::F64 dt_glb){
+		void fullDrift(const PS::F64 dt_glb) {
 			//if(type == FREEZE) return;
 			pos += dt_glb * vel_half;
-			if(PARAM::FLAG_R00 == true){
+			if(PARAM::FLAG_R00 == true) {
 				AVa += ((2.0 - AVa) * std::max(- div_v, 0.0) - (AVa - 0.1) / (smth / snds)) * dt_glb;
 			}
 		}
-		void predict(const PS::F64 dt_glb){
+		void predict(const PS::F64 dt_glb) {
 			//if(type == FREEZE) return;
 			vel += dt_glb * (acc + grav);
 			eng += dt_glb * eng_dot;
 		}
-		void finalKick(const PS::F64 dt_glb){
+		void finalKick(const PS::F64 dt_glb) {
 			//if(type == FREEZE) return;
 			vel = vel_half + 0.5 * dt_glb * (acc + grav);
 			eng = eng_half + 0.5 * dt_glb * eng_dot;
 		}
-        void dampMotion(double damping) {
-            vel *= damping;
-        }
+		void dampMotion(double damping) {
+			vel *= damping;
+		}
 	};
 
-	namespace EPI{
-		class Dens{
-			public:
+	namespace EPI {
+		class Dens {
+		  public:
 			PS::F64vec pos;
 			PS::F64    mass;
 			PS::F64    smth;
 			PS::S64    id;
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->pos  = rp.pos;
 				this->mass = rp.mass;
 				this->smth = rp.smth;
 				this->id   = rp.id;
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
-			void setPos(const PS::F64vec& pos){
+			void setPos(const PS::F64vec& pos) {
 				this->pos = pos;
 			}
 		};
 
-		class Drvt{
-			public:
+		class Drvt {
+		  public:
 			PS::F64vec pos;
 			PS::F64vec vel;
 			PS::F64    smth;
 			PS::F64    dens;
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				pos  = rp.pos;
 				vel  = rp.vel;
 				dens = rp.dens;
 				smth = rp.smth;
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
 		};
 
-		class Hydro{
-			public:
+		class Hydro {
+		  public:
 			PS::F64vec pos;
 			PS::F64vec vel;
 			PS::F64    smth;
@@ -251,7 +251,7 @@ namespace STD{
 			PS::F64    Bal;
 			PS::F64    AVa;
 			PS::S64    id;///DEBUG
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->pos       = rp.pos;
 				this->vel       = rp.vel;
 				this->smth      = rp.smth;
@@ -263,26 +263,26 @@ namespace STD{
 				this->AVa       = rp.AVa;
 				this->id        = rp.id;///DEBUG
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
 		};
-		class Grav{
-			public:
+		class Grav {
+		  public:
 			PS::F64vec pos;
 			PS::F64    eps2;
 			PS::S64    id;
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getEps2(void) const{
+			PS::F64 getEps2(void) const {
 				//return (1.0e-4 * 6400.0e+3) * (1.0e-4 * 6400.0e+3);//GI Unit
 				return eps2;
 			}
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				pos = rp.pos;
 				id  = rp.id;
 				//eps2 = 1.0e-2 * rp.smth * rp.smth;
@@ -291,52 +291,52 @@ namespace STD{
 		};
 	}
 
-	namespace EPJ{
-		class Dens{
-		public:
+	namespace EPJ {
+		class Dens {
+		  public:
 			PS::F64    mass;
 			PS::F64vec pos;
 			PS::F64    smth;
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->mass = rp.mass;
 				this->pos  = rp.pos;
 				this->smth = rp.smth;
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			void setPos(const PS::F64vec& pos){
+			void setPos(const PS::F64vec& pos) {
 				this->pos = pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
 		};
-		class Drvt{
-			public:
+		class Drvt {
+		  public:
 			PS::F64    mass;
 			PS::F64vec pos;
 			PS::F64vec vel;
 			PS::F64    smth;
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->mass = rp.mass;
 				this->pos  = rp.pos;
 				this->vel  = rp.vel;
 				this->smth = rp.smth;
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			void setPos(const PS::F64vec& pos){
+			void setPos(const PS::F64vec& pos) {
 				this->pos = pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
 		};
 
-		class Hydro{
-			public:
+		class Hydro {
+		  public:
 			PS::F64vec pos;
 			PS::F64vec vel;
 			PS::F64    dens;
@@ -348,7 +348,7 @@ namespace STD{
 			PS::F64    Bal;
 			PS::F64    AVa;
 			PS::S64    id;///DEBUG
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->pos   = rp.pos;
 				this->vel   = rp.vel;
 				this->dens  = rp.dens;
@@ -361,29 +361,29 @@ namespace STD{
 				this->AVa   = rp.AVa;
 				this->id    = rp.id;///DEBUG
 			}
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getRSearch() const{
+			PS::F64 getRSearch() const {
 				return kernel_t::supportRadius() * this->smth;
 			}
-			void setPos(const PS::F64vec& pos){
+			void setPos(const PS::F64vec& pos) {
 				this->pos = pos;
 			}
 		};
 
-		class Grav{
-			public:
+		class Grav {
+		  public:
 			PS::F64vec pos;
 			PS::F64    mass;
 			PS::S64    id;
-			PS::F64vec getPos() const{
+			PS::F64vec getPos() const {
 				return this->pos;
 			}
-			PS::F64 getCharge(void) const{
+			PS::F64 getCharge(void) const {
 				return this->mass;
 			}
-			void copyFromFP(const RealPtcl& rp){
+			void copyFromFP(const RealPtcl& rp) {
 				this->mass = rp.mass;
 				this->pos  = rp.pos;
 				this->id   = rp.id;
@@ -392,20 +392,20 @@ namespace STD{
 	}
 }
 
-template <class Ptcl> class Problem{
-	Problem(){
+template <class Ptcl> class Problem {
+	Problem() {
 	}
-	public:
-	static void setupIC(PS::ParticleSystem<Ptcl>&, system_t&, PS::DomainInfo&){
+  public:
+	static void setupIC(PS::ParticleSystem<Ptcl>&, system_t&, PS::DomainInfo&) {
 	}
-	static void setEoS(PS::ParticleSystem<Ptcl>&){
+	static void setEoS(PS::ParticleSystem<Ptcl>&) {
 	}
-	static void setDomain(PS::DomainInfo&){
+	static void setDomain(PS::DomainInfo&) {
 	}
-	static void addExternalForce(PS::ParticleSystem<Ptcl>&, system_t&){
+	static void addExternalForce(PS::ParticleSystem<Ptcl>&, system_t&) {
 		std::cout << "No Ext. Force" << std::endl;
 	}
-	static void postTimestepProcess(PS::ParticleSystem<Ptcl>&, system_t&){
+	static void postTimestepProcess(PS::ParticleSystem<Ptcl>&, system_t&) {
 	}
 };
 
diff --git a/src/force.h b/src/force.h
index 14b9579..1700dcb 100644
--- a/src/force.h
+++ b/src/force.h
@@ -1,39 +1,39 @@
 #pragma once
-namespace STD{
-	class CalcDensity{
+namespace STD {
+	class CalcDensity {
 		kernel_t kernel;
-		public:
-		void operator () (const EPI::Dens* const ep_i, const PS::S32 Nip, const EPJ::Dens* const ep_j, const PS::S32 Njp, RESULT::Dens* const dens){
-			for(PS::S32 i = 0 ; i < Nip ; ++ i){
+	  public:
+		void operator () (const EPI::Dens* const ep_i, const PS::S32 Nip, const EPJ::Dens* const ep_j, const PS::S32 Njp, RESULT::Dens* const dens) {
+			for(PS::S32 i = 0 ; i < Nip ; ++ i) {
 				const EPI::Dens& ith = ep_i[i];
-				for(PS::S32 j = 0 ; j < Njp ; ++ j){
+				for(PS::S32 j = 0 ; j < Njp ; ++ j) {
 					const EPJ::Dens& jth = ep_j[j];
 					const PS::F64vec dr = jth.pos - ith.pos;
 					dens[i].dens += jth.mass * kernel.W(dr, ith.smth);
 				}
-				#ifdef FLAG_GI
+#ifdef FLAG_GI
 				dens[i].dens = std::max(5.0, dens[i].dens);
-				#endif
+#endif
 				dens[i].smth = PARAM::SMTH * pow(ith.mass / dens[i].dens, 1.0/(PS::F64)(PARAM::Dim));
 			}
 		}
 	};
 
-	void CalcPressure(PS::ParticleSystem<STD::RealPtcl>& sph_system){
+	void CalcPressure(PS::ParticleSystem<STD::RealPtcl>& sph_system) {
 		#pragma omp parallel for
-		for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+		for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 			sph_system[i].pres = sph_system[i].EoS->Pressure(sph_system[i].dens, sph_system[i].eng);
 			sph_system[i].snds = sph_system[i].EoS->SoundSpeed(sph_system[i].dens, sph_system[i].eng);
 		}
 	}
 
-	class CalcDerivative{
+	class CalcDerivative {
 		kernel_t kernel;
-		public:
-		void operator () (const EPI::Drvt* ep_i, const PS::S32 Nip, const EPJ::Drvt* ep_j, const PS::S32 Njp, RESULT::Drvt* const drvt){
-			for(PS::S32 i = 0; i < Nip ; ++ i){
+	  public:
+		void operator () (const EPI::Drvt* ep_i, const PS::S32 Nip, const EPJ::Drvt* ep_j, const PS::S32 Njp, RESULT::Drvt* const drvt) {
+			for(PS::S32 i = 0; i < Nip ; ++ i) {
 				const EPI::Drvt& ith = ep_i[i];
-				for(PS::S32 j = 0; j < Njp ; ++ j){
+				for(PS::S32 j = 0; j < Njp ; ++ j) {
 					const EPJ::Drvt& jth = ep_j[j];
 					const PS::F64vec dr = ith.pos - jth.pos;
 					const PS::F64vec dv = ith.vel - jth.vel;
@@ -48,14 +48,14 @@ namespace STD{
 		}
 	};
 
-	class CalcHydroForce{
+	class CalcHydroForce {
 		const kernel_t kernel;
-		public:
-		void operator () (const EPI::Hydro* const ep_i, const PS::S32 Nip, const EPJ::Hydro* const ep_j, const PS::S32 Njp, RESULT::Hydro* const hydro){
-			for(PS::S32 i = 0; i < Nip ; ++ i){
+	  public:
+		void operator () (const EPI::Hydro* const ep_i, const PS::S32 Nip, const EPJ::Hydro* const ep_j, const PS::S32 Njp, RESULT::Hydro* const hydro) {
+			for(PS::S32 i = 0; i < Nip ; ++ i) {
 				PS::F64 v_sig_max = 0.0;
 				const EPI::Hydro& ith = ep_i[i];
-				for(PS::S32 j = 0; j < Njp ; ++ j){
+				for(PS::S32 j = 0; j < Njp ; ++ j) {
 					const EPJ::Hydro& jth = ep_j[j];
 					const PS::F64vec dr = ith.pos - jth.pos;
 					const PS::F64vec dv = ith.vel - jth.vel;
@@ -63,31 +63,31 @@ namespace STD{
 					const PS::F64 v_sig = ith.snds + jth.snds - 3.0 * w_ij;
 					v_sig_max = std::max(v_sig_max, v_sig);
 					PS::F64 AV = - PARAM::AV_STRENGTH * 0.5 * v_sig * w_ij / (0.5 * (ith.dens + jth.dens)) * 0.5 * (ith.Bal + jth.Bal);
-					if(PARAM::FLAG_R00 == true){
+					if(PARAM::FLAG_R00 == true) {
 						AV *= 0.5 * (ith.AVa + jth.AVa);
 					}
-					#if 1
+#if 1
 					const PS::F64vec gradW = 0.5 * (kernel.gradW(dr, ith.smth) * ith.grad_smth + kernel.gradW(dr, jth.smth) * jth.grad_smth);
 					hydro[i].acc     -= jth.mass * (ith.grad_smth * ith.pres / (ith.dens * ith.dens) * kernel.gradW(dr, ith.smth) + jth.grad_smth * jth.pres / (jth.dens * jth.dens) * kernel.gradW(dr, jth.smth) + AV * gradW);
 					hydro[i].eng_dot += jth.mass * (ith.grad_smth * ith.pres / (ith.dens * ith.dens) + 0.5 * AV) * dv * gradW;
-					#else
+#else
 					const PS::F64vec gradW = 0.5 * (kernel.gradW(dr, ith.smth) + kernel.gradW(dr, jth.smth));
 					hydro[i].acc     -= jth.mass * (ith.pres / (ith.dens * ith.dens) + jth.pres / (jth.dens * jth.dens) + AV) * gradW;
 					hydro[i].eng_dot += jth.mass * (ith.pres / (ith.dens * ith.dens) + 0.5 * AV) * dv * gradW;
-					#endif
+#endif
 				}
 				hydro[i].dt = PARAM::C_CFL * 2.0 * ith.smth / v_sig_max;
 			}
 		}
 	};
 
-	template <class TPtclJ> class CalcGravityForce{
+	template <class TPtclJ> class CalcGravityForce {
 		static const double G;
-		public:
-		void operator () (const EPI::Grav* const __restrict ep_i, const PS::S32 Nip, const TPtclJ* const __restrict ep_j, const PS::S32 Njp, RESULT::Grav* const grav){
-			for(PS::S32 i = 0; i < Nip ; ++ i){
+	  public:
+		void operator () (const EPI::Grav* const __restrict ep_i, const PS::S32 Nip, const TPtclJ* const __restrict ep_j, const PS::S32 Njp, RESULT::Grav* const grav) {
+			for(PS::S32 i = 0; i < Nip ; ++ i) {
 				const EPI::Grav& ith = ep_i[i];
-				for(PS::S32 j = 0; j < Njp ; ++ j){
+				for(PS::S32 j = 0; j < Njp ; ++ j) {
 					const TPtclJ& jth = ep_j[j];
 					const PS::F64vec dr = ith.pos - jth.pos;
 					const PS::F64 dr2 = dr * dr;
diff --git a/src/init/test.h b/src/init/test.h
index 21edc9a..f50f14a 100644
--- a/src/init/test.h
+++ b/src/init/test.h
@@ -49,10 +49,10 @@ void SetupIC(PS::ParticleSystem<RealPtcl>& sph_system, PS::F64* end_time, PS::Do
 }
 */
 
-template <class Ptcl> class GI : public Problem<Ptcl>{
-	public:
+template <class Ptcl> class GI : public Problem<Ptcl> {
+  public:
 	static const double END_TIME;
-	static void setupIC(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo, PS::DomainInfo& dinfo){
+	static void setupIC(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo, PS::DomainInfo& dinfo) {
 
 		/////////
 		//place ptcls
@@ -64,9 +64,9 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		const PS::F64 Grav = 1.0;
 		const PS::F64 dx = 1.0 / 8.0;
 		PS::S32 id = 0;
-		for(PS::F64 x = -Radi ; x <= Radi ; x += dx){
-			for(PS::F64 y = -Radi ; y <= Radi ; y += dx){
-				for(PS::F64 z = -Radi ; z <= Radi ; z += dx){
+		for(PS::F64 x = -Radi ; x <= Radi ; x += dx) {
+			for(PS::F64 y = -Radi ; y <= Radi ; y += dx) {
+				for(PS::F64 z = -Radi ; z <= Radi ; z += dx) {
 					const PS::F64 r = sqrt(x*x + y*y + z*z);
 					if(r > Radi) continue;
 					Ptcl ith;
@@ -83,17 +83,17 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 				}
 			}
 		}
-		for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i){
+		for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i) {
 			ptcl[i].mass = ptcl[i].mass / (PS::F64)(ptcl.size());
 		}
 
-		if(PS::Comm::getRank() == 0){
+		if(PS::Comm::getRank() == 0) {
 			const PS::S32 numPtclLocal = ptcl.size();
 			sph_system.setNumberOfParticleLocal(numPtclLocal);
-			for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i){
+			for(PS::U32 i = 0 ; i < ptcl.size() ; ++ i) {
 				sph_system[i] = ptcl[i];
 			}
-		}else{
+		} else {
 			sph_system.setNumberOfParticleLocal(0);
 		}
 		//Fin.
@@ -101,13 +101,13 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 
 	}
 
-	static void setEoS(PS::ParticleSystem<Ptcl>& sph_system){
-		for(PS::U64 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+	static void setEoS(PS::ParticleSystem<Ptcl>& sph_system) {
+		for(PS::U64 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 			sph_system[i].setPressure(&Monoatomic);
 		}
 	}
 
-	static void addExternalForce(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo){
+	static void addExternalForce(PS::ParticleSystem<Ptcl>& sph_system, system_t& sysinfo) {
 	}
 };
 
diff --git a/src/integral.h b/src/integral.h
index 41dbc1f..da91e97 100644
--- a/src/integral.h
+++ b/src/integral.h
@@ -1,18 +1,18 @@
 #pragma once
-template <class ThisPtcl> PS::F64 getTimeStepGlobal(const PS::ParticleSystem<ThisPtcl>& sph_system){
+template <class ThisPtcl> PS::F64 getTimeStepGlobal(const PS::ParticleSystem<ThisPtcl>& sph_system) {
 	PS::F64 dt = 1.0e+30;//set VERY LARGE VALUE
-	for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+	for(PS::S32 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 		PS::F64 dt_tmp = 1.0e+30;
-		#if 1
+#if 1
 		const PS::F64 acc = sqrt((sph_system[i].grav + sph_system[i].acc) * (sph_system[i].grav + sph_system[i].acc));
 		//dt_tmp = std::min(dt_tmp, sqrt(sph_system[i].smth / (acc + 1.0e-16)));
 		//dt_tmp = std::min(dt_tmp, std::abs(sph_system[i].eng) / std::abs(sph_system[i].eng_dot + 1.0e-16));
 		//dt_tmp = std::min(dt_tmp, 1.0 / std::abs(sph_system[i].div_v + 1.0e-16));
 		dt = std::min(dt, std::min(sph_system[i].dt, dt_tmp));
-		#else
+#else
 		dt_tmp = std::max(sph_system[i].dt, PARAM::C_CFL * sph_system[i].smth / (sph_system[i].snds + 1.2 * (sph_system[i].snds + 2.0 * sph_system[i].smth * std::min(sph_system[i].div_v, 0.0))));
 		dt = std::min(dt, dt_tmp);
-		#endif
+#endif
 	}
 	return PS::Comm::getMinValue(dt);
 }
diff --git a/src/io.h b/src/io.h
index 3536374..9dd301e 100644
--- a/src/io.h
+++ b/src/io.h
@@ -4,37 +4,37 @@
 #include <dirent.h>
 #include <cassert>
 
-template <class ThisPtcl> void OutputBinary(PS::ParticleSystem<ThisPtcl>& sph_system, const system_t& sysinfo, std::string &out_dir){
+template <class ThisPtcl> void OutputBinary(PS::ParticleSystem<ThisPtcl>& sph_system, const system_t& sysinfo, std::string &out_dir) {
 	//Binary
 	char filename[256];
 	std::ofstream fout;
-    sprintf(filename, "%s/%05d_%05d.bin", out_dir.c_str(), PS::Comm::getNumberOfProc(), PS::Comm::getRank());
+	sprintf(filename, "%s/%05d_%05d.bin", out_dir.c_str(), PS::Comm::getNumberOfProc(), PS::Comm::getRank());
 	fout.open(filename, std::ios::out | std::ios::binary | std::ios::trunc);
-	if(!fout){
+	if(!fout) {
 		std::cout << "cannot write restart file." << std::endl;
 		exit(1);
 	}
 	fout.write(reinterpret_cast<const char * const>(&sysinfo), sizeof(system_t));
-	for(std::size_t i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
+	for(std::size_t i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i) {
 		ThisPtcl ith = sph_system[i];
 		fout.write((char*)&ith, sizeof(ThisPtcl));
 	}
 	fout.close();
 }
 
-template <class ThisPtcl> void OutputFileWithTimeInterval(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo, const PS::F64 end_time, std::string &out_dir){
-	if(sysinfo.time > sysinfo.output_time){
+template <class ThisPtcl> void OutputFileWithTimeInterval(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo, const PS::F64 end_time, std::string &out_dir) {
+	if(sysinfo.time > sysinfo.output_time) {
 		FileHeader header;
 		header.time = sysinfo.time;
 		header.Nbody = sph_system.getNumberOfParticleLocal();
 		char filename[20];
-        sprintf(filename, "results.%05d", sysinfo.output_id);
+		sprintf(filename, "results.%05d", sysinfo.output_id);
 		std::string full_filename = out_dir + filename;
 
 		sph_system.writeParticleAscii(full_filename.c_str(), "%s_%05d_%05d.dat", header);
-		if(PS::Comm::getRank() == 0){
-            std::cout << "//================================" << std::endl;
-            std::cout << "output " << full_filename << "." << std::endl;
+		if(PS::Comm::getRank() == 0) {
+			std::cout << "//================================" << std::endl;
+			std::cout << "output " << full_filename << "." << std::endl;
 			std::cout << "//================================" << std::endl;
 		}
 		sysinfo.output_time += end_time / PARAM::NUMBER_OF_SNAPSHOTS;
@@ -42,28 +42,28 @@ template <class ThisPtcl> void OutputFileWithTimeInterval(PS::ParticleSystem<Thi
 	}
 }
 
-template <class ThisPtcl> void InputFileWithTimeInterval(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo){
+template <class ThisPtcl> void InputFileWithTimeInterval(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo) {
 	FileHeader header;
 	char filename[256];
-    sprintf(filename, "results/%05d", sysinfo.step);
+	sprintf(filename, "results/%05d", sysinfo.step);
 	sph_system.readParticleAscii(filename, "%s_%05d_%05d.dat", header);
 	sysinfo.time = header.time;
 	std::cout << header.time << std::endl;
 }
 
-template <class ThisPtcl> void InputBinary(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo){
+template <class ThisPtcl> void InputBinary(PS::ParticleSystem<ThisPtcl>& sph_system, system_t& sysinfo) {
 	char filename[256];
 	//sprintf(filename, "result/%05d_%05d_%05d.bin", sysinfo->step, PS::Comm::getNumberOfProc(), PS::Comm::getRank());
 	sprintf(filename, "./%05d_%05d.bin", PS::Comm::getNumberOfProc(), PS::Comm::getRank());
 	std::ifstream fin(filename, std::ios::in | std::ios::binary);
-	if(!fin){
+	if(!fin) {
 		std::cout << "cannot open restart file." << std::endl;
 		exit(1);
 	}
 	std::vector<ThisPtcl> ptcl;
 	fin.read((char*)&sysinfo, sizeof(system_t));
 	//sysinfo.step ++;
-	while(1){
+	while(1) {
 		ThisPtcl ith;
 		fin.read((char*)&ith, sizeof(ThisPtcl));
 		if(fin.eof() == true) break;
@@ -71,19 +71,18 @@ template <class ThisPtcl> void InputBinary(PS::ParticleSystem<ThisPtcl>& sph_sys
 	}
 	fin.close();
 	sph_system.setNumberOfParticleLocal(ptcl.size());
-		for(std::size_t i = 0 ; i < ptcl.size() ; ++ i){
+	for(std::size_t i = 0 ; i < ptcl.size() ; ++ i) {
 		sph_system[i] = ptcl[i];
 	}
 	std::cout << "read binary" << std::endl;
 }
 
-void createOutputDirectory(const std::string &directory_name){
+void createOutputDirectory(const std::string &directory_name) {
 	// check if the output directory exists, if not create it.
-	if (DIR *output_directory = opendir(directory_name.c_str())){
+	if (DIR *output_directory = opendir(directory_name.c_str())) {
 		int error = closedir(output_directory);
 		assert(error == 0);
-	}
-	else{
+	} else {
 		int error = mkdir(directory_name.c_str(),S_IRWXU | S_IRGRP | S_IXGRP | S_IROTH | S_IXOTH);
 		assert(error == 0);
 	}
diff --git a/src/kernel.h b/src/kernel.h
index cb6b4ae..7fbf906 100644
--- a/src/kernel.h
+++ b/src/kernel.h
@@ -1,35 +1,35 @@
 #pragma once
 
 //Wendland C6
-struct WendlandC6{
-	WendlandC6(){}
+struct WendlandC6 {
+	WendlandC6() {}
 	//W
-	PS::F64 W(const PS::F64vec dr, const PS::F64 h) const{
+	PS::F64 W(const PS::F64vec dr, const PS::F64 h) const {
 		const PS::F64 H = supportRadius() * h;
 		const PS::F64 s = sqrt(dr * dr) / H;
 		PS::F64 r_value;
 		r_value = (1.0 + s * (8.0 + s * (25.0 + s * (32.0)))) * math::pow8(math::plus(1.0 - s));
-		#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
+#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 		r_value *= (78./7.) / (H * H * math::pi);
-		#else
+#else
 		r_value *= (1365./64.) / (H * H * H * math::pi);
-		#endif
+#endif
 		return r_value;
 	}
 	//gradW
-	PS::F64vec gradW(const PS::F64vec dr, const PS::F64 h) const{
+	PS::F64vec gradW(const PS::F64vec dr, const PS::F64 h) const {
 		const PS::F64 H = supportRadius() * h;
 		const PS::F64 s = sqrt(dr * dr) / H;
 		PS::F64 r_value;
 		r_value = math::pow7(math::plus(1.0 - s)) * (math::plus(1.0 - s) * (8.0 + s * (50.0 + s * (96.0))) - 8.0 * (1.0 + s * (8.0 + s * (25.0 + s * (32.0)))));
-		#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
+#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 		r_value *= (78./7.) / (H * H * math::pi);
-		#else
+#else
 		r_value *= (1365./64.) / (H * H * H * math::pi);
-		#endif
+#endif
 		return dr * r_value / (sqrt(dr * dr) * H  + 1.0e-6 * h);
 	}
-	static PS::F64 supportRadius(){
+	static PS::F64 supportRadius() {
 		return 2.5;
 	}
 };
diff --git a/src/mathfunc.h b/src/mathfunc.h
index f7efe9f..fc909e7 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -1,38 +1,38 @@
 #pragma once
 
-namespace math{
+namespace math {
 	const PS::F64 pi = atan(1.0) * 4.0;
 	const PS::F64 NaN = + 0.0 / 0.0;
 	const PS::F64 VERY_LARGE_VALUE = 1.0e+30;
-	template <typename type> inline type plus(const type arg){
+	template <typename type> inline type plus(const type arg) {
 		return (arg > 0) ? arg : 0;
 	}
-	template <typename type> inline type sign(const type arg){
+	template <typename type> inline type sign(const type arg) {
 		return (arg > 0) ? 1.0 : - 1.0;
 	}
-	template <typename type> inline type pow2(const type arg){
+	template <typename type> inline type pow2(const type arg) {
 		return arg * arg;
 	}
-	template <typename type> inline type pow3(const type arg){
+	template <typename type> inline type pow3(const type arg) {
 		return arg * arg * arg;
 	}
-	template <typename type> inline type pow4(const type arg){
+	template <typename type> inline type pow4(const type arg) {
 		const type arg2 = arg * arg;
 		return arg2 * arg2;
 	}
-	template <typename type> inline type pow5(const type arg){
+	template <typename type> inline type pow5(const type arg) {
 		const type arg2 = arg * arg;
 		return arg2 * arg2 * arg;
 	}
-	template <typename type> inline type pow6(const type arg){
+	template <typename type> inline type pow6(const type arg) {
 		const type arg3 = arg * arg * arg;
 		return arg3 * arg3;
 	}
-	template <typename type> inline type pow7(const type arg){
+	template <typename type> inline type pow7(const type arg) {
 		const type arg3 = arg * arg * arg;
 		return arg3 * arg3 * arg;
 	}
-	template <typename type> inline type pow8(const type arg){
+	template <typename type> inline type pow8(const type arg) {
 		const type arg2 = arg * arg;
 		const type arg4 = arg2 * arg2;
 		return arg4 * arg4;
diff --git a/src/param.h b/src/param.h
index ddd4b77..92e7ff9 100644
--- a/src/param.h
+++ b/src/param.h
@@ -1,11 +1,11 @@
 #pragma once
 
-namespace PARAM{
-	#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
+namespace PARAM {
+#ifdef PARTICLE_SIMULATOR_TWO_DIMENSION
 	const short int Dim = 2;
-	#else
+#else
 	const short int Dim = 3;
-	#endif
+#endif
 	const PS::F64 SMTH = 1.2;
 	const PS::F64 C_CFL = 0.3;
 	const PS::U64 NUMBER_OF_SNAPSHOTS = 200;
diff --git a/src/parse.h b/src/parse.h
index 1db02ca..1c65e5e 100644
--- a/src/parse.h
+++ b/src/parse.h
@@ -12,48 +12,43 @@
  * whitespace at the beginning and end of a string.
  */
 std::string
-trim(const std::string &input)
-{
-std::string::size_type left = 0;
-std::string::size_type right = input.size() > 0 ? input.size() - 1 : 0;
+trim(const std::string &input) {
+	std::string::size_type left = 0;
+	std::string::size_type right = input.size() > 0 ? input.size() - 1 : 0;
 
-for (; left < input.size(); ++left)
-  {
-	if (!std::isspace(input[left]))
-	  {
-		break;
-	  }
-  }
+	for (; left < input.size(); ++left) {
+		if (!std::isspace(input[left])) {
+			break;
+		}
+	}
 
-for (; right >= left; --right)
-  {
-	if (!std::isspace(input[right]))
-	  {
-		break;
-	  }
-  }
+	for (; right >= left; --right) {
+		if (!std::isspace(input[right])) {
+			break;
+		}
+	}
 
-return std::string(input, left, right - left + 1);
+	return std::string(input, left, right - left + 1);
 }
 
 /**
  * This class handles reading in model input parameters from input files.
  */
-class ParameterFile{
+class ParameterFile {
 	std::map<std::string, std::string> param;
-	public:
-	ParameterFile(const std::string &file){
+  public:
+	ParameterFile(const std::string &file) {
 		std::ifstream input;
 		input.open(file, std::ios::in);
-		if(input.fail() == true){
+		if(input.fail() == true) {
 			std::cout << "Cannot open file, assuming default values ..." << std::endl;
 			return;
 		}
 		std::string line;
-		while(std::getline(input, line)){
+		while(std::getline(input, line)) {
 			// Remove comments
 			auto trailing_comment = std::find(line.begin(),line.end(),'#');
-			if (trailing_comment != line.end()){
+			if (trailing_comment != line.end()) {
 				line.erase(trailing_comment,line.end());
 			}
 			line = trim(line);
@@ -63,39 +58,37 @@ class ParameterFile{
 			std::istringstream stream(line);
 			std::string tmp;
 			std::vector<std::string> keyval;
-			while(std::getline(stream, tmp, '=')){
+			while(std::getline(stream, tmp, '=')) {
 				keyval.push_back(trim(tmp));
 			}
 			param.insert(std::map<std::string, std::string>::value_type(keyval[0], keyval[1]));
 		}
 		input.close();
 	}
-	void show(){
-		for(std::map<std::string, std::string>::iterator it = param.begin() ; it != param.end() ; ++ it){
+	void show() {
+		for(std::map<std::string, std::string>::iterator it = param.begin() ; it != param.end() ; ++ it) {
 			std::cout << it->first << " = " << it->second << std::endl;
 		}
 	}
 	template <typename T>
-	T getValueOf(const std::string &key, const T &default_parameter_value){
+	T getValueOf(const std::string &key, const T &default_parameter_value) {
 
 		T parameter_value = default_parameter_value;
 
-		if (param.find(key) != param.end()){
+		if (param.find(key) != param.end()) {
 			parameter_value = T();
 			std::stringstream convert_from_string_to_value(param[key]);
 			convert_from_string_to_value >> parameter_value;
 
-		    if(convert_from_string_to_value.fail())
-		    {
-		      std::cout << "Convert error: cannot convert string '" << param[key] << "' to value" << std::endl;
-		      throw;
-		    }
-		}
-		else{
+			if(convert_from_string_to_value.fail()) {
+				std::cout << "Convert error: cannot convert string '" << param[key] << "' to value" << std::endl;
+				throw;
+			}
+		} else {
 			std::cout << "Parameter <" << key << "> not found and no default value provided. Aborting." << std::endl;
 			throw;
 		}
 
 		return parameter_value;
-	  }
+	}
 };