diff --git a/examples/ActivityDrivenInertial/conf.toml b/examples/ActivityDrivenInertial/conf.toml
new file mode 100644
index 0000000..383ba76
--- /dev/null
+++ b/examples/ActivityDrivenInertial/conf.toml
@@ -0,0 +1,34 @@
+[simulation]
+model = "ActivityDrivenInertial"
+# rng_seed = 120 # Leaving this empty will pick a random seed
+
+[io]
+n_output_network = 20 # Write the network every 20 iterations
+n_output_agents = 1   # Write the opinions of agents after every iteration
+print_progress = true # Print the iteration time ; if not set, then does not print
+
+[model]
+max_iterations = 500 # If not set, max iterations is infinite
+
+[ActivityDrivenInertial]
+dt = 0.01               # Timestep for the integration of the coupled ODEs
+m = 10                  # Number of agents contacted, when the agent is active
+eps = 0.01              # Minimum activity epsilon; a_i belongs to [epsilon,1]
+gamma = 2.1             # Exponent of activity power law distribution of activities
+reciprocity = 0.65      # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections
+homophily = 1.0         # aka beta. if zero, agents pick their interaction partners at random
+alpha = 3.0             # Controversialness of the issue, must be greater than 0.
+K = 2.0
+mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents
+mean_weights = false    # Use the meanfield approximation of the network edges
+
+reluctances = true       # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1
+reluctance_mean = 1.0    # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0)
+reluctance_sigma = 0.25  # Width of normal distribution (before truncating)
+reluctance_eps = 0.01    # Minimum such that the normal distribution is truncated at this value
+
+friction_coefficient = 1.0 # Friction coefficient, making agents tend to go to rest without acceleration
+
+[network]
+number_of_agents = 1000
+connections_per_agent = 10
diff --git a/examples/ActivityDrivenReluctance/conf.toml b/examples/ActivityDrivenReluctance/conf.toml
index 668891b..52b1005 100644
--- a/examples/ActivityDrivenReluctance/conf.toml
+++ b/examples/ActivityDrivenReluctance/conf.toml
@@ -24,8 +24,9 @@ mean_weights = false    # Use the meanfield approximation of the network edges
 
 reluctances = true       # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1
 reluctance_mean = 1.0    # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0)
-reluctance_sigma = 0.25  # Width of normal distribution (before truncating)
+reluctance_sigma = 0.15  # Width of normal distribution (before truncating)
 reluctance_eps = 0.01    # Minimum such that the normal distribution is truncated at this value
+covariance_factor = 0.0  # 0.0 means that the reluctances and activities are uncorrelated. Should be in the range of [-1,1]
 
 [network]
 number_of_agents = 1000
diff --git a/examples/DeffuantVector/conf.toml b/examples/DeffuantVector/conf.toml
index 604a40e..1d5f0ca 100644
--- a/examples/DeffuantVector/conf.toml
+++ b/examples/DeffuantVector/conf.toml
@@ -5,7 +5,7 @@ model = "Deffuant"
 [io]
 # n_output_network = 20 # Write the network every 20 iterations
 n_output_agents = 1 # Write the opinions of agents after every iteration
-print_progress = true # Print the iteration time ; if not set, then does not prints
+print_progress = true # Print the iteration time ; if not set, then does not print
 output_initial = true # Print the initial opinions and network file from step 0. If not set, this is true by default.
 start_output = 1 # Start writing out opinions and/or network files from this iteration. If not set, this is 1.
 
diff --git a/include/agents/inertial_agent.hpp b/include/agents/inertial_agent.hpp
new file mode 100644
index 0000000..469376a
--- /dev/null
+++ b/include/agents/inertial_agent.hpp
@@ -0,0 +1,68 @@
+#pragma once
+
+#include "agent.hpp"
+#include "agent_io.hpp"
+#include <util/misc.hpp>
+
+namespace Seldon
+{
+
+struct InertialAgentData
+{
+    double opinion    = 0;   // x_i
+    double activity   = 0;   // alpha_i
+    double reluctance = 1.0; // m_i
+    double velocity   = 0.0; // d(x_i)/dt
+};
+
+using InertialAgent = Agent<InertialAgentData>;
+
+template<>
+inline std::string agent_to_string<InertialAgent>( const InertialAgent & agent )
+{
+    return fmt::format(
+        "{}, {}, {}, {}", agent.data.opinion, agent.data.velocity, agent.data.activity, agent.data.reluctance );
+}
+
+template<>
+inline std::string opinion_to_string<InertialAgent>( const InertialAgent & agent )
+{
+    return fmt::format( "{}", agent.data.opinion );
+}
+
+template<>
+inline InertialAgent agent_from_string<InertialAgent>( const std::string & str )
+{
+    InertialAgent res{};
+
+    auto callback = [&]( int idx_list, std::string & substr )
+    {
+        if( idx_list == 0 )
+        {
+            res.data.opinion = std::stod( substr );
+        }
+        else if( idx_list == 1 )
+        {
+            res.data.velocity = std::stod( substr );
+        }
+        else if( idx_list == 2 )
+        {
+            res.data.activity = std::stod( substr );
+        }
+        else if( idx_list == 3 )
+        {
+            res.data.reluctance = std::stod( substr );
+        }
+    };
+
+    Seldon::parse_comma_separated_list( str, callback );
+
+    return res;
+};
+
+template<>
+inline std::vector<std::string> agent_to_string_column_names<InertialAgent>()
+{
+    return { "opinion", "velocity", "activity", "reluctance" };
+}
+} // namespace Seldon
\ No newline at end of file
diff --git a/include/config_parser.hpp b/include/config_parser.hpp
index 208e27d..326dc59 100644
--- a/include/config_parser.hpp
+++ b/include/config_parser.hpp
@@ -24,7 +24,8 @@ namespace Seldon::Config
 enum class Model
 {
     DeGroot,
-    ActivityDrivenModel,
+    ActivityDrivenModel, // @TODO : no need for model here
+    ActivityDrivenInertial,
     DeffuantModel
 };
 
@@ -81,6 +82,11 @@ struct ActivityDrivenSettings
     double covariance_factor          = 0.0;
 };
 
