diff --git a/demos/scenario_24_bidirection/bidirection_geo_SimplestModel.xml b/demos/scenario_24_bidirection/bidirection_geo_SimplestModel.xml
new file mode 100644
index 0000000000000000000000000000000000000000..56c62550002ae94d5dbd8cce0f69a1c79b4f73b7
--- /dev/null
+++ b/demos/scenario_24_bidirection/bidirection_geo_SimplestModel.xml
@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<geometry version="0.8" caption="FD1d" unit="m">
+  <rooms>
+	<room id="0" caption="corridor">
+      <subroom id="0" closed="0" class="subroom">
+        <polygon>	
+          <vertex px="26" py="4.0" />
+          <vertex px="0" py="4.0" />
+        </polygon>
+        <polygon>	
+          <vertex px="26" py="-4.0" />
+          <vertex px="0" py="-4.0" />
+        </polygon>	  
+      </subroom>
+    </room>
+  </rooms>
+
+  <transitions>
+    <transition id="0" caption="left" type="emergency"
+                room1_id="0" subroom1_id="0" room2_id="-1" subroom2_id="-1">
+      <vertex px="0" py="-4.0" />
+      <vertex px="0" py="4.0" />
+    </transition>
+
+    <transition id="1" caption="right" type="emergency"
+                room1_id="0" subroom1_id="0" room2_id="-1" subroom2_id="-1">
+      <vertex px="26" py="-4.0" />
+      <vertex px="26" py="4.0" />
+    </transition>
+  </transitions>
+</geometry>
diff --git a/demos/scenario_24_bidirection/bidirection_ini_SimplestModel.xml b/demos/scenario_24_bidirection/bidirection_ini_SimplestModel.xml
new file mode 100644
index 0000000000000000000000000000000000000000..daafc2980c22a9d7e587360563c19826bd48b5a0
--- /dev/null
+++ b/demos/scenario_24_bidirection/bidirection_ini_SimplestModel.xml
@@ -0,0 +1,87 @@
+<JuPedSim xmlns:ns0="http://xsd.jupedsim.org/jps_ini_core.xsd" project="JPS-Project" version="0.8" ns0:noNamespaceSchemaLocation="ini.xsd">
+  <seed>1.0</seed>
+  
+  <num_threads>8</num_threads>
+  <max_sim_time unit="sec">200</max_sim_time>
+  <logfile>log.txt</logfile>
+  <geometry>bidirection_geo_SimplestModel.xml</geometry>
+  
+  <trajectories embed_mesh="false" format="xml-plain" fps="20" color_mode ="group">
+    <file location="bidirection_traj_SimplestModel.xml" /> 
+  </trajectories>
+  
+  <show_statistics>false</show_statistics>
+  
+  <routing>
+	<goals>
+		<goal id="0" final="true" caption="goal 0">
+			<polygon>
+				<vertex px="29" py="-2" />
+				<vertex px="29" py="2" />
+				<vertex px="30" py="2" />
+				<vertex px="30" py="-2" />
+				<vertex px="29" py="-2" />
+			</polygon>
+		</goal>
+		<goal id="1" final="true" caption="goal 1">
+			<polygon>
+				<vertex px="-4" py="-2" />
+				<vertex px="-4" py="2" />
+				<vertex px="-3" py="2" />
+				<vertex px="-3" py="-2" />
+				<vertex px="-4" py="-2" />
+			</polygon>
+		</goal>
+	</goals>
+  </routing>
+  
+  
+  
+  
+  
+  
+  <agents operational_model_id="7">
+    <agents_distribution>
+		<group agent_parameter_id="1" goal_id="0" group_id="0" number="100" room_id="0" router_id="1" subroom_id="0" />
+		<group agent_parameter_id="1" goal_id="1" group_id="1" number="100" room_id="0" router_id="1" subroom_id="0" />			
+	</agents_distribution>
+  </agents>
+ 
+  <operational_models>
+	<model description="simplest" operational_model_id="7">
+      <model_parameters>
+        <solver>euler</solver>
+        <stepsize>0.05</stepsize> <!-- Set time steps -->
+        <exit_crossing_strategy>3</exit_crossing_strategy>
+	<periodic>1</periodic>
+        <linkedcells cell_size="300" enabled="true" />
+        <force_ped D="0.1" a="3" />
+        <force_wall D="0.05" a="6" />
+		<time_parameters Td="0.3" Ts="0.3" />
+		<update_method parallel="1"/> <!-- Update method is parallel when parallel=1, while unparallel when parallel=0 -->
+		<waiting_time Tw="0"/> <!-- Deleting pedestrian Tw seconds after the clogging -->
+		<clogging_area size="100" /> <!-- Only deleting pedestrian when the distance between clogging to exit is smaller than the size-->
+		<model_submodel direction="1" speed="1"/> <!-- When direction=1 means using direction part, when speed=1 means using speed part -->
+		<GCVM using="1"/> <!-- When using=1 menas we using GCVM in the simulation -->
+      </model_parameters>
+	  <agent_parameters agent_parameter_id="1">
+        <v0 mu="1.34" sigma="0.26" />
+        <bmax mu="0.20" sigma="0.00000" /> 
+        <bmin mu="0.20" sigma="0.00000" />
+        <amin mu="0.20" sigma="0.00000" />
+        <atau mu="0.00" sigma="0.00000" />
+		<tau mu="0.5" sigma="0.000" />
+		<shape circle="1" /> 
+      </agent_parameters>
+    </model>
+  </operational_models>
+
+  <route_choice_models>
+    <router description="local_shortest" router_id="1">
+      <parameters>
+        
+      </parameters>
+    </router>
+  </route_choice_models>
+
+</JuPedSim>
\ No newline at end of file
diff --git a/math/SimplestModel.cpp b/math/SimplestModel.cpp
index 65fdba942ef17aa1f8628fb00395f0cf7ac11d6a..6d708107b81bc9e733ec70d666f6079a3447a231 100644
--- a/math/SimplestModel.cpp
+++ b/math/SimplestModel.cpp
@@ -368,12 +368,21 @@ void SimplestModel::ComputeNextTimeStep(double current, double deltaT, Building*
 			e.SetSinPhi(sinPhi);
 			ped->SetEllipse(e);
 			ped->SetV(speed);
-			ped->SetPos(pos_neu);
+			
 			if (periodic) {
-				if (ped->GetPos()._x >= xRight_simplest) {
-					ped->SetPos(Point(ped->GetPos()._x - (xRight_simplest - xLeft_simplest), ped->GetPos()._y));
+				if ((ped->GetPos()._x < xRight_simplest)&&(pos_neu._x>=xRight_simplest)) {
+					ped->SetPos(Point(pos_neu._x - (xRight_simplest - xLeft_simplest), pos_neu._y));
+				}
+				else if ((ped->GetPos()._x > xLeft_simplest) && (pos_neu._x <= xLeft_simplest)) {
+					ped->SetPos(Point(pos_neu._x + (xRight_simplest - xLeft_simplest), pos_neu._y));
+				}
+				else {
+					ped->SetPos(pos_neu);
 				}
 			}
+			else {
+				ped->SetPos(pos_neu);
+			}
 		}
 
 	}
