diff --git a/requirements.txt b/requirements.txt
index ede9c66..6c6d661 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -6,4 +6,4 @@ isort
 black
 pre-commit
 fmpy==0.3.15
-cosapp_fmu @ git+ssh://git@github.com/twiinIT/cosapp-fmu.git@master
+#cosapp_fmu @ git+ssh://git@github.com/twiinIT/cosapp-fmu.git@master
diff --git a/rocket-twin b/rocket-twin
new file mode 160000
index 0000000..105bfcc
--- /dev/null
+++ b/rocket-twin
@@ -0,0 +1 @@
+Subproject commit 105bfcc6a87803b023b21871fc3e4f0794d9aad5
diff --git a/rocket_twin/notebooks/Demonstration.ipynb b/rocket_twin/notebooks/Demonstration.ipynb
index 6901dcf..0d84243 100644
--- a/rocket_twin/notebooks/Demonstration.ipynb
+++ b/rocket_twin/notebooks/Demonstration.ipynb
@@ -23,7 +23,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from rocket_twin.systems import Station, StationControllerFMU, RocketControllerFMU, StageControllerFMU\n",
+    "from rocket_twin.systems import Station\n",
     "from rocket_twin.drivers import FuelingRocket, VerticalFlyingRocket, Mission\n",
     "from rocket_twin.utils import run_sequences\n",
     "from cosapp.utils import swap_system\n",
@@ -596,7 +596,9 @@
    "source": [
     "## Choice of control mode\n",
     "\n",
-    "The control of the system can be done through either cosapp or .fmu files. In order to use the fmu control, either a .mo file (from OpenModelica) or an .fmu file should be passed to the system."
+    "The control of the system can be done through either cosapp or .fmu files. In order to use the fmu control, either a .mo file (from OpenModelica) or an .fmu file should be passed to the system.\n",
+    "\n",
+    "In this release, only the cosapp control is available. "
    ]
   },
   {
@@ -606,22 +608,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "use_fmu = False\n",
-    "\n",
-    "model_path = r\"systems\\control\\station_controller.mo\"\n",
     "model_name = \"station_controller\"\n",
-    "\n",
-    "model_path_r = r\"systems\\control\\rocket_controller.mo\"\n",
     "model_name_r = \"rocket_controller\"\n",
-    "\n",
-    "model_path_s = r\"systems\\control\\stage_controller.mo\"\n",
-    "model_name_s = \"stage_controller\"\n",
-    "\n",
-    "if use_fmu:\n",
-    "    swap_system(sys.controller, StationControllerFMU(\"controller\", model_path=model_path, model_name=model_name))\n",
-    "    swap_system(sys.rocket.controller, RocketControllerFMU(\"controller\", model_path=model_path_r, model_name=model_name_r, n_stages=n_stages))\n",
-    "    for i in range(1, n_stages + 1):\n",
-    "        swap_system(sys.rocket[f\"stage_{i}\"].controller, StageControllerFMU(\"controller\", model_path=model_path_s, model_name=model_name_s))"
+    "model_name_s = \"stage_controller\""
    ]
   },
   {
@@ -1020,7 +1009,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.10.5"
   }
  },
  "nbformat": 4,
diff --git a/rocket_twin/systems/__init__.py b/rocket_twin/systems/__init__.py
index e65d754..ac4ba8f 100644
--- a/rocket_twin/systems/__init__.py
+++ b/rocket_twin/systems/__init__.py
@@ -1,19 +1,16 @@
 from rocket_twin.systems.control import (
     RocketControllerCoSApp,
-    RocketControllerFMU,
     StageControllerCoSApp,
-    StageControllerFMU,
     StationControllerCoSApp,
-    StationControllerFMU,
 )
 from rocket_twin.systems.engine import Engine, EngineGeom, EnginePerfo
-from rocket_twin.systems.ground import Ground
 from rocket_twin.systems.physics import Dynamics
 from rocket_twin.systems.structure import NoseGeom, TubeGeom, WingsGeom
 from rocket_twin.systems.tank import Pipe, Tank, TankFuel, TankGeom
 
 from rocket_twin.systems.rocket import OCCGeometry, Stage, Rocket  # isort: skip
 from rocket_twin.systems.station import Station  # isort: skip