+struct ActivityDrivenInertialSettings : public ActivityDrivenSettings
+{
+    double friction_coefficient = 1.0;
+};
+
 struct InitialNetworkSettings
 {
     std::optional<std::string> file;
@@ -90,7 +96,8 @@ struct InitialNetworkSettings
 
 struct SimulationOptions
 {
-    using ModelVariantT = std::variant<DeGrootSettings, ActivityDrivenSettings, DeffuantSettings>;
+    using ModelVariantT
+        = std::variant<DeGrootSettings, ActivityDrivenSettings, ActivityDrivenInertialSettings, DeffuantSettings>;
     Model model;
     std::string model_string;
     int rng_seed = std::random_device()();
diff --git a/include/model_factory.hpp b/include/model_factory.hpp
index d5864ac..a87f9dd 100644
--- a/include/model_factory.hpp
+++ b/include/model_factory.hpp
@@ -3,6 +3,7 @@
 #include "models/ActivityDrivenModel.hpp"
 #include "models/DeGroot.hpp"
 #include "models/DeffuantModel.hpp"
+#include "models/InertialModel.hpp"
 #include "network.hpp"
 #include <memory>
 #include <random>
@@ -61,6 +62,23 @@ create_model_activity_driven( Network<AgentT> & network, const ModelVariantT & m
     }
 }
 
+template<typename AgentT>
+inline auto create_model_activity_driven_inertial(
+    Network<AgentT> & network, const ModelVariantT & model_settings, std::mt19937 & gen )
+{
+    if constexpr( std::is_same_v<AgentT, InertialModel::AgentT> )
+    {
+        auto settings = std::get<Config::ActivityDrivenInertialSettings>( model_settings );
+        auto model    = std::make_unique<InertialModel>( settings, network, gen );
+        return model;
+    }
+    else
+    {
+        throw std::runtime_error( "Incompatible agent and model type!" );
+        return std::unique_ptr<Model<AgentT>>{};
+    }
+}
+
 template<typename AgentT>
 inline auto create_model_deffuant( Network<AgentT> & network, const ModelVariantT & model_settings, std::mt19937 & gen )
 {
diff --git a/include/models/ActivityDrivenModel.hpp b/include/models/ActivityDrivenModel.hpp
index 06d830f..64a1223 100644
--- a/include/models/ActivityDrivenModel.hpp
+++ b/include/models/ActivityDrivenModel.hpp
@@ -1,10 +1,11 @@
 #pragma once
 
 #include "agents/activity_agent.hpp"
+#include "agents/inertial_agent.hpp"
 #include "config_parser.hpp"
-
 #include "model.hpp"
 #include "network.hpp"
+#include "network_generation.hpp"
 #include <cstddef>
 #include <random>
 #include <set>
@@ -15,23 +16,63 @@
 namespace Seldon
 {
 
-class ActivityDrivenModel : public Model<ActivityAgent>
+template<typename AgentT_>
+class ActivityDrivenModelAbstract : public Model<AgentT_>
 {
 public:
-    using AgentT   = ActivityAgent;
+    using AgentT   = AgentT_;
     using NetworkT = Network<AgentT>;
+    using WeightT  = typename NetworkT::WeightT;
 
-    ActivityDrivenModel( const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen );
+    ActivityDrivenModelAbstract(
+        const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen )
+            : Model<AgentT>( settings.max_iterations ),
+              network( network ),
+              contact_prob_list( std::vector<std::vector<WeightT>>( network.n_agents() ) ),
+              gen( gen ),
+              dt( settings.dt ),
+              m( settings.m ),
+              eps( settings.eps ),
+              gamma( settings.gamma ),
+              alpha( settings.alpha ),
+              homophily( settings.homophily ),
+              reciprocity( settings.reciprocity ),
+              K( settings.K ),
+              mean_activities( settings.mean_activities ),
+              mean_weights( settings.mean_weights ),
+              use_reluctances( settings.use_reluctances ),
+              reluctance_mean( settings.reluctance_mean ),
+              reluctance_sigma( settings.reluctance_sigma ),
+              reluctance_eps( settings.reluctance_eps ),
+              covariance_factor( settings.covariance_factor ),
+              n_bots( settings.n_bots ),
+              bot_m( settings.bot_m ),
+              bot_activity( settings.bot_activity ),
+              bot_opinion( settings.bot_opinion ),
+              bot_homophily( settings.bot_homophily )
+    {
+        get_agents_from_power_law();
 
-    void iteration() override;
+        if( mean_weights )
+        {
+            auto agents_copy = network.agents;
+            network          = NetworkGeneration::generate_fully_connected<AgentT>( network.n_agents() );
+            network.agents   = agents_copy;
+        }
+    }
 
-private:
+    void iteration() override {};
+
+protected:
     NetworkT & network;
-    std::vector<std::vector<NetworkT::WeightT>> contact_prob_list; // Probability of choosing i in 1 to m rounds
+
+private:
+    std::vector<std::vector<WeightT>> contact_prob_list; // Probability of choosing i in 1 to m rounds
     // Random number generation
     std::mt19937 & gen; // reference to simulation Mersenne-Twister engine
     std::set<std::pair<size_t, size_t>> reciprocal_edge_buffer{};
 
+protected:
     // Model-specific parameters
     double dt{}; // Timestep for the integration of the coupled ODEs
     // Various free parameters
@@ -66,13 +107,222 @@ class ActivityDrivenModel : public Model<ActivityAgent>
     std::vector<double> k3_buffer{};
     std::vector<double> k4_buffer{};
 
-    void get_agents_from_power_law();
+private:
+    void get_agents_from_power_law()
+    {
+        std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values
+        power_law_distribution<> dist_activity( eps, gamma );
+        truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps );
+
+        bivariate_gaussian_copula copula( covariance_factor, dist_activity, dist_reluctance );
+
+        auto mean_activity = dist_activity.mean();
 
+        // Initial conditions for the opinions, initialize to [-1,1]
+        // The activities should be drawn from a power law distribution
+        for( size_t i = 0; i < network.agents.size(); i++ )
+        {
+            network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value
+
+            auto res                        = copula( gen );
+            network.agents[i].data.activity = res[0];
+
+            if( use_reluctances )
+            {
+                network.agents[i].data.reluctance = res[1];
+            }
+            if( mean_activities )
+            {
+                network.agents[i].data.activity = mean_activity;
+            }
+        }
+
+        if( bot_present() )
+        {
+            for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
+            {
+                network.agents[bot_idx].data.opinion  = bot_opinion[bot_idx];
+                network.agents[bot_idx].data.activity = bot_activity[bot_idx];
+            }
+        }
+    }
+
+    // The weight for contact between two agents
+    double homophily_weight( size_t idx_contacter, size_t idx_contacted )
+    {
+        double homophily = this->homophily;
+
+        if( idx_contacted == idx_contacter )
+            return 0.0;
+
+        if( bot_present() && idx_contacter < n_bots )
+            homophily = this->bot_homophily[idx_contacter];
+
+        constexpr double tolerance = 1e-10;
+        auto opinion_diff
+            = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion );
+        opinion_diff = std::max( tolerance, opinion_diff );
+
+        return std::pow( opinion_diff, -homophily );
+    }
+
+    void update_network_probabilistic()
+    {
+        network.switch_direction_flag();
+
+        std::uniform_real_distribution<> dis_activation( 0.0, 1.0 );
+        std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 );
+        std::vector<size_t> contacted_agents{};
+        reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Test if the agent is activated
+            bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity;
+
+            if( activated )
+            {
+                // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
+                // Not normalised since this is taken care of by the reservoir sampling
+
+                int m_temp = this->m;
+
+                if( bot_present() && idx_agent < n_bots )
+                {
+                    m_temp = bot_m[idx_agent];
+                }
+
+                reservoir_sampling_A_ExpJ(
+                    m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); },
+                    contacted_agents, gen );
+
+                // Fill the outgoing edges into the reciprocal edge buffer
+                for( const auto & idx_outgoing : contacted_agents )
+                {
+                    reciprocal_edge_buffer.insert(
+                        { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing
+                }
+
+                // Set the *outgoing* edges
+                network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 );
+            }
+            else
+            {
+                network.set_neighbours_and_weights( idx_agent, {}, {} );
+            }
+        }
+
+        // Reciprocity check
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Get the outgoing edges
+            auto contacted_agents = network.get_neighbours( idx_agent );
+            // For each outgoing edge we check if the reverse edge already exists
+            for( const auto & idx_outgoing : contacted_agents )
+            {
+                // If the edge is not reciprocated
+                if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) )
+                {
+                    if( dis_reciprocation( gen ) < reciprocity )
+                    {
+                        network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 );
+                    }
+                }
+            }
+        }
+
+        network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges
+    }
+
+    void update_network_mean()
+    {
+        std::vector<WeightT> weights( network.n_agents(), 0.0 );
+
+        // Set all weights to zero in the beginning
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            network.set_weights( idx_agent, weights );
+            contact_prob_list[idx_agent] = weights; // set to zero
+        }
+
+        auto probability_helper = []( double omega, size_t m )
+        {
+            double p = 0;
+            for( size_t i = 1; i <= m; i++ )
+                p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 );
+            return p;
+        };
+
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
+            // Not normalised since this is taken care of by the reservoir sampling
+
+            double normalization = 0;
+            for( size_t k = 0; k < network.n_agents(); k++ )
+            {
+                normalization += homophily_weight( idx_agent, k );
+            }
+
+            // Go through all the neighbours of idx_agent
+            // Calculate the probability of i contacting j (in 1 to m rounds, assuming
+            // the agent is activated
+            int m_temp = m;
+            if( bot_present() && idx_agent < n_bots )
+            {
+                m_temp = bot_m[idx_agent];
+            }
+
+            double activity = std::max( 1.0, network.agents[idx_agent].data.activity );
+            for( size_t j = 0; j < network.n_agents(); j++ )
+            {
+                double omega = homophily_weight( idx_agent, j ) / normalization;
+                // We can calculate the probability of i contacting j ( i->j )
+                // Update contact prob_list (outgoing)
+                contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp );
+            }
+        }
+
+        for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+        {
+            // Calculate the actual weights and reciprocity
+            for( size_t j = 0; j < network.n_agents(); j++ )
+            {
+                double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites
+                double prob_contact_ji = contact_prob_list[j][idx_agent];
+
+                // Set the incoming agent weight, j-i in weight list
+                double & win_ji = network.get_weights( j )[idx_agent];
+                win_ji += prob_contact_ij;
+
+                // Handle the reciprocity for j->i
+                // Update incoming weight i-j
+                double & win_ij = network.get_weights( idx_agent )[j];
+
+                // The probability of reciprocating is
+                win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij;
+            }
+        }
+    }
+
+protected:
     [[nodiscard]] bool bot_present() const
     {
         return n_bots > 0;
     }
 