@@ -400,12 +409,21 @@ void SimplestModel::ComputeNextTimeStep(double current, double deltaT, Building*
 			e.SetSinPhi(sinPhi);
 			ped->SetEllipse(e);
 			ped->SetV(v_neu);
-			ped->SetPos(pos_neu);
+
 			if (periodic) {
-				if (ped->GetPos()._x >= xRight_simplest) {
-					ped->SetPos(Point(ped->GetPos()._x - (xRight_simplest - xLeft_simplest), ped->GetPos()._y));
+				if ((ped->GetPos()._x < xRight_simplest)&&(pos_neu._x>=xRight_simplest)) {
+					ped->SetPos(Point(pos_neu._x - (xRight_simplest - xLeft_simplest), pos_neu._y));
+				}
+				else if ((ped->GetPos()._x > xLeft_simplest) && (pos_neu._x <= xLeft_simplest)) {
+					ped->SetPos(Point(pos_neu._x + (xRight_simplest - xLeft_simplest), pos_neu._y));
+				}
+				else {
+					ped->SetPos(pos_neu);
 				}
 			}
+			else {
+				ped->SetPos(pos_neu);
+			}
 		}
 	}
 	
@@ -429,7 +447,8 @@ void SimplestModel::ComputeNextTimeStep(double current, double deltaT, Building*
 		}
 		if (converse != relations.end())
 		{
-			*converse = ID_pair(first_ID, second_ID);
+			//using this when delete pedestrian
+			//*converse = ID_pair(first_ID, second_ID);
 		}
 		for (int p = start; p <= end; ++p) {
 			Pedestrian* ped = allPeds[p];
@@ -441,7 +460,39 @@ void SimplestModel::ComputeNextTimeStep(double current, double deltaT, Building*
 				}
 				else {
 					ped->SetInCloggingTime(0);
-					pedsToRemove.push_back(ped);
+					double velocity_x=ped->GetEllipse().GetCosPhi();
+					double velocity_y=ped->GetEllipse().GetSinPhi();
+					Point position=ped->GetPos();
+					
+					int random = rand() % 10000;
+					if (random<2500)
+					{
+						// moving forward
+						Point velocity(velocity_x,velocity_y);
+						ped->SetPos(position+velocity*1.34*deltaT);
+					}
+					else if (random < 5000)
+					{
+						// Moving backward
+						Point velocity(velocity_x,velocity_y);
+						ped->SetPos(position+velocity*-1.34*deltaT);
+					}
+					else if (random < 7500)
+					{
+						// Moving left
+						Point velocity(velocity_y,velocity_x);
+						ped->SetPos(position+velocity*1.34*deltaT);
+					}
+					else
+					{
+						// moving right
+						Point velocity(velocity_y,velocity_x);
+						ped->SetPos(position+velocity*-1.34*deltaT);
+					}
+					
+					// Clogging experiment
+					//pedsToRemove.push_back(ped);
+
 					clogging_times++;
 					std::ofstream ofile;
 					string ProjectFileName = building->GetProjectFilename();
@@ -548,12 +599,17 @@ my_pair SimplestModel::GetSpacing(Pedestrian* ped1, Pedestrian* ped2, Point ei,
 	double x1 = ped1->GetPos()._x;
 	double x2_real = ped2->GetPos()._x;
 	double y2 = ped2->GetPos()._y;
+	Point ped2_current = ped2->GetPos();
 
 	if (periodic) {
 		if ((xRight_simplest - x1) + (x2_real - xLeft_simplest) <= cutoff_simplest) {
 			double x2_periodic = x2_real + xRight_simplest - xLeft_simplest;
 			ped2->SetPos(Point(x2_periodic, y2));
 		}
+		if ((x1 - xLeft_simplest) + (xRight_simplest - x2_real) <= cutoff_simplest) {
+			double x2_periodic = xLeft_simplest - xRight_simplest + x2_real;
+			ped2->SetPos(Point(x2_periodic, y2));
+		}
 	}
 
 	Point distp12 = ped2->GetPos() - ped1->GetPos(); // inversed sign
@@ -574,6 +630,7 @@ my_pair SimplestModel::GetSpacing(Pedestrian* ped1, Pedestrian* ped2, Point ei,
 	double eff_dist = eped1.EffectiveDistanceToEllipse(eped2, &dist);
 	double condition1 = ei.ScalarProduct(ep12); // < e_i , e_ij > should be positive
 	if (condition1 <= 0) {
+		ped2->SetPos(ped2_current);
 		return  my_pair(FLT_MAX, ped2->GetID());
 	}
 	if (!ped1->GetEllipse().DoesStretch())
@@ -581,6 +638,7 @@ my_pair SimplestModel::GetSpacing(Pedestrian* ped1, Pedestrian* ped2, Point ei,
 		double l = ped1->GetLargerAxis() + ped2->GetLargerAxis();
 		double condition2 = ei.Rotate(0, 1).ScalarProduct(ep12); // theta = pi/2. condition2 should <= than l/Distance
 		condition2 = (condition2 > 0) ? condition2 : -condition2; // abs
+		ped2->SetPos(ped2_current);
 		if ((condition1 >= 0) && (condition2 <= l / Distance))
 			// return a pair <dist, condition1>. Then take the smallest dist. In case of equality the biggest condition1
 			return  my_pair(distp12.Norm() - l, ped2->GetID());
@@ -621,6 +679,7 @@ my_pair SimplestModel::GetSpacing(Pedestrian* ped1, Pedestrian* ped2, Point ei,
 		double d2 = -sinphi1 * (x2 - x1) + cosphi1 * (y2 - y1) - b1;
 		if (d1*d2 <= 0) {
 			//if the center between two lines, collision
+			ped2->SetPos(ped2_current);
 			return  my_pair(eff_dist, ped2->GetID());
 		}
 		//If the center not between two lines, Judge if ped2 contact with two lines
@@ -704,6 +763,10 @@ Point SimplestModel::ForceRepPed(Pedestrian* ped1, Pedestrian* ped2, Point e0, i
 			double x2_periodic = x_j + xRight_simplest - xLeft_simplest;
 			ped2->SetPos(Point(x2_periodic, y_j));
 		}
+		if ((x - xLeft_simplest) + (xRight_simplest - x_j) <= cutoff_simplest) {
+			double x2_periodic = xLeft_simplest - xRight_simplest + x_j;
+			ped2->SetPos(Point(x2_periodic, y_j));
+		}
 	}
 
 	Point distp12 = ped2->GetPos() - ped1->GetPos();
@@ -745,6 +808,7 @@ Point SimplestModel::ForceRepPed(Pedestrian* ped1, Pedestrian* ped2, Point e0, i
 		R_ij = -_aPed * exp((-dist) / _DPed);
 		F_rep = ep12 * R_ij;
 	}
+	ped2->SetPos(Point(x_j, y_j));
 	return F_rep;
 }//END Velocity:ForceRepPed()