+from rocket_twin.systems.ground import Ground  # isort: skip
 
 __all__ = [
     "Engine",
@@ -31,9 +28,6 @@
     "StageControllerCoSApp",
     "StationControllerCoSApp",
     "RocketControllerCoSApp",
-    "StageControllerFMU",
-    "StationControllerFMU",
-    "RocketControllerFMU",
     "NoseGeom",
     "TubeGeom",
     "WingsGeom",
diff --git a/rocket_twin/systems/control/__init__.py b/rocket_twin/systems/control/__init__.py
index bfd8c01..a182373 100644
--- a/rocket_twin/systems/control/__init__.py
+++ b/rocket_twin/systems/control/__init__.py
@@ -2,15 +2,8 @@
 from rocket_twin.systems.control.stage_controller_cosapp import StageControllerCoSApp
 from rocket_twin.systems.control.station_controller_cosapp import StationControllerCoSApp
 
-from rocket_twin.systems.control.rocket_controller_fmu import RocketControllerFMU  # isort: skip
-from rocket_twin.systems.control.stage_controller_fmu import StageControllerFMU  # isort: skip
-from rocket_twin.systems.control.station_controller_fmu import StationControllerFMU  # isort: skip
-
 __all__ = [
     "StageControllerCoSApp",
     "StationControllerCoSApp",
     "RocketControllerCoSApp",
-    "StageControllerFMU",
-    "StationControllerFMU",
-    "RocketControllerFMU",
 ]
diff --git a/rocket_twin/systems/control/rocket_controller.mo b/rocket_twin/systems/control/rocket_controller.mo
deleted file mode 100644
index ee9c50d..0000000
--- a/rocket_twin/systems/control/rocket_controller.mo
+++ /dev/null
@@ -1,34 +0,0 @@
-model rocket_controller
-  input Boolean flying;    // Whether the rocket is mid-flight or not
-  input Real weight_1;     // Stage 1 fuel mass
-  input Real weight_2;     // Stage 2 fuel mass
-  input Real weight_3;     // Stage 3 fuel mass
-  output Real is_on_1;     // Whether the stage 1 is on or not
-  output Real is_on_2;     // Whether the stage 2 is on or not
-  output Real is_on_3;     // Whether the stage 3 is on or not
-equation
-  if (flying == false) then
-    is_on_1 = 0.;
-    is_on_2 = 0.;
-    is_on_3 = 0.;
-  else
-    if (weight_1 >= 0.1) then
-      is_on_1 = 1.;
-      is_on_2 = 0.;
-      is_on_3 = 0.;
-    elseif (weight_2 >= 0.1) then
-      is_on_1 = 0.;
-      is_on_2 = 1.;
-      is_on_3 = 0.;
-    elseif (weight_3 >= 0.1) then
-      is_on_1 = 0.;
-      is_on_2 = 0.;
-      is_on_3 = 1.;
-    else
-      is_on_1 = 0.;
-      is_on_2 = 0.;
-      is_on_3 = 0.;
-    end if;
-  end if;
-
-end rocket_controller;
diff --git a/rocket_twin/systems/control/rocket_controller/rocket_controller.fmu b/rocket_twin/systems/control/rocket_controller/rocket_controller.fmu
index 217677e..8bf28a5 100644
Binary files a/rocket_twin/systems/control/rocket_controller/rocket_controller.fmu and b/rocket_twin/systems/control/rocket_controller/rocket_controller.fmu differ
diff --git a/rocket_twin/systems/control/rocket_controller_fmu.py b/rocket_twin/systems/control/rocket_controller_fmu.py
deleted file mode 100644
index 68e8623..0000000
--- a/rocket_twin/systems/control/rocket_controller_fmu.py
+++ /dev/null
@@ -1,95 +0,0 @@
-import os
-
-from cosapp.base import System
-from cosapp_fmu.FMUsystem import FMUSystem
-from OMPython import ModelicaSystem
-
-import rocket_twin.systems.control
-
-
-class RocketControllerFMU(System):
-    """Controller of the command variables.
-
-    Inputs
-    ------
-    model_path: string,
-        the path to the .mo file, if any
-    model_name: string,
-        the .fmu file name
-    flying: boolean,
-        whether the rocket is mid-flight or not
-    n_stages: int,
-        rocket's number of stages
-    weight_prop_i: float,
-        i-th stage fuel weight
-
-    Outputs
-    ------
-    is_on_i: boolean,
-        whether the i-th stage controller is active or not
-    """
-
-    def setup(self, n_stages, model_path, model_name):
-
-        self.add_inward("n_stages", n_stages, desc="number of stages")
-        self.add_inward("stage", 1, desc="Current active stage")
-
-        pulling = {"flying": "flying"}
-
-        for i in range(1, n_stages + 1):
-            self.add_outward(f"is_on_{i}", 0, desc=f"Whether the stage {i} is on or not")
-            pulling[f"weight_{i}"] = f"weight_prop_{i}"
-            pulling[f"is_on_{i}"] = f"is_on_{i}"
-
-        fmu_path = self.create_fmu(model_path, model_name)
-        self.add_child(
-            FMUSystem("fmu_controller", fmu_path=fmu_path),
-            pulling=pulling,
-        )
-
-        self.add_event("drop", trigger="weight_prop_1 < 0.1")
-
-    def compute(self):
-
-        for i in range(1, self.n_stages):
-            self[f"is_on_{i}"] = bool(self[f"is_on_{i}"])
-
-    def transition(self):
-
-        if self.drop.present:
-            if self.stage < self.n_stages:
-                self.stage += 1
-                self.drop.trigger = f"weight_prop_{self.stage} < 0.1"
-
-    def create_fmu(self, model_path, model_name):
-        """Create an fmu file in the control folder from an mo file.
-
-        Inputs
-        ------
-        model_path: string
-            the path of the .mo file
-        model_name: string
-            the name of the .fmu file to be created
-
-        Outputs
-        ------
-
-        fmu: string
-            the path to the .fmu file
-        """
-
-        fmu_path = os.path.join(rocket_twin.systems.control.__path__[0], model_name)
-        model_path = os.path.join(rocket_twin.__path__[0], model_path)
-        model_path = model_path.replace("\\", "/")
-        try:
-            os.mkdir(fmu_path)
-        except OSError:
-            pass
-        os.chdir(fmu_path)
-        mod = ModelicaSystem(model_path, model_name)
-        fmu = mod.convertMo2Fmu()
-        for filename in os.listdir(fmu_path):
-            if filename != (model_name + ".fmu"):
-                os.remove(filename)
-
-        return fmu
diff --git a/rocket_twin/systems/control/stage_controller.mo b/rocket_twin/systems/control/stage_controller.mo
deleted file mode 100644
index c0fd11f..0000000
--- a/rocket_twin/systems/control/stage_controller.mo
+++ /dev/null
@@ -1,10 +0,0 @@
-model stage_controller
-  input Boolean is_on;   // Whether the controller is on or not
-  output Real w;         // Command flow
-equation
-  if (is_on == true) then
-    w = 1.;
-  else
-    w = 0.;
-  end if;
-end stage_controller;
diff --git a/rocket_twin/systems/control/stage_controller/stage_controller.fmu b/rocket_twin/systems/control/stage_controller/stage_controller.fmu
index 2fdaadd..20979fc 100644
Binary files a/rocket_twin/systems/control/stage_controller/stage_controller.fmu and b/rocket_twin/systems/control/stage_controller/stage_controller.fmu differ
diff --git a/rocket_twin/systems/control/stage_controller_cosapp.py b/rocket_twin/systems/control/stage_controller_cosapp.py
index 8724095..c6af650 100644
--- a/rocket_twin/systems/control/stage_controller_cosapp.py
+++ b/rocket_twin/systems/control/stage_controller_cosapp.py
@@ -27,7 +27,7 @@ def setup(self):
 
         self.add_outward("w", 1.0, desc="Command flow", unit="")
 
-        self.add_event("full", trigger="weight_prop == weight_max")
+        self.add_event("full", trigger="weight_prop == weight_max") # transition is defined in station
 
     def compute(self):
 
diff --git a/rocket_twin/systems/control/stage_controller_fmu.py b/rocket_twin/systems/control/stage_controller_fmu.py
deleted file mode 100644
index 73d783c..0000000
--- a/rocket_twin/systems/control/stage_controller_fmu.py
+++ /dev/null
@@ -1,77 +0,0 @@
-import os
-
-from cosapp.base import System
-from cosapp_fmu.FMUsystem import FMUSystem
-from OMPython import ModelicaSystem
-
-import rocket_twin.systems.control
-
-
-class StageControllerFMU(System):
-    """Controller of the command variables.
-
-    Inputs
-    ------
-    model_path: string,
-        the path to the .mo file, if any
-    model_name: string,
-        the .fmu file name
-    is_on: boolean,
-        whether the system is in fueling phase or not
-    weight_prop: float,
-        stage fuel weight
-    weight_max: float,
-        stage maximum fuel weight
-
-    Outputs
-    ------
-    w: float,
-        command flux
-    """
-
-    def setup(self, model_path, model_name):
-
-        self.add_inward("weight_prop", 0.0, desc="Stage propellant weight", unit="kg")
-        self.add_inward("weight_max", 1.0, desc="Stage maximum propellant weight", unit="kg")
-        self.add_inward("is_on", False, desc="Whether the stage is on or not")
-
-        fmu_path = self.create_fmu(model_path, model_name)
-        self.add_child(
-            FMUSystem("fmu_controller", fmu_path=fmu_path),
-            pulling=["is_on", "w"],
-        )
-
-        self.add_event("full", trigger="weight_prop == weight_max")
-
-    def create_fmu(self, model_path, model_name):
-        """Create an fmu file in the control folder from an mo file.
-
-        Inputs
-        ------
-        model_path: string
-            the path of the .mo file
-        model_name: string
-            the name of the .fmu file to be created
-
-        Outputs
-        ------
-
-        fmu: string
-            the path to the .fmu file
-        """
-
-        fmu_path = os.path.join(rocket_twin.systems.control.__path__[0], model_name)
-        model_path = os.path.join(rocket_twin.__path__[0], model_path)
-        model_path = model_path.replace("\\", "/")
-        try:
-            os.mkdir(fmu_path)
-        except OSError:
-            pass
-        os.chdir(fmu_path)
-        mod = ModelicaSystem(model_path, model_name)
-        fmu = mod.convertMo2Fmu()
-        for filename in os.listdir(fmu_path):
-            if filename != (model_name + ".fmu"):
-                os.remove(filename)
-
-        return fmu
diff --git a/rocket_twin/systems/control/station_controller.mo b/rocket_twin/systems/control/station_controller.mo
deleted file mode 100644
index 54f50cb..0000000
--- a/rocket_twin/systems/control/station_controller.mo
+++ /dev/null
@@ -1,10 +0,0 @@
-model station_controller
-  input Boolean fueling;       // Whether the system is in its fueling phase
-  output Real w;               // Command flow
-equation
-  if (fueling == true) then
-    w = 1.;
-  else
-    w = 0.;
-  end if;
-end station_controller;
diff --git a/rocket_twin/systems/control/station_controller/station_controller.fmu b/rocket_twin/systems/control/station_controller/station_controller.fmu
index 76fee4e..e5366f6 100644
Binary files a/rocket_twin/systems/control/station_controller/station_controller.fmu and b/rocket_twin/systems/control/station_controller/station_controller.fmu differ
diff --git a/rocket_twin/systems/control/station_controller_fmu.py b/rocket_twin/systems/control/station_controller_fmu.py
deleted file mode 100644
index 096a029..0000000
--- a/rocket_twin/systems/control/station_controller_fmu.py
+++ /dev/null
@@ -1,67 +0,0 @@
-import os
-
-from cosapp.base import System
-from cosapp_fmu.FMUsystem import FMUSystem
-from OMPython import ModelicaSystem
-
-import rocket_twin.systems.control
-
-
-class StationControllerFMU(System):
-    """Controller of the command variables.
-
-    Inputs
-    ------
-    model_path: string,
-        the path to the .mo file, if any
-    model_name: string,
-        the .fmu file name
-    fueling: boolean,
-        whether the system is in fueling phase or not
-
-    Outputs
-    ------
-    w: float,
-        command flux
-    """
-
-    def setup(self, model_path, model_name):
-
-        fmu_path = self.create_fmu(model_path, model_name)
-        self.add_child(
-            FMUSystem("fmu_controller", fmu_path=fmu_path),
-            pulling=["fueling", "w"],
-        )
-
-    def create_fmu(self, model_path, model_name):
-        """Create an fmu file in the control folder from an mo file.
-
-        Inputs
-        ------
-        model_path: string
-            the path of the .mo file
-        model_name: string
-            the name of the .fmu file to be created
-
-        Outputs
-        ------
-
-        fmu: string
-            the path to the .fmu file
-        """
-
-        fmu_path = os.path.join(rocket_twin.systems.control.__path__[0], model_name)
-        model_path = os.path.join(rocket_twin.__path__[0], model_path)
-        model_path = model_path.replace("\\", "/")
-        try:
-            os.mkdir(fmu_path)
-        except OSError:
-            pass
-        os.chdir(fmu_path)
-        mod = ModelicaSystem(model_path, model_name)
-        fmu = mod.convertMo2Fmu()
-        for filename in os.listdir(fmu_path):
-            if filename != (model_name + ".fmu"):
-                os.remove(filename)
-
-        return fmu
diff --git a/rocket_twin/systems/engine/engine.py b/rocket_twin/systems/engine/engine.py
index 3203f7d..dd76f41 100644
--- a/rocket_twin/systems/engine/engine.py
+++ b/rocket_twin/systems/engine/engine.py
@@ -24,4 +24,4 @@ class Engine(System):
     def setup(self):
 
         self.add_child(EngineGeom("geom"), pulling=["shape", "props"])
-        self.add_child(EnginePerfo("perfo"), pulling=["w_out", "force"])
+        self.add_child(EnginePerfo("perfo"), pulling=["w_out", "force","v"])
diff --git a/rocket_twin/systems/engine/engine_perfo.py b/rocket_twin/systems/engine/engine_perfo.py
index 75d2867..e255b5b 100644
--- a/rocket_twin/systems/engine/engine_perfo.py
+++ b/rocket_twin/systems/engine/engine_perfo.py
@@ -1,6 +1,9 @@
+import numpy as np
+from numpy import linalg
 from cosapp.base import System
 
 
+
 class EnginePerfo(System):
     """Simple model of an engine's thrust force.
 
@@ -15,17 +18,18 @@ class EnginePerfo(System):
         thrust force
     """
 
-    def setup(self):
+    def setup(self, stations= None):
 
         # Inputs
         self.add_inward("w_out", 0.0, desc="Fuel consumption rate", unit="kg/s")
 
         # Parameters
+        self.add_inward("v", 1, desc = "velocity", unit= "m/s")
         self.add_inward("isp", 20.0, desc="Specific impulsion in vacuum", unit="s")
         self.add_inward("g_0", 10.0, desc="Gravity at Earth's surface", unit="m/s**2")
 
-        self.add_outward("force", 1.0, desc="Thrust force", unit="N")
+        self.add_outward("force", np.array([0.0, 0.0, 1.0]), desc="Thrust force", unit="N")
 
     def compute(self):
 
-        self.force = self.isp * self.w_out * self.g_0
+        self.force = (self.v/np.linalg.norm(self.v))*self.isp * self.w_out * self.g_0
\ No newline at end of file
diff --git a/rocket_twin/systems/ground.py b/rocket_twin/systems/ground.py
index 6e201b9..36c4f8e 100644
--- a/rocket_twin/systems/ground.py
+++ b/rocket_twin/systems/ground.py
@@ -1,4 +1,8 @@
 from cosapp.base import System
+from math import pi
+import numpy as np
+
+from rocket_twin.systems import Station
 
 
 class Ground(System):
@@ -18,6 +22,18 @@ def setup(self, stations=None):
         if stations is None:
             stations = []
 
+        self.add_property("stations", stations)  # to keep the station information in compute
+
         for station in stations:
 
-            self.add_child(station)
+            self.add_child(
+                Station(station, n_stages=1), pulling={"a": f"a_{station}","v":f"v_{station}", "pos": f"pos_{station}"}
+            )
+            self.add_transient(
+                f"v_{station}", der=f"a_{station}"
+            )  # integrate acceleration to get velocity
+            self.add_transient(
+                f"pos_{station}", der=f"v_{station}"
+            )  # integrate velocity to get position
+
+
diff --git a/rocket_twin/systems/physics/dynamics.py b/rocket_twin/systems/physics/dynamics.py
index 1932387..ed43130 100644
--- a/rocket_twin/systems/physics/dynamics.py
+++ b/rocket_twin/systems/physics/dynamics.py
@@ -1,3 +1,4 @@
+import numpy as np
 from cosapp.base import System
 
 
@@ -32,19 +33,29 @@ def setup(self, forces=None, weights=None):
         self.add_property("forces", forces)
         self.add_property("weights", weights)
 
-        self.add_inward("g", -10.0, desc="Gravity", unit="m/s**2")
+        # inwards
+        self.add_inward("pos", np.array([0.0, 0.0, 6400.0e3]), desc="local cartesian position", unit="m")
+        self.add_inward("g0", -10.0, desc="Gravity on earth", unit="m/s**2")
+        self.add_inward("r0", 6400.0e3, desc="Radius of the Earth", unit="m")
 
-        self.add_outward("a", 0.0, desc="Acceleration", unit="m/s**2")
+        # outwards
+        self.add_outward("g", np.array([0.0, 0.0, -10.0]), desc="Gravity", unit="m/s**2")
+        self.add_outward("a", np.zeros(3), desc="Acceleration", unit="m/s**2")
 
         for weight in self.weights:
             self.add_inward(weight, 0.0, desc=f"Weight called {weight}", unit="kg")
         for force in self.forces:
-            self.add_inward(force, 0.0, desc=f"Force called {force}", unit="N")
+            self.add_inward(force, np.zeros(3), desc=f"Force called {force}", unit="N")
 
-        self.add_outward("force", 1.0, desc="Force", unit="N")
+        self.add_outward("force", np.zeros(3), desc="Force", unit="N")
         self.add_outward("weight", 1.0, desc="Weight", unit="kg")
 
     def compute(self):
+        # gravity computation
+        r = np.linalg.norm(self.pos)
+        u = self.pos / r
+        self.g = self.g0 * (self.r0 / r) ** 2 * u
+
         self.weight = 0
 
         for weight in self.weights:
diff --git a/rocket_twin/systems/rocket/rocket.py b/rocket_twin/systems/rocket/rocket.py
index dc5c585..0ec85c4 100644
--- a/rocket_twin/systems/rocket/rocket.py
+++ b/rocket_twin/systems/rocket/rocket.py
@@ -1,5 +1,7 @@
+import numpy as np
 from cosapp.base import System
 from OCC.Core.GProp import GProp_GProps
+from math import pi
 
 from rocket_twin.systems import Dynamics, RocketControllerCoSApp
 from rocket_twin.systems.rocket import OCCGeometry, Stage
@@ -34,6 +36,11 @@ def setup(self, n_stages=1):
 
         shapes, properties, forces = ([None] * n_stages for i in range(3))
 
+        #parameter for rotating the rocket
+        self.add_inward("teta", 2*pi/180, desc="rocket rotation angle")
+        self.add_inward("Hm", 6500.0e3 , desc ="maneuvering altitude at which we want the rocket to rotate")
+        self.add_inward("rotate", True , desc ="condition for rotating")
+
         self.add_inward("n_stages", n_stages, desc="Number of stages")
         self.add_outward("stage", 1, desc="Current stage")
         self.add_inward("flying", False, desc="Whether the rocket is flying or not")
@@ -50,7 +57,7 @@ def setup(self, n_stages=1):
                 Stage(f"stage_{i}", nose=nose, wings=wings),
                 pulling={
                     "w_in": f"w_in_{i}",
-                    "weight_prop": f"weight_prop_{i}",
+                    "weight_prop": f"weight_prop_{i}", "v":"v"
                 },
             )
             shapes[i - 1] = f"stage_{i}_s"
