diff --git a/src/analysis/runup_models.py b/src/analysis/runup_models.py
index 403768c..650c36b 100644
--- a/src/analysis/runup_models.py
+++ b/src/analysis/runup_models.py
@@ -1,6 +1,7 @@
 import numpy as np
+import pandas as pd
 
-def stockdon06(Hs0, Tp, beta):
+def sto06_individual(Hs0, Tp, beta):
 
     Lp = 9.8 * Tp ** 2 / 2 / np.pi
 
@@ -19,5 +20,32 @@ def stockdon06(Hs0, Tp, beta):
 
     return R2, setup, S_total, S_inc, S_ig
 
+def sto06(df, Hs0_col, Tp_col, beta_col):
+    """
+    Vectorized version of Stockdon06 which can be used with dataframes
+    :param df:
+    :param Hs0_col:
+    :param Tp_col:
+    :param beta_col:
+    :return:
+    """
+
+    Lp  = 9.8 * df[Tp_col] ** 2 / 2 / np.pi
+
+    # General equation
+    S_ig = pd.to_numeric(0.06 * np.sqrt(df[Hs0_col] * Lp), errors='coerce')
+    S_inc = pd.to_numeric(0.75 * df[beta_col] * np.sqrt(df[Hs0_col] * Lp), errors='coerce')
+    setup = pd.to_numeric(0.35 * df[beta_col] * np.sqrt(df[Hs0_col] * Lp), errors='coerce')
+    S_total = np.sqrt(S_inc ** 2 + S_ig ** 2)
+    R2 = 1.1 * (setup + S_total / 2)
+
+    # Dissipative conditions
+    dissipative = df[beta_col] / (df[Hs0_col] / Lp)**(0.5) <= 0.3
+    setup.loc[dissipative,:] = 0.016 * (df[Hs0_col] * Lp) ** (0.5)  # eqn 16
+    S_total.loc[dissipative,:] = 0.046 * (df[Hs0_col] * Lp) ** (0.5)  # eqn 17
+    R2.loc[dissipative,:] = 0.043 * (df[Hs0_col] * Lp) ** (0.5)  # eqn 18
+
+    return R2, setup, S_total, S_inc, S_ig
+
 if __name__ == '__main__':
     pass