diff --git a/collective_effects/__init__.py b/collective_effects/__init__.py
index 39ac40292aebe5ae7dadd9859d8a99a319dfb1a2..b680bc03a42065f9f5f210a1a6a070edfee7b716 100644
--- a/collective_effects/__init__.py
+++ b/collective_effects/__init__.py
@@ -5,7 +5,7 @@ Created on Tue Jan 14 12:25:30 2020
 @author: gamelina
 """
 
-from mbtrack2.collective_effects.instabilities import mbi_threshold, cbi_threshold
+from mbtrack2.collective_effects.instabilities import mbi_threshold, cbi_threshold, lcbi_growth_rate_mode, lcbi_growth_rate, plot_critical_mass
 from mbtrack2.collective_effects.resistive_wall import (skin_depth, CircularResistiveWall, Coating)
 from mbtrack2.collective_effects.resonator import Resonator, PureInductive, PureResistive
 from mbtrack2.collective_effects.tapers import StupakovRectangularTaper, StupakovCircularTaper
diff --git a/collective_effects/instabilities.py b/collective_effects/instabilities.py
index 12598cae3576dde196f51d92b812735be6285bc8..1fe58eca44183583da2580d3ef19b09945ca690a 100644
--- a/collective_effects/instabilities.py
+++ b/collective_effects/instabilities.py
@@ -3,7 +3,6 @@
 General calculations about instabilities
 
 @author: Alexis Gamelin
-@date: 15/01/2020
 """
 
 import numpy as np
@@ -88,6 +87,120 @@ def cbi_threshold(ring, I, Vrf, f, beta, Ncav=1):
     
     return (Zlong, Zxdip, Zydip)
 
+def lcbi_growth_rate_mode(ring, I, Vrf, M, mu, fr=None, Rs=None, QL=None, Z=None):
+    """
+    Compute the longitudinal coupled bunch instability growth rate driven by
+    HOMs for a given coupled bunch mode mu [1].
+    
+    Use either the resonator model (fr,Rs,QL) or an Impedance object (Z).
+
+    Parameters
+    ----------
+    ring : Synchrotron object
+    I : float
+        Total beam current in [A].
+    Vrf : float
+        Total RF voltage in [V].
+    M : int
+        Nomber of bunches in the beam.
+    mu : int
+        Coupled bunch mode number (= 0, ..., M-1).
+    fr : float, optional
+        Frequency of the HOM in [Hz].
+    Rs : float, optional
+        Shunt impedance of the HOM in [Ohm].
+    QL : float, optional
+        Loaded quality factor of the HOM.
+    Z : Impedance, optional
+        Longitunial impedance to consider.
+
+    Returns
+    -------
+    float
+        Coupled bunch instability growth rate for the mode mu.
+        
+    References
+    ----------
+    [1] : Eq. 51 p139 of Akai, Kazunori. "RF System for Electron Storage 
+    Rings." Physics And Engineering Of High Performance Electron Storage Rings 
+    And Application Of Superconducting Technology. 2002. 118-149.
+
+    """
+
+    nu_s = ring.synchrotron_tune(Vrf)
+    factor = ring.eta * I / (4 * np.pi * ring.E0 * nu_s)
+    
+    if Z is None:
+        omega_r = 2 * np.pi * fr
+        n_max = int(10 * omega_r / (ring.omega0 * M))
+        def Zr(omega):
+            return np.real(Rs / (1 + 1j * QL * (omega_r/omega - omega/omega_r)))
+    else:
+        fmax = Z.data.index.max()
+        n_max = int(2 * np.pi * fmax / (ring.omega0 * M))
+        def Zr(omega):
+            return np.real( Z( omega / (2*np.pi) ) )
+        
+    n0 = np.arange(n_max)
+    n1 = np.arange(1, n_max)
+    omega_p = ring.omega0 * (n0 * M + mu + nu_s)
+    omega_m = ring.omega0 * (n1 * M - mu - nu_s)
+        
+    sum_val = np.sum(omega_p*Zr(omega_p)) - np.sum(omega_m*Zr(omega_m))
+
+    return factor * sum_val
+    
+def lcbi_growth_rate(ring, I, Vrf, M, fr=None, Rs=None, QL=None, Z=None):
+    """
+    Compute the maximum growth rate for longitudinal coupled bunch instability 
+    driven by HOMs [1].
+    
+    Use either the resonator model (fr,Rs,QL) or an Impedance object (Z).
+
+    Parameters
+    ----------
+    ring : Synchrotron object
+    I : float
+        Total beam current in [A].
+    Vrf : float
+        Total RF voltage in [V].
+    M : int
+        Nomber of bunches in the beam.
+    fr : float, optional
+        Frequency of the HOM in [Hz].
+    Rs : float, optional
+        Shunt impedance of the HOM in [Ohm].
+    QL : float, optional
+        Loaded quality factor of the HOM.
+    Z : Impedance, optional
+        Longitunial impedance to consider.
+
+    Returns
+    -------
+    growth_rate : float
+        Maximum coupled bunch instability growth rate.
+    mu : int
+        Coupled bunch mode number corresponding to the maximum coupled bunch 
+        instability growth rate.
+    growth_rates : array
+        Coupled bunch instability growth rates for the different mode numbers.
+        
+    References
+    ----------
+    [1] : Eq. 51 p139 of Akai, Kazunori. "RF System for Electron Storage 
+    Rings." Physics And Engineering Of High Performance Electron Storage Rings 
+    And Application Of Superconducting Technology. 2002. 118-149.
+
+    """
+    growth_rates = np.zeros(M)
+    for i in range(M):
+        growth_rates[i] = lcbi_growth_rate_mode(ring, I, Vrf, M, i, fr=fr, Rs=Rs, QL=QL, Z=Z)
+    
+    growth_rate = np.max(growth_rates)
+    mu = np.argmax(growth_rates)
+    
+    return growth_rate, mu, growth_rates
+
 def plot_critical_mass(ring, bunch_charge, bunch_spacing, n_points=1e4):
     """
     Plot ion critical mass, using Eq. (7.70) p147 of [1]