@@ -63,7 +70,9 @@ def setup(self, n_stages=1):
             pulling=["flying"],
         )
         self.add_child(OCCGeometry("geom", shapes=shapes, properties=properties))
-        self.add_child(Dynamics("dyn", forces=forces, weights=["weight_rocket"]), pulling=["a"])
+        self.add_child(
+            Dynamics("dyn", forces=forces, weights=["weight_rocket"]), pulling=["a", "pos"]
+        )  # pulling acceleration and position to dynamics
 
         for i in range(1, n_stages + 1):
             self.connect(
@@ -79,15 +88,20 @@ def setup(self, n_stages=1):
 
         self.connect(self.geom.outwards, self.dyn.inwards, {"weight": "weight_rocket"})
 
+        self.add_event("rotation", trigger = "pos[2] >= Hm")  # if we reach the altitude Hm, we want to rotate the rocket
+
     def compute(self):
-        self.a *= self.flying
+        self.a *= self.flying #if the whole rocket is not flying, its acceleration is False = 0
 
-    def transition(self):
+    def transition(self):  
+        
+        if self.rotation.present:  # rotate the speed vector around y-axis by multiplying by rotation matrix at a given height (Hm)
+            self.v = np.dot(np.array([[np.cos(self.teta), 0, np.sin(self.teta)], [0, 1, 0], [-np.sin(self.teta), 0, np.cos(self.teta)]]),self.v.T ) 
 
-        if self.controller.drop.present:
+        if self.controller.drop.present: #change active stage when it runs out of propellant (mass = 0)
             if self.stage < self.n_stages:
                 stage = self.pop_child(f"stage_{self.stage}")
                 self.add_child(stage, execution_index=self.stage - 1)
                 self.geom[f"stage_{self.stage}"] = GProp_GProps()
-                self.dyn[f"thrust_{self.stage}"] = 0.0
+                self.dyn[f"thrust_{self.stage}"] = np.zeros(3)
                 self.stage += 1
diff --git a/rocket_twin/systems/rocket/stage.py b/rocket_twin/systems/rocket/stage.py
index a29cf6f..2ddbfe0 100644
--- a/rocket_twin/systems/rocket/stage.py
+++ b/rocket_twin/systems/rocket/stage.py
@@ -31,7 +31,7 @@ def setup(self, nose=False, wings=False):
 
         self.add_child(StageControllerCoSApp("controller"), pulling=["is_on"])
         self.add_child(Tank("tank"), pulling=["w_in", "weight_prop"])
-        self.add_child(Engine("engine"), pulling={"force": "thrust"})
+        self.add_child(Engine("engine"), pulling={"force": "thrust", "v":"v"})
         self.add_child(TubeGeom("tube"))
 
         if nose:
diff --git a/rocket_twin/systems/station/station.py b/rocket_twin/systems/station/station.py
index 7bbc846..fb9969b 100644
--- a/rocket_twin/systems/station/station.py
+++ b/rocket_twin/systems/station/station.py
@@ -33,7 +33,9 @@ def setup(self, n_stages=1):
         self.add_child(StationControllerCoSApp("controller"), pulling=["fueling"])
         self.add_child(Tank("g_tank"))
         self.add_child(Pipe("pipe"))
-        self.add_child(Rocket("rocket", n_stages=n_stages))
+        self.add_child(
+            Rocket("rocket", n_stages=n_stages), pulling=["a","v", "pos"]
+        )  # pulling acceleration and position (from ground to station to rocket to dynamics))
 
         self.connect(self.g_tank.outwards, self.pipe.inwards, {"w_out": "w_in"})
         self.connect(self.pipe.outwards, self.rocket.inwards, {"w_out": "w_in_1"})
diff --git a/rocket_twin/tests/test_controller_cosapp.py b/rocket_twin/tests/test_controller_cosapp.py
index 1920ec6..03e8960 100644
--- a/rocket_twin/tests/test_controller_cosapp.py
+++ b/rocket_twin/tests/test_controller_cosapp.py
@@ -40,4 +40,4 @@ def test_control(self):
         acel = np.asarray(data["rocket.a"])
 
         np.testing.assert_allclose(sys.rocket.geom.weight, 4.0, atol=10 ** (0))
-        np.testing.assert_allclose(acel[-3], 40.0, atol=0.1)
+        np.testing.assert_allclose(acel[-3], np.array([0.0, 0.0, 40.0]), atol=0.1)
diff --git a/rocket_twin/tests/test_controller_fmu.py b/rocket_twin/tests/test_controller_fmu.py
deleted file mode 100644
index 68ea07c..0000000
--- a/rocket_twin/tests/test_controller_fmu.py
+++ /dev/null
@@ -1,73 +0,0 @@
-import numpy as np
-from cosapp.drivers import NonLinearSolver, RungeKutta
-from cosapp.recorders import DataFrameRecorder
-from cosapp.utils import swap_system
-
-from rocket_twin.systems import (
-    RocketControllerFMU,
-    StageControllerFMU,
-    Station,
-    StationControllerFMU,
-)
-
-
-class TestControllerFMU:
-    """Tests for the FMU controller."""
-
-    def test_controller_fmu(self):
-
-        n_stages = 3
-
-        model_path = r"systems\control\station_controller.mo"
-        model_name = "station_controller"
-
-        model_path_r = r"systems\control\rocket_controller.mo"
-        model_name_r = "rocket_controller"
-
-        model_path_s = r"systems\control\stage_controller.mo"
-        model_name_s = "stage_controller"
-
-        sys = Station("sys", n_stages=n_stages)
-        swap_system(
-            sys.controller,
-            StationControllerFMU("controller", model_path=model_path, model_name=model_name),
-        )
-        swap_system(
-            sys.rocket.controller,
-            RocketControllerFMU(
-                "controller", n_stages=n_stages, model_path=model_path_r, model_name=model_name_r
-            ),
-        )
-        for i in range(1, n_stages + 1):
-            swap_system(
-                sys.rocket[f"stage_{i}"].controller,
-                StageControllerFMU("controller", model_path=model_path_s, model_name=model_name_s),
-            )
-
-        driver = sys.add_driver(RungeKutta(order=4, time_interval=[0, 35], dt=1.0))
-        driver.add_child(NonLinearSolver("solver"))
-        init = {
-            "g_tank.fuel.weight_p": 20.0,
-            "rocket.stage_1.tank.fuel.weight_p": 0.0,
-            "rocket.stage_2.tank.fuel.weight_p": 0.0,
-            "rocket.stage_3.tank.fuel.weight_p": 0.0,
-        }
-        values = {
-            "g_tank.fuel.w_out_max": 1.0,
-            "time_int": 5.0,
-            "rocket.stage_1.tank.fuel.w_out_max": 1.0,
-            "rocket.stage_2.tank.fuel.w_out_max": 1.0,
-            "rocket.stage_3.tank.fuel.w_out_max": 1.0,
-        }
-        driver.set_scenario(init=init, values=values)
-        driver.add_recorder(
-            DataFrameRecorder(includes=["rocket.a"]),
-            period=1.0,
-        )
-        sys.run_drivers()
-        data = driver.recorder.export_data()
-
-        acel = np.asarray(data["rocket.a"])
-
-        np.testing.assert_allclose(sys.rocket.geom.weight, 4.0, atol=10 ** (0))
-        np.testing.assert_allclose(acel[-2], 40.0, rtol=0.1)
diff --git a/rocket_twin/tests/test_drivers.py b/rocket_twin/tests/test_drivers.py
index af98b18..8e5e921 100644
--- a/rocket_twin/tests/test_drivers.py
+++ b/rocket_twin/tests/test_drivers.py
@@ -67,6 +67,6 @@ def test_flight(self):
 
         acel = np.asarray(data["rocket.a"])
 
-        np.testing.assert_allclose(acel[-2], 65.0, atol=10 ** (-10))
+        np.testing.assert_allclose(acel[-2], np.array([0.0, 0.0, 65.0]), atol=10 ** (-10))
         np.testing.assert_allclose(self.sys.rocket.stage_1.tank.weight_prop, 0.0, atol=10 ** (-10))
         np.testing.assert_allclose(self.sys.g_tank.weight_prop, 0.0, atol=10 ** (-10))
diff --git a/rocket_twin/tests/test_dynamics.py b/rocket_twin/tests/test_dynamics.py
index aee8258..8d347ec 100644
--- a/rocket_twin/tests/test_dynamics.py
+++ b/rocket_twin/tests/test_dynamics.py
@@ -8,18 +8,18 @@ class TestDynamics:
 
     def test_is_on(self):
         sys = Dynamics("sys", forces=["F"], weights=["w"])
-        sys.F = 100.0
+        sys.F = np.array([0.0, 0.0, 100.0])
         sys.w = 5.0
 
         sys.run_once()
 
-        np.testing.assert_allclose(sys.a, 10.0, atol=10 ** (-10))
+        np.testing.assert_allclose(sys.a, np.array([0.0, 0.0, 10.0]), atol=10 ** (-10))
 
     def test_is_off(self):
         sys = Dynamics("sys", forces=["F"], weights=["w"])
-        sys.F = 10.0
+        sys.F = np.array([0.0, 0.0, 10.0])
         sys.w = 5.0
 
         sys.run_once()
 
-        np.testing.assert_allclose(sys.a, -8.0, atol=10 ** (-10))
+        np.testing.assert_allclose(sys.a, np.array([0.0, 0.0, -8.0]), atol=10 ** (-10))
diff --git a/rocket_twin/tests/test_engine.py b/rocket_twin/tests/test_engine.py
index e4350da..2d1f53a 100644
--- a/rocket_twin/tests/test_engine.py
+++ b/rocket_twin/tests/test_engine.py
@@ -56,4 +56,10 @@ def test_perfo(self):
 
         sys.run_drivers()
 
-        np.testing.assert_allclose(sys.force, 2000.0, atol=10 ** (-5))
+        np.testing.assert_allclose(sys.force, np.array([0.0, 0.0, 2000.0]), atol=10 ** (-5))
+
+    def test_velocity(self):
+        sys = Engine("sys")
+        sys.run_once()
+        
+        np.testing.assert_allclose(sys.force/np.linalg.norm(sys.force), sys.v/np.linalg.norm(sys.v)) #test that velocity is parallel to thrust
diff --git a/rocket_twin/tests/test_mission.py b/rocket_twin/tests/test_mission.py
index f6b7323..8568028 100644
--- a/rocket_twin/tests/test_mission.py
+++ b/rocket_twin/tests/test_mission.py
@@ -33,6 +33,6 @@ def test_run_once(self):
 
         acel = np.asarray(data["rocket.a"])
 
-        np.testing.assert_allclose(acel[-2], 65.0, atol=10 ** (-10))
+        np.testing.assert_allclose(acel[-2], np.array([0.0, 0.0, 65.0]), atol=10 ** (-10))
         np.testing.assert_allclose(sys.rocket.stage_1.tank.weight_prop, 0.0, atol=10 ** (-10))
         np.testing.assert_allclose(sys.g_tank.weight_prop, 5.0, atol=10 ** (-10))
diff --git a/rocket_twin/tests/test_sequences.py b/rocket_twin/tests/test_sequences.py
index 737c7d8..d5789a4 100644
--- a/rocket_twin/tests/test_sequences.py
+++ b/rocket_twin/tests/test_sequences.py
@@ -81,7 +81,7 @@ def test_flight(self):
             atol=10 ** (-6),
         )
         np.testing.assert_allclose(self.sys.rocket.stage_1.tank.weight_prop, 0.0, atol=10 ** (-6))