+    void update_network()
+    {
+
+        if( !mean_weights )
+        {
+            update_network_probabilistic();
+        }
+        else
+        {
+            update_network_mean();
+        }
+    }
+
     template<typename Opinion_Callback>
     void get_euler_slopes( std::vector<double> & k_buffer, Opinion_Callback opinion )
     {
@@ -93,16 +343,12 @@ class ActivityDrivenModel : public Model<ActivityAgent>
                 k_buffer[idx_agent] += 1.0 / network.agents[idx_agent].data.reluctance * K * weight_buffer[j]
                                        * std::tanh( alpha * opinion( j_index ) );
             }
-            // Multiply by the timestep
-            k_buffer[idx_agent] *= dt;
+            // Here, we won't multiply by the timestep.
+            // Instead multiply in the update rule
         }
     }
-
-    // The weight for contact between two agents
-    double homophily_weight( size_t idx_contacter, size_t idx_contacted );
-    void update_network_probabilistic();
-    void update_network_mean();
-    void update_network();
 };
 
+using ActivityDrivenModel = ActivityDrivenModelAbstract<ActivityAgent>;
+
 } // namespace Seldon
\ No newline at end of file
diff --git a/include/models/InertialModel.hpp b/include/models/InertialModel.hpp
new file mode 100644
index 0000000..ba95a0b
--- /dev/null
+++ b/include/models/InertialModel.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+#include "agents/activity_agent.hpp"
+#include "agents/inertial_agent.hpp"
+#include "config_parser.hpp"
+#include "model.hpp"
+#include "models/ActivityDrivenModel.hpp"
+#include "network.hpp"
+#include "network_generation.hpp"
+#include <cstddef>
+#include <random>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Seldon
+{
+
+class InertialModel : public ActivityDrivenModelAbstract<InertialAgent>
+{
+public:
+    using AgentT   = InertialAgent;
+    using NetworkT = Network<AgentT>;
+    using WeightT  = typename NetworkT::WeightT;
+
+    InertialModel( const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen );
+
+    void iteration() override;
+
+private:
+    double friction_coefficient = 1.0;
+    std::vector<double> drift_t_buffer{};
+    std::vector<double> drift_next_t_buffer{};
+
+    void calc_velocity();
+    void calc_position();
+};
+
+} // namespace Seldon
\ No newline at end of file
diff --git a/include/simulation.hpp b/include/simulation.hpp
index 4eca2fc..118b824 100644
--- a/include/simulation.hpp
+++ b/include/simulation.hpp
@@ -70,6 +70,10 @@ class Simulation : public SimulationInterface
         {
             model = ModelFactory::create_model_activity_driven( network, options.model_settings, gen );
         }
