diff --git a/Data_structure.py b/Data_structure.py
new file mode 100755
index 0000000000000000000000000000000000000000..f1741d28e99a8594646a4e1a8e6bae95675b9be8
--- /dev/null
+++ b/Data_structure.py
@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun May 12 19:51:18 2019
+
+@author: MAC
+"""
+
+import Engineering_tools as et
+
+class Rocket():
+    """
+    """
+    
+    def __init__(self,stage1,stage2,payload,launch_coordinates):
+        """
+        stage1 and stage2 are instances of the class Stage,
+        payload is an instance of the class Payload,
+        launch_coordiantes = {latitude(deg),logitude(deg)}
+        """
+        
+        self.stage1 = stage1
+        self.stage2 = stage2
+        self.payload = payload
+        self.injection_velocity = et.final_target_vel(self.payload.perigee_altitude,
+                                                      self.payload.target_orbit['a'])
+        self.launch_coordinates = launch_coordinates
+        if self.payload.target_orbit['i'] > 90: # retrograade orbit
+            self.initial_velocity = et.initial_vel(self.launch_coordinates['latitude'],retrograde=True)
+        else:
+            self.initial_velocity = et.initial_vel(self.launch_coordinates['latitude'])
+        
+    def target_Vc(self,losses):
+        """
+        Compute the target critical velocity
+        losses : velocity losses approximation (m.s-1)
+        """
+        
+        self.target_Vc = self.injection_velocity-self.initial_velocity+losses
+        
+    
+class Stage():
+    """
+    """
+    
+    def __init__(self,stage_number,number_of_engines):
+        """
+        stage_number : 1 or 2
+        
+        """
+        self.stage_number = stage_number
+                
+    
+    avionics_mass = 0
+    
+    number_of_engines = 9
+    
+    def tank_mass(self,)
+        
+class Payload():
+    """
+    """
+    
+    def __init__(self,mass,target_orbit):
+        """
+        mass in kg
+        target_orbit = {'a','e','i','capital_omega','small_omega','nu'}
+        where:
+            a : semi_major axis (m)
+            e : eccentricity (adimentional)
+            i : inclination (deg)                                
+            capital_omega : argument of the ascending node (deg) 
+            small_omega   : argument of the perigee (deg)
+            nu : true anomaly (deg)
+        """
+        
+        self.mass = mass
+        self.target_orbit = target_orbit
+        self.perigee_altitude = target_orbit['a']*(1-target_orbit['e'])
+        
+#    def propagate_orbit():
+#        """
+#        """
+    
+class Engine():
+    """
+    """
+    
+    cycle_type=''
+    
+    def __init__(self,oxidiser,fuel,thermo_cycle,exit_section_area):
+        """
+        """
+        self.oxidiser = oxidiser
+        self.fuel = fuel
+        self.thermo_cycle = thermo_cycle
+        self.exit_section_area = exit_section_area
\ No newline at end of file
diff --git a/Engineering_tools.py b/Engineering_tools.py
new file mode 100755
index 0000000000000000000000000000000000000000..21d3397d92b20ca51def1c26af741ed51814870d
--- /dev/null
+++ b/Engineering_tools.py
@@ -0,0 +1,395 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun May 12 20:08:52 2019
+
+@author: MAC
+"""
+
+import numpy as np
+from numpy import linalg as LA
+    
+def initial_vel(latitude,retrograde=False):
+    """
+    Initial velocity (m.s-1) depending on the launch point latitude (°)
+    """
+    
+    T_side = 86164.10   # sidereal day (s)
+    R_T = 6378e3        # earth raduis (m)
+    
+    initial_velocity = (2*np.pi/T_side)*R_T*np.cos(np.deg2rad(np.abs(latitude)))
+    if retrograde:
+        return -initial_velocity
+    else:
+        return initial_velocity
+    
+def final_target_vel(h, a):
+    """
+    inputs:
+        h: injection altitude (m)
+        a: orbit semi major axis (m)
+    output: 
+        orbit absolute injection velocity (m.s-1)
+    """
+    
+    mu = 3.9860e14      # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
+    R_T = 6378e3        # earth raduis (m)
+    
+    return np.sqrt(2*mu*((1/(R_T+h))-(1/(2*a))))
+    
+def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp):
+    """
+    Optimise mass loading in the case of two stage rocket
+    inputs:
+        isp1, isp2  : specific impulse of the 1st and 2nd stage (s)
+        eps1, eps2  : structural coefficient of the 1st and 2nd stage (adimentional)
+        Vc          : target critical velocity (m.s-1)
+        Mp          : payload mass (Kg)
+    outputs:
+        M1  : mass of the 1st stage
+        M2  : mass of the 2nd stage
+        
+    NB: 
+                                    structure mass
+    structural coefficient = --------------------------------  
+                             structure mass + propellant mass
+                             
+    ctritical velocity = final velocity - initial velocity + losses
+    """
+    
+#    Notation :
+#        
+#        Mt = total mass
+#        q1 = M1/Mt
+#        q2 = M2/Mt
+#        qp = Mp/Mt
+#        omega = qp + q2
+    
+    g0 = 9.81 # standard gravity (m.s-2)
+    
+    omega_tab=np.linspace(0.001,0.9,100)  # array of omega to test
+    qp_tab=[]  
+    for omega in omega_tab:
+        q2 = omega+eps1*(1-omega)
+        q2 = np.power(q2,isp1/isp2)*np.exp(Vc/(g0*isp2))
+        q2 = 1/q2
+        q2 = 1-q2
+        q2 *= omega/(1-eps2)
+        qp = omega-q2
+        qp_tab.append(qp)
+    
+    qp=max(qp_tab)
+    omega=omega_tab[qp_tab.index(qp)]
+    q2=omega-qp
+    q1=1-q2-qp
+    M1=(Mp/qp)*q1
+    M2=(Mp/qp)*q2
+    
+    return M1, M2
+    
+def standard_atmosphere(altitude):
+    """
+    Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model    
+    input:
+        altitude (m)
+    outputs:
+        temperature (K)
+        pressure (Pa)
+        density (Kg.m-3)
+    """
+    
+    assert altitude < 84852 , str(altitude)+' is higher than max altitude of the ISA model' 
+    
+    g0 = 9.81 # standard gravity (m.s-2)
+    r = 287 # specific gas constant (J.Kg-1.K-1)
+    # ISA model
+    h0_tab = np.array([0,      11e3,    20e3,    32e3,    47e3,    51e3,    71e3,    84852]) # base altitude (m)
+    a_tab = np.array([-6.5e-3, 0,       1e-3,    2.8e-3,  0,       -2.8e-3, -2.0e-3]) # lapse rate (K.m-1)
+    # the base properties of each range are tabulated to reduce computational cost
+    T0_tab = np.array([288.15, 216.65,  216.65,  228.65,  270.65,  270.65,  214.65]) # base temperature (K)
+    P0_tab = np.array([101325, 22632,   5474.9,  868.02,  110.91,  66.939,  3.9564]) # base pressure (Pa)
+    
+    index = 0
+    while altitude > h0_tab[index+1]:
+        index += 1
+    
+    a = a_tab[index]
+    T0 = T0_tab[index]
+    P0 = P0_tab[index]
+    T = T0+a*(altitude-h0_tab[index])     # Toussaint equation
+    
+    if T == T0 : # isohermal region
+        P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T))
+    else : # gradient region
+        P = P0*np.power(T/T0,-g0/(a*r))
+    
+    rho = P/(r*T)    
+    
+    return T,P,rho
+    
+def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5):
+    """
+    Compute the thrust assuming a started nozzle
+    inputs:
+        Pt      : total pressure in the combustion chambre (Pa)
+        Tt      : total temperature in the combustion chambre (K)
+        Pamb    : ambient pressure (Pa)
+        Ath     : throat section area (m2)
+        Ae      : exit section area (m2)
+        gamma   : ratio of specific heats (adimentional)
+        r       : specific gas constant (J.Kg-1.K-1)
+        tolerance : 
+    outputs:
+        Thrust (N)
+        Specific impulse (s) 
+        
+    N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area
+    """
+    
+    assert gamma > 1, 'gamma must be > 1'
+    
+    g0 = 9.81 # standard gravity (m.s-2)
+    
+    A_ratio = lambda M,gamma: (1/M)*np.power((2/(gamma+1))*(1+((gamma-1)/2)*M*M),(gamma+1)/(2*gamma-2))
+    
+    # compute the nozzle exit section area corresponding to the adapted regime
+    Madapted = np.sqrt((2/(gamma-1))*(np.power(Pamb/Pt,(1-gamma)/gamma)-1))
+    Aadapted = Ath*A_ratio(Madapted,gamma)
+        
+    if Aadapted > Ae: # the nozzle is under-expanded
+        # compute the exit Mach number iteratively using dichotomy algorithm
+        A_ratio_target = Ae/Ath
+        Me_upper = Madapted
+        Me_lower = 1
+        Me = (Me_upper+Me_lower)/2
+        error = A_ratio_target-A_ratio(Me,gamma)
+        while np.abs(error) > tolerance:
+            if  error > 0:
+                Me_lower = Me
+            else:
+                Me_upper = Me
+            Me = (Me_upper+Me_lower)/2
+            error = A_ratio_target-A_ratio(Me,gamma)
+        Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
+        Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure
+        
+    elif Aadapted < Ae: # the nozzle is over-expanded
+        Mshock = 1
+        step = 1
+        A_ratio_target = Ae/Ath
+        Me_upper = 1
+        Me_lower = 0
+        Me = (Me_upper+Me_lower)/2
+        error = A_ratio_target-A_ratio(Me,gamma)
+        while np.abs(error) > tolerance:
+            if  error < 0:
+                Me_lower = Me
+            else:
+                Me_upper = Me
+            Me = (Me_upper+Me_lower)/2
+            error = A_ratio_target-A_ratio(Me,gamma)
+        Te = Tt/(1+((gamma-1)/2)*Me*Me)
+        Pe = Pt*np.power(Tt/Te,gamma/(1-gamma))
+        while (Pe-Pamb)*1e-5 > tolerance:
+            Mshock_old = Mshock
+            Mshock += step
+            Pe_old = Pe
+            Pt_ratio = np.power(((gamma+1)*Mshock*Mshock)/((gamma-1)*Mshock*Mshock+2),gamma/(gamma-1))*np.power((gamma+1)/(2*gamma*Mshock*Mshock-gamma+1),1/(gamma-1))
+            Pi = Pt*Pt_ratio
+            Ai = Ath/Pt_ratio
+            A_ratio_target = Ae/Ai
+            Me_upper = 1
+            Me_lower = 0
+            Me = (Me_upper+Me_lower)/2
+            error = A_ratio_target-A_ratio(Me,gamma)
+            while np.abs(error) > tolerance:
+                if  error < 0:
+                    Me_lower = Me
+                else:
+                    Me_upper = Me
+                Me = (Me_upper+Me_lower)/2
+                error = A_ratio_target-A_ratio(Me,gamma)
+            Te = Tt/(1+((gamma-1)/2)*Me*Me)
+            Pe = Pi*np.power(Tt/Te,gamma/(1-gamma))
+            if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae/Ath):
+                # the normal shock is outside the divergent and all the flow inside is supersonic
+                A_ratio_target = Ae/Ath
+                Me_upper = Mshock
+                Me_lower = 1
+                Me = (Me_upper+Me_lower)/2
+                error = A_ratio_target-A_ratio(Me,gamma)
+                while np.abs(error) > tolerance:
+                    if  error > 0:
+                        Me_lower = Me
+                    else:
+                        Me_upper = Me
+                    Me = (Me_upper+Me_lower)/2
+                    error = A_ratio_target-A_ratio(Me,gamma)
+                Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
+                Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure
+                break
+            if Pe < Pamb:
+                Pe = Pe_old
+                Mshock = Mshock_old
+                step /= 2
+        
+    else: # the nozzle is adapted
+        Me = Madapted
+        Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature
+        Pe = Pamb
+    
+    Ve = Me*np.sqrt(gamma*r*Te) # exit velocity
+    
+    # maximum mass flow rate in the case of a started nozzle
+    Q = np.power(2/(gamma+1),(gamma+1)/(2*(gamma-1)))*np.sqrt(gamma/(Tt*r))*Pt*Ath
+    Thrust = Q*Ve+(Pe-Pamb)*Ae
+    isp = Thrust/(Q*g0)
+
+    return Thrust,isp
+    
+def adiabatic_flame(mech,T,P,X):
+    """
+    Compute adiatabic flame temperature
+    inputs:
+        mech : chemical mechanism in cti format
+        T : initial temperature (K)
+        P : pressure (Pa)
+        X : initial chemical composition
+    """
+    # work in progress
+    
+
+#%% Trajectory tools
+
+def solve_traj():
+    """
+    """
+    
+# .......
+    
+#%% Astrodynamics tools
+    
+def state_to_orbital_elem(r,v):
+    """
+    Convert state vector to orbital elements
+    inputs:
+        r : coordinate vector (m)
+        v : velocity vector (m.s-1)
+    output: {'a','e','i','capital_omega','small_omega','nu'}
+    where:
+        a : semi_major axis (m)
+        e : eccentricity (adimentional)
+        i : inclination (deg)                                
+        capital_omega : argument of the ascending node (deg) 
+        small_omega   : argument of the perigee (deg)
+        nu : true anomaly (deg)
+    NB: - An internal quantity (small) is considered for the case of undefined 
+          classical orbital elements. When an element is undefined it is set to
+          zero and other(s) following it are measured according to this convention
+        - The function stops when r and v are aligned
+        - This function is widely inspired by the Matlab function written by 
+          Prof. Josep Masdemart from UPC 
+    """
+    
+    orbel={'a':0,'e':0,'i':0,'capital_omega':0,'small_omega':0,'nu':0}
+    
+    mu = 3.9860e14      # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2)
+    small = 1e-10
+    
+    h = np.cross(r,v) # angular momentum
+    rm = LA.norm(r) # module of the coordinate vector
+    vm = LA.norm(v) # module of the velocity vector
+    hm = LA.norm(h) # module of the angular momentum vector
+    if hm < small:
+        print('rectilinar or collision trajectory')
+        return
+    n = np.array([-h[1],h[0],0])
+    nm = LA.norm(n)
+    if nm < small:
+        n = np.array([1,0,0])
+        nm = 0
+    else:
+        n = n/nm
+        nm = 1
+    dotrv = np.dot(r,v)
+    e = ((vm*vm-mu/rm)*r-dotrv*v)/mu
+    em = LA.norm(e)
+    if em < small:
+        e = n
+        em = 0
+    else:
+        e = e/em
+    orbel['e'] = em
+    orbel['a'] = (hm*hm/mu)/(1-em*em)
+    aux = h[2]/hm
+    if np.abs(aux) > 1: aux = np.sign(aux)
+    # all the angles are first computed in rad
+    orbel['i'] = np.arccos(aux)
+    if nm > small:
+        aux = n[0]
+        if np.abs(aux) > 1: aux = np.sign(aux)
+        orbel['capital_omega'] = np.arccos(aux)
+    if n[1] < 0: orbel['capital_omega'] = 2*np.pi-orbel['capital_omega']
+    if em > small:
+        aux = np.dot(e,n)
+        if np.abs(aux) > 1: aux = np.sign(aux)
+        orbel['small_omega'] = np.arccos(aux)
+        if e[2] < 0: orbel['small_omega'] = 2*np.pi-orbel['small_omega']
+    aux = np.dot(e,r)/rm
+    if np.abs(aux) > 1: aux = np.sign(aux)
+    orbel['nu'] = np.arccos(aux)
+    if em > small:
+        if dotrv < 0: orbel['nu'] = 2*np.pi-orbel['nu']
+    else:
+        ha = np.cross(e,r)
+        if np.dot(h,ha) < 0: orbel['nu'] = 2*np.pi-orbel['nu']
+    # convert the angles to deg
+    orbel['i'] = np.rad2deg(orbel['i'])
+    orbel['capital_omega'] = np.rad2deg(orbel['capital_omega'])
+    orbel['small_omega'] = np.rad2deg(orbel['small_omega'])
+    orbel['nu'] = np.rad2deg(orbel['nu'])
+    return orbel
+
+def propagate_orbit():
+    """
+    """
+
+#%% Thermodynamic cycles 
+# the following functions provide the pressure, the mixuture composition 
+# and the initial temperature in the combustion chambre and properties of some 
+# components of the cycle.
+
+def pressure_fed():
+    """
+    """
+    
+def electric_turbo():
+    """
+    """
+    
+def open_expander():
+    """
+    """
+    
+def closed_expander():
+    """
+    """
+    
+def dual_expander():
+    """
+    """
+    
+def tap_off():
+    """
+    """
+    
+def gas_generator():
+    """
+    """
+    
+def ox_rich_staged_combust():
+    """
+    """
+    
+def full_flow_staged_combust():
+    """
+    """
diff --git a/test_laval_nozzle.py b/test_laval_nozzle.py
new file mode 100755
index 0000000000000000000000000000000000000000..a6d925a5801d5644018120f6d7f0222e59b1f4bb
--- /dev/null
+++ b/test_laval_nozzle.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed May 29 10:35:27 2019
+
+@author: MAC
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from Engineering_tools import laval_nozzle
+
+Pt = 20e5   # total pressure in the combustion chambre (Pa)
+Tt = 2000   # total temperature in the combustion chambre (K)
+Ath = 0.4   # throat section area (m2)
+Ae = 1.6    # exit section area (m2)
+gamma = 1.4 # ratio of specific heats
+r = 287     # specific gas constant (J.Kg-1.K-1)
+
+data = pd.DataFrame(columns=['Pamb','Thrust','isp'])
+Pamb_tab = np.linspace(1e5,100,100) # array of the ambient pressure (Pa)
+
+for Pamb in Pamb_tab:
+    data.loc[len(data)] = (Pamb,)+laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r)
+    
+plt.figure('Laval Nozzle Thrust')
+plt.plot(data['Pamb']*1e-5,data['Thrust']*1e-3)
+plt.title('Thrust Vs Ambient Pressure')
+plt.xlabel('Ambient Pressure (bar)')
+plt.ylabel('Thrust (kN)')
+
+plt.figure('Laval Nozzle Isp')
+plt.plot(data['Pamb']*1e-5,data['isp'])
+plt.title('Isp Vs Ambient Pressure')
+plt.xlabel('Ambient Pressure (bar)')
+plt.ylabel('Isp (s)')
diff --git a/test_standard_atmosphere.py b/test_standard_atmosphere.py
new file mode 100755
index 0000000000000000000000000000000000000000..2ca63bd160bc91e83378a0c6558b4d68f0b6cfad
--- /dev/null
+++ b/test_standard_atmosphere.py
@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun May 26 13:40:33 2019
+
+@author: MAC
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from Engineering_tools import standard_atmosphere
+
+altitude = np.linspace(0,80e3,801)
+data = pd.DataFrame(columns=['altitude','T','P','rho'])
+
+for h in altitude:
+    data.loc[len(data)] = (h,)+standard_atmosphere(h)
+    
+plt.figure('ISA Temperature')
+plt.plot(data['T'],data['altitude']*1e-3)
+plt.title('Temperature Vs Altitude')
+plt.xlabel('Temperature (K)')
+plt.ylabel('Altitude (km)')
+
+plt.figure('ISA Pressure')
+plt.plot(data['P'],data['altitude']*1e-3)
+plt.xscale('log')
+plt.title('Pressure Vs Altitude')
+plt.xlabel('Pressure (Pa)')
+plt.ylabel('Altitude (km)')
+
+plt.figure('ISA Density')
+plt.plot(data['rho'],data['altitude']*1e-3)
+plt.xscale('log')
+plt.title('Density Vs Altitude')
+plt.xlabel('Density (Kg.m-3)')
+plt.ylabel('Altitude (km)')
\ No newline at end of file