-        np.testing.assert_allclose(acel[-2], 2.5, atol=10 ** (-6))
+        np.testing.assert_allclose(acel[-2], np.array([0.0, 0.0, 2.5]), atol=10 ** (-6))
 
     def test_all(self):
 
@@ -124,4 +124,4 @@ def test_all(self):
             atol=10 ** (-6),
         )
         np.testing.assert_allclose(sys2.rocket.stage_1.tank.weight_prop, 0.0, atol=10 ** (-6))
-        np.testing.assert_allclose(acel[-2], 2.5, atol=10 ** (-6))
+        np.testing.assert_allclose(acel[-2], np.array([0.0, 0.0, 2.5]), atol=10 ** (-6))
diff --git a/rocket_twin/tests/test_stage.py b/rocket_twin/tests/test_stage.py
index fc32cc7..44e41a5 100644
--- a/rocket_twin/tests/test_stage.py
+++ b/rocket_twin/tests/test_stage.py
@@ -66,6 +66,7 @@ def test_flight(self):
         stop = "rocket.weight_prop_3 == 0."
 
         includes = [
+            "rocket.a",
             "g_tank.weight_prop",
             "rocket.weight_prop_1",
             "rocket.weight_prop_2",
@@ -93,11 +94,11 @@ def test_flight(self):
         print(data1)
         print(data2)
 
-        # acel = np.asarray(data["rocket.a"])
+        acel = np.asarray(data["rocket.a"])
         # print(data)
 
         np.testing.assert_allclose(self.sys.rocket.weight_prop_1, 0.0, rtol=10 ** (-1))
         np.testing.assert_allclose(self.sys.rocket.weight_prop_2, 0.0, rtol=10 ** (-1))
         np.testing.assert_allclose(self.sys.rocket.weight_prop_3, 0.0, rtol=10 ** (-1))
         np.testing.assert_allclose(self.sys.rocket.geom.weight, 4.0, rtol=10 ** (-1))
-        # np.testing.assert_allclose(acel[-2], 40.0, atol=10 ** (-2))
+        np.testing.assert_allclose(acel[-2], np.array([0.0, 0.0, 40.0]), atol=10 ** (-2))