+        else if( options.model == Config::Model::ActivityDrivenInertial )
+        {
+            model = ModelFactory::create_model_activity_driven_inertial( network, options.model_settings, gen );
+        }
         else if( options.model == Config::Model::DeffuantModel )
         {
             auto deffuant_settings = std::get<Config::DeffuantSettings>( options.model_settings );
diff --git a/include/util/erfinv.hpp b/include/util/erfinv.hpp
new file mode 100644
index 0000000..6096e54
--- /dev/null
+++ b/include/util/erfinv.hpp
@@ -0,0 +1,112 @@
+#pragma once
+#include <cmath>
+#include <limits>
+
+namespace Seldon::Math
+{
+
+constexpr long double pi      = 3.1415926535897932384626433832795028841971693993751L;
+constexpr long double sqrt_pi = 1.7724538509055160272981674833411451827975494561224L;
+
+// Implementation adapted from https://github.com/lakshayg/erfinv same as used in golang math library
+template<typename T>
+T erfinv( T x )
+{
+    if( x < -1 || x > 1 )
+    {
+        return std::numeric_limits<T>::quiet_NaN();
+    }
+    else if( x == 1.0 )
+    {
+        return std::numeric_limits<T>::infinity();
+    }
+    else if( x == -1.0 )
+    {
+        return -std::numeric_limits<T>::infinity();
+    }
+
+    const T LN2 = 6.931471805599453094172321214581e-1L;
+
+    const T A0 = 1.1975323115670912564578e0L;
+    const T A1 = 4.7072688112383978012285e1L;
+    const T A2 = 6.9706266534389598238465e2L;
+    const T A3 = 4.8548868893843886794648e3L;
+    const T A4 = 1.6235862515167575384252e4L;
+    const T A5 = 2.3782041382114385731252e4L;
+    const T A6 = 1.1819493347062294404278e4L;
+    const T A7 = 8.8709406962545514830200e2L;
+
+    const T B0 = 1.0000000000000000000e0L;
+    const T B1 = 4.2313330701600911252e1L;
+    const T B2 = 6.8718700749205790830e2L;
+    const T B3 = 5.3941960214247511077e3L;
+    const T B4 = 2.1213794301586595867e4L;
+    const T B5 = 3.9307895800092710610e4L;
+    const T B6 = 2.8729085735721942674e4L;
+    const T B7 = 5.2264952788528545610e3L;
+
+    const T C0 = 1.42343711074968357734e0L;
+    const T C1 = 4.63033784615654529590e0L;
+    const T C2 = 5.76949722146069140550e0L;
+    const T C3 = 3.64784832476320460504e0L;
+    const T C4 = 1.27045825245236838258e0L;
+    const T C5 = 2.41780725177450611770e-1L;
+    const T C6 = 2.27238449892691845833e-2L;
+    const T C7 = 7.74545014278341407640e-4L;
+
+    const T D0 = 1.4142135623730950488016887e0L;
+    const T D1 = 2.9036514445419946173133295e0L;
+    const T D2 = 2.3707661626024532365971225e0L;
+    const T D3 = 9.7547832001787427186894837e-1L;
+    const T D4 = 2.0945065210512749128288442e-1L;
+    const T D5 = 2.1494160384252876777097297e-2L;
+    const T D6 = 7.7441459065157709165577218e-4L;
+    const T D7 = 1.4859850019840355905497876e-9L;
+
+    const T E0 = 6.65790464350110377720e0L;
+    const T E1 = 5.46378491116411436990e0L;
+    const T E2 = 1.78482653991729133580e0L;
+    const T E3 = 2.96560571828504891230e-1L;
+    const T E4 = 2.65321895265761230930e-2L;
+    const T E5 = 1.24266094738807843860e-3L;
+    const T E6 = 2.71155556874348757815e-5L;
+    const T E7 = 2.01033439929228813265e-7L;
+
+    const T F0 = 1.414213562373095048801689e0L;
+    const T F1 = 8.482908416595164588112026e-1L;
+    const T F2 = 1.936480946950659106176712e-1L;
+    const T F3 = 2.103693768272068968719679e-2L;
+    const T F4 = 1.112800997078859844711555e-3L;
+    const T F5 = 2.611088405080593625138020e-5L;
+    const T F6 = 2.010321207683943062279931e-7L;
+    const T F7 = 2.891024605872965461538222e-15L;
+
+    T abs_x = std::abs( x );
+
+    T r, num, den;
+
+    if( abs_x <= 0.85 )
+    {
+        r   = 0.180625 - 0.25 * x * x;
+        num = ( ( ( ( ( ( ( A7 * r + A6 ) * r + A5 ) * r + A4 ) * r + A3 ) * r + A2 ) * r + A1 ) * r + A0 );
+        den = ( ( ( ( ( ( ( B7 * r + B6 ) * r + B5 ) * r + B4 ) * r + B3 ) * r + B2 ) * r + B1 ) * r + B0 );
+        return x * num / den;
+    }
+
+    r = std::sqrt( LN2 - std::log1p( -abs_x ) );
+    if( r <= 5.0 )
+    {
+        r   = r - 1.6;
+        num = ( ( ( ( ( ( ( C7 * r + C6 ) * r + C5 ) * r + C4 ) * r + C3 ) * r + C2 ) * r + C1 ) * r + C0 );
+        den = ( ( ( ( ( ( ( D7 * r + D6 ) * r + D5 ) * r + D4 ) * r + D3 ) * r + D2 ) * r + D1 ) * r + D0 );
+    }
+    else
+    {
+        r   = r - 5.0;
+        num = ( ( ( ( ( ( ( E7 * r + E6 ) * r + E5 ) * r + E4 ) * r + E3 ) * r + E2 ) * r + E1 ) * r + E0 );
+        den = ( ( ( ( ( ( ( F7 * r + F6 ) * r + F5 ) * r + F4 ) * r + F3 ) * r + F2 ) * r + F1 ) * r + F0 );
+    }
+
+    return std::copysign<T>( num / den, x );
+}
+} // namespace Seldon::Math
\ No newline at end of file
diff --git a/include/util/math.hpp b/include/util/math.hpp
index e6041fc..493e6ab 100644
--- a/include/util/math.hpp
+++ b/include/util/math.hpp
@@ -1,4 +1,6 @@
 #pragma once
+#include "fmt/core.h"
+#include "util/erfinv.hpp"
 #include <algorithm>
 #include <cstddef>
 #include <optional>
@@ -129,9 +131,19 @@ class power_law_distribution
 
     template<typename Generator>
     ScalarT operator()( Generator & gen )
+    {
+        return inverse_cdf( dist( gen ) );
+    }
+
+    ScalarT pdf( ScalarT x )
+    {
+        return ( 1.0 - gamma ) / ( 1.0 - std::pow( eps, ( 1 - gamma ) ) * std::pow( x, ( -gamma ) ) );
+    }
+
+    ScalarT inverse_cdf( ScalarT x )
     {
         return std::pow(
-            ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * dist( gen ) + std::pow( eps, ( 1.0 - gamma ) ),
+            ( 1.0 - std::pow( eps, ( 1.0 - gamma ) ) ) * x + std::pow( eps, ( 1.0 - gamma ) ),
             ( 1.0 / ( 1.0 - gamma ) ) );
     }
 
@@ -154,25 +166,114 @@ class truncated_normal_distribution
     ScalarT mean{};
     ScalarT sigma{};
     ScalarT eps{};
-    std::normal_distribution<ScalarT> normal_dist{};
-    size_t max_tries = 5000;
+    std::uniform_real_distribution<ScalarT> uniform_dist{};
+
+    ScalarT inverse_cdf_gauss( ScalarT y )
+    {
+        return Math::erfinv( 2.0 * y - 1 ) * std::sqrt( 2.0 ) * sigma + mean;
+    }
+
+    ScalarT cdf_gauss( ScalarT x )
+    {
+        return 0.5 * ( 1 + std::erf( ( x - mean ) / ( sigma * std::sqrt( 2.0 ) ) ) );
+    }
+
+    ScalarT pdf_gauss( ScalarT x )
+    {
+        return 1.0 / ( sigma * std::sqrt( 2 * Math::pi ) ) * std::exp( -0.5 * std::pow( ( ( x - mean ) / sigma ), 2 ) );
+    }
 
 public:
     truncated_normal_distribution( ScalarT mean, ScalarT sigma, ScalarT eps )
-            : mean( mean ), sigma( sigma ), eps( eps ), normal_dist( std::normal_distribution<ScalarT>( mean, sigma ) )
+            : mean( mean ), sigma( sigma ), eps( eps ), uniform_dist( 0, 1 )
     {
     }
 
     template<typename Generator>
     ScalarT operator()( Generator & gen )
     {
-        for( size_t i = 0; i < max_tries; i++ )
-        {
-            auto sample = normal_dist( gen );
-            if( sample > eps )
-                return sample;
-        }
-        return eps;
+        return inverse_cdf( uniform_dist( gen ) );
+    }
+
+    ScalarT inverse_cdf( ScalarT y )
+    {
+        return inverse_cdf_gauss( y * ( 1.0 - cdf_gauss( eps ) ) + cdf_gauss( eps ) );
+    }
+
+    ScalarT pdf( ScalarT x )
+    {
+        if( x < eps )
+            return 0.0;
+        else
+            return 1.0 / ( 1.0 - cdf_gauss( eps ) ) * pdf_gauss( x );
+    }
+};
+
+/**
+ * @brief Bivariate normal distribution
+ * with mean mu = [0,0]
+ * and covariance matrix Sigma = [[1, cov], [cov, 1]]
+ * |cov| < 1 is required
+ */
+template<typename ScalarT = double>
+class bivariate_normal_distribution
+{
+private:
+    ScalarT covariance;
+    std::normal_distribution<ScalarT> normal_dist{};
+
+public:
+    bivariate_normal_distribution( ScalarT covariance ) : covariance( covariance ) {}
+
+    template<typename Generator>
+    std::array<ScalarT, 2> operator()( Generator & gen )
+    {
+        ScalarT n1 = normal_dist( gen );
+        ScalarT n2 = normal_dist( gen );
+
+        ScalarT r1 = n1;
+        ScalarT r2 = covariance * n1 + std::sqrt( 1.0 - covariance * covariance ) * n2;
+
+        return { r1, r2 };
+    }
+};
+
+template<typename ScalarT, typename dist1T, typename dist2T>
+class bivariate_gaussian_copula
+{
+private:
+    ScalarT covariance;
+    bivariate_normal_distribution<ScalarT> biv_normal_dist{};
+    // std::normal_distribution<ScalarT> normal_dist{};
+
+    // Cumulative probability function for gaussian with mean 0 and variance 1
+    ScalarT cdf_gauss( ScalarT x )
+    {
+        return 0.5 * ( 1 + std::erf( ( x ) / std::sqrt( 2.0 ) ) );
+    }
+
+    dist1T dist1;
+    dist2T dist2;
+
+public:
+    bivariate_gaussian_copula( ScalarT covariance, dist1T dist1, dist2T dist2 )
+            : covariance( covariance ),
+              biv_normal_dist( bivariate_normal_distribution( covariance ) ),
+              dist1( dist1 ),
+              dist2( dist2 )
+    {
+    }
+
+    template<typename Generator>
+    std::array<ScalarT, 2> operator()( Generator & gen )
+    {
+        // 1. Draw from bivariate gaussian
+        auto z = biv_normal_dist( gen );
+        // 2. Transform marginals to unit interval
+        std::array<ScalarT, 2> z_unit = { cdf_gauss( z[0] ), cdf_gauss( z[1] ) };
+        // 3. Apply inverse transform sampling
+        std::array<ScalarT, 2> res = { dist1.inverse_cdf( z_unit[0] ), dist2.inverse_cdf( z_unit[1] ) };
+        return res;
     }
 };
 
diff --git a/meson.build b/meson.build
index e7017c3..c9fd40d 100644
--- a/meson.build
+++ b/meson.build
@@ -12,6 +12,7 @@ sources_seldon = [
   'src/models/ActivityDrivenModel.cpp',
   'src/models/DeffuantModel.cpp',
   'src/models/DeffuantModelVector.cpp',
+  'src/models/InertialModel.cpp',
   'src/util/tomlplusplus.cpp',
 ]
 
@@ -26,12 +27,14 @@ tests = [
   ['Test_Tarjan', 'test/test_tarjan.cpp'],
   ['Test_DeGroot', 'test/test_deGroot.cpp'],
   ['Test_Activity_Driven', 'test/test_activity.cpp'],
+  ['Test_Activity_Driven_Inertial', 'test/test_activity_inertial.cpp'],
   ['Test_Deffuant', 'test/test_deffuant.cpp'],
   ['Test_Network', 'test/test_network.cpp'],
   ['Test_Network_Generation', 'test/test_network_generation.cpp'],
   ['Test_Sampling', 'test/test_sampling.cpp'],
   ['Test_IO', 'test/test_io.cpp'],
   ['Test_Util', 'test/test_util.cpp'],
+  ['Test_Prob', 'test/test_probability_distributions.cpp'],
 ]
 
 Catch2 = dependency('Catch2', method : 'cmake', modules : ['Catch2::Catch2WithMain', 'Catch2::Catch2'])
diff --git a/src/config_parser.cpp b/src/config_parser.cpp
index b65d83a..a492e49 100644
--- a/src/config_parser.cpp
+++ b/src/config_parser.cpp
@@ -29,9 +29,83 @@ Model model_string_to_enum( std::string_view model_string )
     {
         return Model::DeffuantModel;
     }
+    else if( model_string == "ActivityDrivenInertial" )
+    {
+        return Model::ActivityDrivenInertial;
+    }
     throw std::runtime_error( fmt::format( "Invalid model string {}", model_string ) );
 }
 
+void set_if_specified( auto & opt, const auto & toml_opt )
+{
+    using T    = typename std::remove_reference<decltype( opt )>::type;
+    auto t_opt = toml_opt.template value<T>();
+    if( t_opt.has_value() )
+        opt = t_opt.value();
+}
+
+// Convenience function to parse activity settings
+void parse_activity_settings( auto & model_settings, const auto & toml_model_opt, const auto & tbl )
+{
+    set_if_specified( model_settings.dt, toml_model_opt["dt"] );
+    set_if_specified( model_settings.m, toml_model_opt["m"] );
+    set_if_specified( model_settings.eps, toml_model_opt["eps"] );
+    set_if_specified( model_settings.gamma, toml_model_opt["gamma"] );
+    set_if_specified( model_settings.homophily, toml_model_opt["homophily"] );
+    set_if_specified( model_settings.reciprocity, toml_model_opt["reciprocity"] );
+    set_if_specified( model_settings.alpha, toml_model_opt["alpha"] );
+    set_if_specified( model_settings.K, toml_model_opt["K"] );
+    // Mean activity model options
+    set_if_specified( model_settings.mean_activities, toml_model_opt["mean_activities"] );
+    set_if_specified( model_settings.mean_weights, toml_model_opt["mean_weights"] );
+    // Reluctances
+    set_if_specified( model_settings.covariance_factor, toml_model_opt["covariance_factor"] );
+    set_if_specified( model_settings.use_reluctances, toml_model_opt["reluctances"] );
+    set_if_specified( model_settings.reluctance_mean, toml_model_opt["reluctance_mean"] );
+    set_if_specified( model_settings.reluctance_sigma, toml_model_opt["reluctance_sigma"] );
+    set_if_specified( model_settings.reluctance_eps, toml_model_opt["reluctance_eps"] );
+
+    model_settings.max_iterations = tbl["model"]["max_iterations"].template value<int>();
+
+    // bot
+    set_if_specified( model_settings.n_bots, toml_model_opt["n_bots"] );
+
+    auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value )
+    {
+        if( toml_node.is_array() )
+        {
+            toml::array * toml_arr = toml_node.as_array();
+
+            toml_arr->for_each(
+                [&]( auto && elem )
+                {
+                    if( elem.is_integer() )
+                    {
+                        options_array.push_back( elem.as_integer()->get() );
+                    }
+                    else if( elem.is_floating_point() )
+                    {
+                        options_array.push_back( elem.as_floating_point()->get() );
+                    }
+                } );
+        }
+        else
+        {
+            options_array = std::vector<decltype( default_value )>( model_settings.n_bots, default_value );
+        }
+    };
+
+    auto bot_opinion   = toml_model_opt["bot_opinion"];
+    auto bot_m         = toml_model_opt["bot_m"];
+    auto bot_activity  = toml_model_opt["bot_activity"];
+    auto bot_homophily = toml_model_opt["bot_homophily"];
+
+    push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m );
+    push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 );
+    push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 );
+    push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily );
+}
+
 SimulationOptions parse_config_file( std::string_view config_file_path )
 {
     SimulationOptions options;
@@ -41,14 +115,6 @@ SimulationOptions parse_config_file( std::string_view config_file_path )
 
     options.rng_seed = tbl["simulation"]["rng_seed"].value_or( int( options.rng_seed ) );
 
-    auto set_if_specified = [&]( auto & opt, const auto & toml_opt )
-    {
-        using T    = typename std::remove_reference<decltype( opt )>::type;
-        auto t_opt = toml_opt.template value<T>();
-        if( t_opt.has_value() )
-            opt = t_opt.value();
-    };
-
     // Parse output settings
     options.output_settings.n_output_network = tbl["io"]["n_output_network"].value<size_t>();
     options.output_settings.n_output_agents  = tbl["io"]["n_output_agents"].value<size_t>();
@@ -87,63 +153,15 @@ SimulationOptions parse_config_file( std::string_view config_file_path )
     {
         auto model_settings = ActivityDrivenSettings();
 
-        set_if_specified( model_settings.dt, tbl[options.model_string]["dt"] );
-        set_if_specified( model_settings.m, tbl[options.model_string]["m"] );
-        set_if_specified( model_settings.eps, tbl[options.model_string]["eps"] );
-        set_if_specified( model_settings.gamma, tbl[options.model_string]["gamma"] );
-        set_if_specified( model_settings.homophily, tbl[options.model_string]["homophily"] );
-        set_if_specified( model_settings.reciprocity, tbl[options.model_string]["reciprocity"] );
-        set_if_specified( model_settings.alpha, tbl[options.model_string]["alpha"] );
-        set_if_specified( model_settings.K, tbl[options.model_string]["K"] );
-        // Mean activity model options
-        set_if_specified( model_settings.mean_activities, tbl[options.model_string]["mean_activities"] );
-        set_if_specified( model_settings.mean_weights, tbl[options.model_string]["mean_weights"] );
-        // Reluctances
-        set_if_specified( model_settings.use_reluctances, tbl[options.model_string]["reluctances"] );
-        set_if_specified( model_settings.reluctance_mean, tbl[options.model_string]["reluctance_mean"] );
-        set_if_specified( model_settings.reluctance_sigma, tbl[options.model_string]["reluctance_sigma"] );
-        set_if_specified( model_settings.reluctance_eps, tbl[options.model_string]["reluctance_eps"] );
-
-        model_settings.max_iterations = tbl["model"]["max_iterations"].value<int>();
-
-        // bot
-        set_if_specified( model_settings.n_bots, tbl[options.model_string]["n_bots"] );
-
-        auto push_back_bot_array = [&]( auto toml_node, auto & options_array, auto default_value )
-        {
-            if( toml_node.is_array() )
-            {
-                toml::array * toml_arr = toml_node.as_array();
-
-                toml_arr->for_each(
-                    [&]( auto && elem )
-                    {
-                        if( elem.is_integer() )
-                        {
-                            options_array.push_back( elem.as_integer()->get() );
-                        }
-                        else if( elem.is_floating_point() )
-                        {
-                            options_array.push_back( elem.as_floating_point()->get() );
-                        }
-                    } );
-            }
-            else
-            {
-                options_array = std::vector<decltype( default_value )>( model_settings.n_bots, default_value );
-            }
-        };
-
-        auto bot_opinion   = tbl[options.model_string]["bot_opinion"];
-        auto bot_m         = tbl[options.model_string]["bot_m"];
-        auto bot_activity  = tbl[options.model_string]["bot_activity"];
-        auto bot_homophily = tbl[options.model_string]["bot_homophily"];
-
-        push_back_bot_array( bot_m, model_settings.bot_m, model_settings.m );
-        push_back_bot_array( bot_opinion, model_settings.bot_opinion, 0.0 );
-        push_back_bot_array( bot_activity, model_settings.bot_activity, 1.0 );
-        push_back_bot_array( bot_homophily, model_settings.bot_homophily, model_settings.homophily );
+        parse_activity_settings( model_settings, tbl[options.model_string], tbl );
+        options.model_settings = model_settings;
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = ActivityDrivenInertialSettings();
 
+        parse_activity_settings( model_settings, tbl[options.model_string], tbl );
+        set_if_specified( model_settings.friction_coefficient, tbl[options.model_string]["friction_coefficient"] );
         options.model_settings = model_settings;
     }
 
@@ -183,10 +201,8 @@ void validate_settings( const SimulationOptions & options )
     // @TODO: Check that start_output is less than the max_iterations?
     check( name_and_var( options.output_settings.start_output ), g_zero );
 
-    if( options.model == Model::ActivityDrivenModel )
+    auto validate_activity = [&]( const auto & model_settings )
     {
-        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
-
         check( name_and_var( model_settings.dt ), g_zero );
         check( name_and_var( model_settings.m ), geq_zero );
         check( name_and_var( model_settings.eps ), g_zero );
@@ -198,7 +214,7 @@ void validate_settings( const SimulationOptions & options )
         check( name_and_var( model_settings.reluctance_mean ), g_zero );
         check( name_and_var( model_settings.reluctance_sigma ), g_zero );
         check( name_and_var( model_settings.reluctance_eps ), g_zero );
-        check( name_and_var( model_settings.covariance_factor ), geq_zero );
+        check( name_and_var( model_settings.covariance_factor ), []( auto x ) { return x >= -1.0 && x <= 1.0; } );
         // Bot options
         size_t n_bots             = model_settings.n_bots;
         auto check_bot_size       = [&]( auto x ) { return x.size() >= n_bots; };
@@ -207,6 +223,19 @@ void validate_settings( const SimulationOptions & options )
         check( name_and_var( model_settings.bot_activity ), check_bot_size, bot_msg );
         check( name_and_var( model_settings.bot_opinion ), check_bot_size, bot_msg );
         check( name_and_var( model_settings.bot_homophily ), check_bot_size, bot_msg );
+    };
+
+    if( options.model == Model::ActivityDrivenModel )
+    {
+        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
+
+        validate_activity( model_settings );
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = std::get<ActivityDrivenInertialSettings>( options.model_settings );
+        check( name_and_var( model_settings.friction_coefficient ), geq_zero );
+        validate_activity( model_settings );
     }
     else if( options.model == Model::DeGroot )
     {
@@ -237,11 +266,8 @@ void print_settings( const SimulationOptions & options )
     fmt::print( "[Model]\n" );
     fmt::print( "    type {}\n", options.model_string );
 
-    // @TODO: Optionally print *all* settings to the console, including defaults that were set
-    if( options.model == Model::ActivityDrivenModel )
+    auto print_activity_settings = [&]( auto model_settings )
     {
-        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
-
         fmt::print( "    max_iterations {}\n", model_settings.max_iterations );
         fmt::print( "    dt {} \n", model_settings.dt );
         fmt::print( "    m {} \n", model_settings.m );
@@ -270,6 +296,20 @@ void print_settings( const SimulationOptions & options )
             fmt::print( "    reluctance_eps {}\n", model_settings.reluctance_eps );
             fmt::print( "    covariance_factor {}\n", model_settings.covariance_factor );
         }
+    };
+
+    // @TODO: Optionally print *all* settings to the console, including defaults that were set
+    if( options.model == Model::ActivityDrivenModel )
+    {
+        auto model_settings = std::get<ActivityDrivenSettings>( options.model_settings );
+
+        print_activity_settings( model_settings );
+    }
+    else if( options.model == Model::ActivityDrivenInertial )
+    {
+        auto model_settings = std::get<ActivityDrivenInertialSettings>( options.model_settings );
+        print_activity_settings( model_settings );
+        fmt::print( "    friction_coefficient {}\n", model_settings.friction_coefficient );
     }
     else if( options.model == Model::DeGroot )
     {
diff --git a/src/main.cpp b/src/main.cpp
index 9343730..ac57310 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,7 @@
 #include "config_parser.hpp"
 #include "models/DeGroot.hpp"
 #include "models/DeffuantModel.hpp"
+#include "models/InertialModel.hpp"
 #include "simulation.hpp"
 #include <fmt/format.h>
 #include <fmt/ostream.h>
@@ -70,11 +71,25 @@ int main( int argc, char * argv[] )
         simulation = std::make_unique<Seldon::Simulation<Seldon::ActivityDrivenModel::AgentT>>(
             simulation_options, network_file, agent_file );
     }
-    else if( simulation_options.model == Seldon::Config::Model::DeffuantModel )
+    else if( simulation_options.model == Seldon::Config::Model::ActivityDrivenInertial )
     {
-        simulation = std::make_unique<Seldon::Simulation<Seldon::DeffuantModel::AgentT>>(
+        simulation = std::make_unique<Seldon::Simulation<Seldon::InertialModel::AgentT>>(
             simulation_options, network_file, agent_file );
     }
+    else if( simulation_options.model == Seldon::Config::Model::DeffuantModel )
+    {
+        auto model_settings = std::get<Seldon::Config::DeffuantSettings>( simulation_options.model_settings );
+        if( model_settings.use_binary_vector )
+        {
+            simulation = std::make_unique<Seldon::Simulation<Seldon::DeffuantModelVector::AgentT>>(
+                simulation_options, network_file, agent_file );
+        }
+        else
+        {
+            simulation = std::make_unique<Seldon::Simulation<Seldon::DeffuantModel::AgentT>>(
+                simulation_options, network_file, agent_file );
+        }
+    }
     else
     {
         throw std::runtime_error( "Model has not been created" );
diff --git a/src/models/ActivityDrivenModel.cpp b/src/models/ActivityDrivenModel.cpp
index 0a0ba06..0b0dd9c 100644
--- a/src/models/ActivityDrivenModel.cpp
+++ b/src/models/ActivityDrivenModel.cpp
@@ -1,6 +1,5 @@
 #include "models/ActivityDrivenModel.hpp"
 #include "network.hpp"
-#include "network_generation.hpp"
 #include "util/math.hpp"
 #include <cstddef>
 #include <random>
@@ -9,252 +8,8 @@
 namespace Seldon
 {
 
-ActivityDrivenModel::ActivityDrivenModel(
-    const Config::ActivityDrivenSettings & settings, NetworkT & network, std::mt19937 & gen )
-        : Model<ActivityDrivenModel::AgentT>( settings.max_iterations ),
-          network( network ),
-          contact_prob_list( std::vector<std::vector<NetworkT::WeightT>>( network.n_agents() ) ),
-          gen( gen ),
-          dt( settings.dt ),
-          m( settings.m ),
-          eps( settings.eps ),
-          gamma( settings.gamma ),
-          alpha( settings.alpha ),
-          homophily( settings.homophily ),
-          reciprocity( settings.reciprocity ),
-          K( settings.K ),
-          mean_activities( settings.mean_activities ),
-          mean_weights( settings.mean_weights ),
-          use_reluctances( settings.use_reluctances ),
-          reluctance_mean( settings.reluctance_mean ),
-          reluctance_sigma( settings.reluctance_sigma ),
-          reluctance_eps( settings.reluctance_eps ),
-          n_bots( settings.n_bots ),
-          bot_m( settings.bot_m ),
-          bot_activity( settings.bot_activity ),
-          bot_opinion( settings.bot_opinion ),
-          bot_homophily( settings.bot_homophily )
-{
-    get_agents_from_power_law();
-
-    if( mean_weights )
-    {
-        auto agents_copy = network.agents;
-        network          = NetworkGeneration::generate_fully_connected<AgentT>( network.n_agents() );
-        network.agents   = agents_copy;
-    }
-}
-
-double ActivityDrivenModel::homophily_weight( size_t idx_contacter, size_t idx_contacted )
-{
-    double homophily = this->homophily;
-
-    if( idx_contacted == idx_contacter )
-        return 0.0;
-
-    if( bot_present() && idx_contacter < n_bots )
-        homophily = this->bot_homophily[idx_contacter];
-
-    constexpr double tolerance = 1e-10;
-    auto opinion_diff
-        = std::abs( network.agents[idx_contacter].data.opinion - network.agents[idx_contacted].data.opinion );
-    opinion_diff = std::max( tolerance, opinion_diff );
-
-    return std::pow( opinion_diff, -homophily );
-}
-
-void ActivityDrivenModel::get_agents_from_power_law()
-{
-    std::uniform_real_distribution<> dis_opinion( -1, 1 ); // Opinion initial values
-    power_law_distribution<> dist_activity( eps, gamma );
-    truncated_normal_distribution<> dist_reluctance( reluctance_mean, reluctance_sigma, reluctance_eps );
-
-    auto mean_activity = dist_activity.mean();
-
-    // Initial conditions for the opinions, initialize to [-1,1]
-    // The activities should be drawn from a power law distribution
-    for( size_t i = 0; i < network.agents.size(); i++ )
-    {
-        network.agents[i].data.opinion = dis_opinion( gen ); // Draw the opinion value
-
-        if( !mean_activities )
-        { // Draw from a power law distribution (1-gamma)/(1-eps^(1-gamma)) * a^(-gamma)
-            network.agents[i].data.activity = dist_activity( gen );
-        }
-        else
-        {
-            network.agents[i].data.activity = mean_activity;
-        }
-
-        if( use_reluctances )
-        {
-            network.agents[i].data.reluctance = dist_reluctance( gen );
-        }
-    }
-
-    if( bot_present() )
-    {
-        for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
-        {
-            network.agents[bot_idx].data.opinion  = bot_opinion[bot_idx];
-            network.agents[bot_idx].data.activity = bot_activity[bot_idx];
-        }
-    }
-}
-
-void ActivityDrivenModel::update_network_probabilistic()
-{
-    network.switch_direction_flag();
-
-    std::uniform_real_distribution<> dis_activation( 0.0, 1.0 );
-    std::uniform_real_distribution<> dis_reciprocation( 0.0, 1.0 );
-    std::vector<size_t> contacted_agents{};
-    reciprocal_edge_buffer.clear(); // Clear the reciprocal edge buffer
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Test if the agent is activated
-        bool activated = dis_activation( gen ) < network.agents[idx_agent].data.activity;
-
-        if( activated )
-        {
-            // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
-            // Not normalised since this is taken care of by the reservoir sampling
-
-            int m_temp = this->m;
-
-            if( bot_present() && idx_agent < n_bots )
-            {
-                m_temp = bot_m[idx_agent];
-            }
-
-            reservoir_sampling_A_ExpJ(
-                m_temp, network.n_agents(), [&]( int j ) { return homophily_weight( idx_agent, j ); }, contacted_agents,
-                gen );
-
-            // Fill the outgoing edges into the reciprocal edge buffer
-            for( const auto & idx_outgoing : contacted_agents )
-            {
-                reciprocal_edge_buffer.insert(
-                    { idx_agent, idx_outgoing } ); // insert the edge idx_agent -> idx_outgoing
-            }
-
-            // Set the *outgoing* edges
-            network.set_neighbours_and_weights( idx_agent, contacted_agents, 1.0 );
-        }
-        else
-        {
-            network.set_neighbours_and_weights( idx_agent, {}, {} );
-        }
-    }
-
-    // Reciprocity check
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Get the outgoing edges
-        auto contacted_agents = network.get_neighbours( idx_agent );
-        // For each outgoing edge we check if the reverse edge already exists
-        for( const auto & idx_outgoing : contacted_agents )
-        {
-            // If the edge is not reciprocated
-            if( !reciprocal_edge_buffer.contains( { idx_outgoing, idx_agent } ) )
-            {
-                if( dis_reciprocation( gen ) < reciprocity )
-                {
-                    network.push_back_neighbour_and_weight( idx_outgoing, idx_agent, 1.0 );
-                }
-            }
-        }
-    }
-
-    network.toggle_incoming_outgoing(); // switch direction, so that we have incoming edges
-}
-
-void ActivityDrivenModel::update_network_mean()
-{
-    using WeightT = NetworkT::WeightT;
-    std::vector<WeightT> weights( network.n_agents(), 0.0 );
-
-    // Set all weights to zero in the beginning
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        network.set_weights( idx_agent, weights );
-        contact_prob_list[idx_agent] = weights; // set to zero
-    }
-
-    auto probability_helper = []( double omega, size_t m )
-    {
-        double p = 0;
-        for( size_t i = 1; i <= m; i++ )
-            p += ( std::pow( -omega, i + 1 ) + omega ) / ( omega + 1 );
-        return p;
-    };
-
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Implement the weight for the probability of agent `idx_agent` contacting agent `j`
-        // Not normalised since this is taken care of by the reservoir sampling
-
-        double normalization = 0;
-        for( size_t k = 0; k < network.n_agents(); k++ )
-        {
-            normalization += homophily_weight( idx_agent, k );
-        }
-
-        // Go through all the neighbours of idx_agent
-        // Calculate the probability of i contacting j (in 1 to m rounds, assuming
-        // the agent is activated
-        int m_temp = m;
-        if( bot_present() && idx_agent < n_bots )
-        {
-            m_temp = bot_m[idx_agent];
-        }
-
-        double activity = std::max( 1.0, network.agents[idx_agent].data.activity );
-        for( size_t j = 0; j < network.n_agents(); j++ )
-        {
-            double omega = homophily_weight( idx_agent, j ) / normalization;
-            // We can calculate the probability of i contacting j ( i->j )
-            // Update contact prob_list (outgoing)
-            contact_prob_list[idx_agent][j] = activity * probability_helper( omega, m_temp );
-        }
-    }
-
-    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
-    {
-        // Calculate the actual weights and reciprocity
-        for( size_t j = 0; j < network.n_agents(); j++ )
-        {
-            double prob_contact_ij = contact_prob_list[idx_agent][j]; // outgoing probabilites
-            double prob_contact_ji = contact_prob_list[j][idx_agent];
-
-            // Set the incoming agent weight, j-i in weight list
-            double & win_ji = network.get_weights( j )[idx_agent];
-            win_ji += prob_contact_ij;
-
-            // Handle the reciprocity for j->i
-            // Update incoming weight i-j
-            double & win_ij = network.get_weights( idx_agent )[j];
-
-            // The probability of reciprocating is
-            win_ij += ( 1.0 - prob_contact_ji ) * reciprocity * prob_contact_ij;
-        }
-    }
-}
-
-void ActivityDrivenModel::update_network()
-{
-
-    if( !mean_weights )
-    {
-        update_network_probabilistic();
-    }
-    else
-    {
-        update_network_mean();
-    }
-}
-
-void ActivityDrivenModel::iteration()
+template<>
+void ActivityDrivenModelAbstract<ActivityAgent>::iteration()
 {
     Model<AgentT>::iteration();
 
@@ -265,19 +20,21 @@ void ActivityDrivenModel::iteration()
     get_euler_slopes( k1_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
     // k_2  =   hf(x_n+1/2h,y_n+1/2k_1)
     get_euler_slopes(
-        k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k1_buffer[i]; } );
+        k2_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k1_buffer[i]; } );
     // k_3  =   hf(x_n+1/2h,y_n+1/2k_2)
     get_euler_slopes(
-        k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * this->k2_buffer[i]; } );
+        k3_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + 0.5 * dt * this->k2_buffer[i]; } );
     // k_4  =   hf(x_n+h,y_n+k_3)
-    get_euler_slopes( k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + this->k3_buffer[i]; } );
+    get_euler_slopes(
+        k4_buffer, [this]( size_t i ) { return network.agents[i].data.opinion + dt * this->k3_buffer[i]; } );
 
     // Update the agent opinions
     for( size_t idx_agent = 0; idx_agent < network.n_agents(); ++idx_agent )
     {
         // y_(n+1) =   y_n+1/6k_1+1/3k_2+1/3k_3+1/6k_4+O(h^5)
         network.agents[idx_agent].data.opinion
-            += ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] )
+            += dt
+               * ( k1_buffer[idx_agent] + 2 * k2_buffer[idx_agent] + 2 * k3_buffer[idx_agent] + k4_buffer[idx_agent] )
                / 6.0;
     }
 
@@ -289,4 +46,6 @@ void ActivityDrivenModel::iteration()
         }
     }
 }
+
+template class ActivityDrivenModelAbstract<ActivityAgent>;
 } // namespace Seldon
\ No newline at end of file
diff --git a/src/models/InertialModel.cpp b/src/models/InertialModel.cpp
new file mode 100644
index 0000000..df46cc1
--- /dev/null
+++ b/src/models/InertialModel.cpp
@@ -0,0 +1,77 @@
+#include "models/InertialModel.hpp"
+#include "agents/inertial_agent.hpp"
+#include "network.hpp"
+#include "util/math.hpp"
+#include <cstddef>
+#include <random>
+#include <vector>
+
+namespace Seldon
+{
+
+InertialModel::InertialModel(
+    const Config::ActivityDrivenInertialSettings & settings, NetworkT & network, std::mt19937 & gen )
+        : ActivityDrivenModelAbstract<InertialAgent>( settings, network, gen ),
+          friction_coefficient( settings.friction_coefficient )
+{
+}
+
+// X(t+dt)
+// Also updates the position
+void InertialModel::calc_position()
+{
+    // Calculating 'drift' = a(t)-friction
+    get_euler_slopes( drift_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
+
+    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+    {
+        auto & agent_data    = network.agents[idx_agent].data;
+        auto accleration     = drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity;
+        double next_position = agent_data.opinion + agent_data.velocity * dt + 0.5 * (accleration)*dt * dt;
+
+        // Update the position to the new position
+        agent_data.opinion = next_position;
+    }
+}
+
+// V(t+dt)
+// Also updates the velocity
+void InertialModel::calc_velocity()
+{
+    // Calculating new 'drift'
+    get_euler_slopes( drift_next_t_buffer, [this]( size_t i ) { return network.agents[i].data.opinion; } );
+
+    for( size_t idx_agent = 0; idx_agent < network.n_agents(); idx_agent++ )
+    {
+        auto & agent_data    = network.agents[idx_agent].data;
+        double next_velocity = agent_data.velocity
+                               + 0.5 * dt
+                                     * ( drift_t_buffer[idx_agent] - friction_coefficient * agent_data.velocity
+                                         + drift_next_t_buffer[idx_agent] );
+        next_velocity /= 1.0 + 0.5 * friction_coefficient * dt;
+
+        // Update velocity
+        agent_data.velocity = next_velocity;
+    }
+}
+
+void InertialModel::iteration()
+{
+    Model<AgentT>::iteration();
+
+    update_network();
+
+    // Use Velocity Verlet algorithm
+    calc_position();
+    calc_velocity();
+
+    if( bot_present() )
+    {
+        for( size_t bot_idx = 0; bot_idx < n_bots; bot_idx++ )
+        {
+            network.agents[bot_idx].data.opinion = bot_opinion[bot_idx];
+        }
+    }
+}
+
+} // namespace Seldon
\ No newline at end of file
diff --git a/test/res/1bot_1agent_inertial.toml b/test/res/1bot_1agent_inertial.toml
new file mode 100644
index 0000000..6520309
--- /dev/null
+++ b/test/res/1bot_1agent_inertial.toml
@@ -0,0 +1,42 @@
+[simulation]
+model = "ActivityDrivenInertial"
+rng_seed = 120           # Leaving this empty will pick a random seed
+
+[io]
+# n_output_network = 1  # Write the network every 20 iterations
+# n_output_agents = 1
+print_progress = false # Print the iteration time ; if not set, then does not print
+
+[model]
+max_iterations = 1000 # If not set, max iterations is infinite
+
+[ActivityDrivenInertial]
+dt = 0.001              # Timestep for the integration of the coupled ODEs
+m = 1                   # Number of agents contacted, when the agent is active
+eps = 1                 # Minimum activity epsilon; a_i belongs to [epsilon,1]
+gamma = 2.1             # Exponent of activity power law distribution of activities
+reciprocity = 1         # probability that when agent i contacts j via weighted reservoir sampling, j also sends feedback to i. So every agent can have more than m incoming connections
+homophily = 0.5         # aka beta. if zero, agents pick their interaction partners at random
+alpha = 1.5             # Controversialness of the issue, must be greater than 0.
+K = 2.0                 # Social interaction strength
+mean_activities = false # Use the mean value of the powerlaw distribution for the activities of all agents
+mean_weights = false    # Use the meanfield approximation of the network edges
+
+reluctances = true     # Assigns a "reluctance" (m_i) to each agent. By default; false and every agent has a reluctance of 1
+reluctance_mean = 1.5  # Mean of distribution before drawing from a truncated normal distribution (default set to 1.0)
+reluctance_sigma = 0.1 # Width of normal distribution (before truncating)
+reluctance_eps = 0.01  # Minimum such that the normal distribution is truncated at this value
+
+n_bots = 1 # The number of bots to be used; if not specified defaults to 0 (which means bots are deactivated)
+# Bots are agents with fixed opinions and different parameters, the parameters are specified in the following lists
+# If n_bots is smaller than the length of any of the lists, the first n_bots entries are used. If n_bots is greater the code will throw an exception.
+bot_m = [1]           # If not specified, defaults to `m`
+bot_homophily = [0.7] # If not specified, defaults to `homophily`
+bot_activity = [1.0]  # If not specified, defaults to 0
+bot_opinion = [2]     # The fixed opinions of the bots
+
+friction_coefficient = 0.5 # Friction coefficient 
+
+[network]
+number_of_agents = 2
+connections_per_agent = 1
diff --git a/test/test_activity_inertial.cpp b/test/test_activity_inertial.cpp
new file mode 100644
index 0000000..60f3eef
--- /dev/null
+++ b/test/test_activity_inertial.cpp
@@ -0,0 +1,76 @@
+#include "catch2/matchers/catch_matchers.hpp"
+#include "models/InertialModel.hpp"
+#include "util/math.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <complex>
+#include <config_parser.hpp>
+#include <filesystem>
+#include <network_io.hpp>
+#include <optional>
+#include <simulation.hpp>
+
+namespace fs = std::filesystem;
+
+TEST_CASE( "Test the probabilistic inertial activity driven model with one bot and one agent", "[inertial1Bot1Agent]" )
+{
+    using namespace Seldon;
+    using namespace Catch::Matchers;
+    using AgentT = InertialModel::AgentT;
+
+    auto proj_root_path = fs::current_path();
+    auto input_file     = proj_root_path / fs::path( "test/res/1bot_1agent_inertial.toml" );
+
+    auto options = Config::parse_config_file( input_file.string() );
+    Config::print_settings( options );
+
+    auto simulation = Simulation<AgentT>( options, std::nullopt, std::nullopt );
+
+    // We need an output path for Simulation, but we won't write anything out there
+    fs::path output_dir_path = proj_root_path / fs::path( "test/output_inertial" );
+
+    // Get the bot opinion (which won't change)
+    auto bot   = simulation.network.agents[0];
+    auto x_bot = bot.data.opinion; // Bot opinion
+    fmt::print( "We have set the bot opinion = {}\n", x_bot );
+
+    // Get the initial agent opinion
+    auto & agent = simulation.network.agents[1];
+    auto x_0     = agent.data.opinion; // Agent opinion
+    fmt::print( "We have set agent x_0 = {}\n", x_0 );
+
+    simulation.run( output_dir_path );
+
+    auto model_settings     = std::get<Seldon::Config::ActivityDrivenInertialSettings>( options.model_settings );
+    auto K                  = model_settings.K;
+    auto alpha              = model_settings.alpha;
+    auto iterations         = model_settings.max_iterations.value();
+    auto dt                 = model_settings.dt;
+    std::complex<double> mu = model_settings.friction_coefficient;
+    auto time_elapsed       = iterations * dt;
+
+    // Final agent and bot opinions after the simulation run
+    auto x_t        = agent.data.opinion;
+    auto x_t_bot    = bot.data.opinion;
+    auto reluctance = agent.data.reluctance;
+
+    fmt::print( "Agent reluctance is = {}\n", reluctance );
+
+    // The bot opinion should not change during the simulation
+    REQUIRE_THAT( x_t_bot, WithinAbs( x_bot, 1e-16 ) );
+
+    // C = K/m tanh (alpha*x_bot)
+    auto C = K / reluctance * std::tanh( alpha * x_bot );
+
+    fmt::print( "C =     {}\n", C );
+    auto a1 = 0.5 * ( -std::sqrt( mu * mu - 4.0 ) - mu );
+    auto a2 = 0.5 * ( std::sqrt( mu * mu - 4.0 ) - mu );
+    auto c1 = ( x_0 - C ) / ( 1.0 - a1 / a2 );
+    auto c2 = -c1 * a1 / a2;
+    // Test that the agent opinion matches the analytical solution for an agent with a bot
+    // Analytical solution is
+    // x_t = c1 * exp(a1*t) + c2 *exp(a2*t) + C
+    auto x_t_analytical = c1 * std::exp( a1 * time_elapsed ) + c2 * std::exp( a2 * time_elapsed ) + C;
+
+    REQUIRE_THAT( x_t, WithinAbs( x_t_analytical.real(), 1e-5 ) );
+}
\ No newline at end of file
diff --git a/test/test_probability_distributions.cpp b/test/test_probability_distributions.cpp
new file mode 100644
index 0000000..e72a3df
--- /dev/null
+++ b/test/test_probability_distributions.cpp
@@ -0,0 +1,55 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
+
+#include "util/math.hpp"
+#include <fmt/ostream.h>
+#include <algorithm>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <random>
+namespace fs = std::filesystem;
+
+template<typename T, std::size_t N>
+std::ostream & operator<<( std::ostream & os, std::array<T, N> const & v1 )
+{
+    std::for_each( begin( v1 ), end( v1 ), [&os]( T val ) { os << val << " "; } );
+    return os;
+}
+
+// Samples the distribution n_samples times and writes results to file
+template<typename distT>
+void write_results_to_file( int N_Samples, distT dist, const std::string & filename )
+{
+    auto proj_root_path = fs::current_path();
+    auto file           = proj_root_path / fs::path( "test/output_probability_distributions/" + filename );
+    fmt::print( "file = {}\n", fmt::streamed( file ) );
+    fs::create_directories( file.parent_path() );
+
+    auto gen = std::mt19937( 0 );
+
+    std::ofstream filestream( file );
+    filestream << std::setprecision( 16 );
+
+    for( int i = 0; i < N_Samples; i++ )
+    {
+        filestream << dist( gen ) << "\n";
+    }
+    filestream.close();
+}
+
+TEST_CASE( "Test the probability distributions", "[prob]" )
+{
+    write_results_to_file( 10000, Seldon::truncated_normal_distribution( 1.0, 0.5, 0.75 ), "truncated_normal.txt" );
+    write_results_to_file( 10000, Seldon::power_law_distribution( 0.01, 2.1 ), "power_law.txt" );
+    write_results_to_file( 10000, Seldon::bivariate_normal_distribution( 0.5 ), "bivariate_normal.txt" );
+}
+
+TEST_CASE( "Test bivariate gaussian copula", "[prob_copula]" )
+{
+    auto dist1  = Seldon::power_law_distribution( 0.02, 2.5 );
+    auto dist2  = Seldon::truncated_normal_distribution( 1.0, 0.75, 0.2 );
+    auto copula = Seldon::bivariate_gaussian_copula( 0.5, dist1, dist2 );
+    write_results_to_file( 10000, copula, "gaussian_copula.txt" );
+}