diff --git a/README.md b/README.md
index 866fdc0..30aa32f 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,10 @@
 # HOA
-SUTD 2020 10.008 Hands-on Activity 1 Plotting Code Dump
+SUTD 2020 10.008 Hands-on Activity Code Dump
 
-This repo serves as a personal code dump of the script used to generate the answers and the resulting graphs for the Thermal Physics HOA 1 assignment.
+This repo serves as a personal code dump of the script used to generate the answers and the resulting graphs/plots for the Thermal Physics HOA 1 and HOA 3 assignments.
 
 ## Usage
 
-Add in your Plotly username and API key into the script and run by executing `python3 hoa1.py`.
\ No newline at end of file
+For the HOA 1 script, add in your Plotly username and API key into the script and run by executing `python3 hoa1.py`.
+
+For the HOA 3 script, just run `python3 hoa3.py`.
\ No newline at end of file
diff --git a/hoa1.py b/hoa1.py
deleted file mode 100644
index 71ffa3a..0000000
--- a/hoa1.py
+++ /dev/null
@@ -1,371 +0,0 @@
-# 10.008 Hands-on Activity 1 Plotting Script
-# This will spawn application window and browser tab popups to display the graphs
-# NOTE: Ignore the pandas.util.testing and plotly.graph_objs.Marker deprecation warnings in the console when running this script
-# Created by James Raphael Tiovalen (2020)
-
-# Setup plotly
-import chart_studio
-chart_studio.tools.set_credentials_file(username='', api_key='')
-chart_studio.tools.set_config_file(world_readable=False, sharing='private')
-import chart_studio.plotly as py
-import plotly.graph_objects as go
-import plotly.io as pio
-
-# Import scientific libraries
-import pandas as pd
-import numpy as np
-from scipy.optimize import curve_fit
-from sklearn import metrics
-from sympy import symbols, diff, integrate, lambdify, exp
-import matplotlib.pyplot as plt
-import seaborn as sns
-
-# Define exponential function to fit to
-def exponential(x, a, b):
-    return a * np.exp(- b * (x - 61.9))
-
-# Define output power function
-def output(rl, rin, vse):
-    return ((vse ** 2) * rl) / ((rin + rl) ** 2)
-
-# Plot exponential graphs
-def plot_exponential_graph_set(time, t_c, t_h, power, resistance, room_temp):
-    data_dict = {"T_c": (t_c, "#17becf", 40, 20), "T_h": (t_h, "#de1738", 40, -20), "\dot{W}": (power, "#7851a9", 40, -20)}
-
-    scatters = []
-    trendlines = []
-    annotations = []
-
-    for dataname, datapair in data_dict.items():
-        # Use dogleg algorithm
-        popt, pcov = curve_fit(exponential, time, datapair[0], method="dogbox")
-        
-        # Prep constants for return values
-        if dataname == "T_c":
-            tc_coeffs = popt
-        elif dataname == "T_h":
-            th_coeffs = popt
-        elif dataname == "\dot{W}":
-            w_coeffs = popt
-
-        # Compute R^2 value
-        r_squared = round(metrics.r2_score(exponential(time, *popt), datapair[0]), 3)
-
-        xx = np.linspace(time[0], time[-1], 2000)
-        yy = exponential(xx, *popt)
-        
-        if dataname == "T_c" or dataname == "T_h":
-            # Creating the dataset, and generating the plot
-            scatters.append(go.Marker(
-                              x=time,
-                              y=datapair[0] + room_temp,
-                              mode='markers',
-                              marker=go.Marker(color=list(range(140)), colorscale="Viridis"),
-                              name="",
-                              showlegend=False
-                              ))
-
-            trendlines.append(go.Marker(
-                              x=xx,
-                              y=yy + room_temp,
-                              mode='lines',
-                              marker=go.Marker(color=datapair[1]),
-                              name='${{{}}}$'.format(dataname),
-                              ))
-
-            annotations.append(go.layout.Annotation(
-                              x=time[-20],
-                              y=exponential(time[-20], *popt) + room_temp,
-                              text='${} = {}e^{{-{}(t - 61.9)}} + {{{}}}$'.format(dataname, round(popt[0], 3), round(popt[1], 3), room_temp),
-                              showarrow=True,
-                              ax=datapair[2],
-                              ay=datapair[3],
-                              ))
-
-            annotations.append(go.layout.Annotation(
-                              x=time[-1],
-                              y=exponential(time[-1], *popt) + room_temp,
-                              text='$R^2 = {}$'.format(r_squared),
-                              showarrow=True,
-                              ax=datapair[2],
-                              ay=datapair[3],
-                              ))
-
-        elif dataname == "\dot{W}":
-            scatters.append(go.Marker(
-                              x=time,
-                              y=datapair[0],
-                              mode='markers',
-                              marker=go.Marker(color=list(range(140)), colorscale="Viridis"),
-                              name="",
-                              showlegend=False
-                              ))
-
-            trendlines.append(go.Marker(
-                              x=xx,
-                              y=yy,
-                              mode='lines',
-                              marker=go.Marker(color=datapair[1]),
-                              name='${{{}}}$'.format(dataname),
-                              ))
-
-            annotations.append(go.layout.Annotation(
-                              x=time[-20],
-                              y=exponential(time[-20], *popt),
-                              text='${} = {}e^{{-{}(t - 61.9)}}$'.format(dataname, round(popt[0], 3), round(popt[1], 3)),
-                              showarrow=True,
-                              ax=datapair[2],
-                              ay=datapair[3],
-                              ))
-
-            annotations.append(go.layout.Annotation(
-                              x=time[-1],
-                              y=exponential(time[-1], *popt),
-                              text='$R^2 = {}$'.format(r_squared),
-                              showarrow=True,
-                              ax=datapair[2],
-                              ay=datapair[3],
-                              ))
-
-    layout = go.Layout(
-                title='$\\textbf{{Graph for }} {} \\textbf{{ }} \Omega$'.format(resistance),
-                plot_bgcolor='rgb(229, 229, 229)',
-                xaxis=go.layout.XAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$t (s)$'),
-                yaxis=go.layout.YAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$T (^\circ C)$'),
-                annotations=annotations,
-                )
-
-    scatter_trendline_pair = list(zip(scatters, trendlines))
-    final_data = [item for sublist in scatter_trendline_pair for item in sublist]
-    fig = go.Figure(data=final_data, layout=layout)
-
-    # Show plot on localhost
-    pio.show(fig)
-    # pio.write_html(fig, file='T_c vs Time', auto_open=True)
-    return tc_coeffs, th_coeffs, w_coeffs
-
-def plot_qc_qh_and_w(time, tc_coeffs, th_coeffs, resistance, room_temp):
-    # Define symbols for differentiation
-    t, Tc, Th = symbols("t Tc Th")
-    Tc = tc_coeffs[0] * exp(-tc_coeffs[1] * (t - 61.9)) + room_temp
-    Th = th_coeffs[0] * exp(-th_coeffs[1] * (t - 61.9)) + room_temp
-
-    dTcdt = diff(Tc, t)
-    dThdt = diff(Th, t)
-
-    # mc = 21.87 J/°C
-    dQcdt = 21.87 * dTcdt
-    dQhdt = 21.87 * dThdt
-    dWdt = dQhdt + dQcdt
-    actual_Qc = lambdify(t, dQcdt)
-    actual_Qh = lambdify(t, dQhdt)
-    actual_W = lambdify(t, dWdt)
-
-    # Start plotting
-    data_dict = {"\dot{Q_c}": (tc_coeffs[1], actual_Qc, "#17becf", dQcdt, -40, 20), "\dot{Q_h}": (th_coeffs[1], actual_Qh, "#de1738", dQhdt, 40, -20), "\dot{W}": (th_coeffs[1] + tc_coeffs[1], actual_W, "#7851a9", dWdt, 40, 20)}
-    plots = []
-    annotations = []
-
-    for dataname, datapair in data_dict.items():
-        xx = np.linspace(time[0], time[-1], 2000)
-        yy = abs(datapair[1](xx))
-
-        plots.append(go.Marker(
-                          x=xx,
-                          y=yy,
-                          mode='lines',
-                          marker=go.Marker(color=datapair[2]),
-                          name='${{{}}}$'.format(dataname),
-                          ))
-
-        if dataname == "\dot{Q_c}" or dataname == "\dot{Q_h}":
-            annotations.append(go.layout.Annotation(
-                          x=time[20],
-                          y=abs(datapair[1](time[20])),
-                          text='${} = {}e^{{-{}t}}$'.format(dataname, abs(round(datapair[3].args[0], 3)), round(datapair[0], 3)),
-                          showarrow=True,
-                          ax=datapair[4],
-                          ay=datapair[5],
-                          ))
-
-        elif dataname == "\dot{W}":
-            annotations.append(go.layout.Annotation(
-                          x=time[-20],
-                          y=datapair[1](time[-20]),
-                          text='${} = {}e^{{-{}t}} - {}e^{{-{}t}}$'.format(dataname, abs(round(dQhdt.args[0], 3)), round(th_coeffs[1], 3), abs(round(dQcdt.args[0], 3)), round(tc_coeffs[1], 3)),
-                          showarrow=True,
-                          ax=datapair[4],
-                          ay=datapair[5],
-                          ))
-
-    layout = go.Layout(
-                title='$\\textbf{{Graph for }} {} \\textbf{{ }} \Omega$'.format(resistance),
-                plot_bgcolor='rgb(229, 229, 229)',
-                xaxis=go.layout.XAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$t (s)$'),
-                yaxis=go.layout.YAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$\dot{Q} (J/s)$'),
-                annotations=annotations,
-                )
-
-    fig = go.Figure(data=plots, layout=layout)
-
-    # Show plot on localhost
-    pio.show(fig)
-
-# Trapezoidal rule method
-def stage_1_sum(stage_1_time, stage_1_power, resistance):
-    total_energy = 0
-    for i in range(len(stage_1_time) - 1):
-        total_energy += ((stage_1_power[i] + stage_1_power[i+1]) / 2) * (stage_1_time[i+1] - stage_1_time[i])
-    
-    return total_energy
-
-# Integration method
-def stage_2_sum(time, w_coeffs, resistance):
-    t, W = symbols("t W")
-    W = w_coeffs[0] * exp(-w_coeffs[1] * (t - 61.9))
-
-    IWdt = integrate(W, (t, time[0], time[-1]))
-    
-    return IWdt
-
-# Plot boxplot to detect anomalies
-def detect_anomalies(output_power_values):
-    sns.set(style="whitegrid")
-    ax = sns.boxplot(x=output_power_values, orient="h", palette="Set3", flierprops=dict(markerfacecolor='r', markersize=5, linestyle='none'))
-    ax.set(xlabel="Output power, $\dot{W}$")
-    plt.show()
-
-# Plot actual boxplot
-def plot_actual_boxplot(actual_output_power):
-    sns.set(style="whitegrid")
-    ax = sns.boxplot(x=actual_output_power, orient="h", palette="Set3", flierprops=dict(markerfacecolor='r', markersize=5, linestyle='none'))
-    ax.set(xlabel="Actual output power, $\dot{W}$")
-    plt.show()
-
-# Manual plotting
-def output_power_graph(resistance_values, output_power_values):
-    scatters = []
-    trendlines = []
-    annotations = []
-    
-    # Use dogleg algorithm
-    popt, pcov = curve_fit(output, resistance_values, output_power_values, method="dogbox")
-
-    # Compute R^2 value
-    r_squared = round(metrics.r2_score(output(resistance_values, *popt), output_power_values), 3)
-
-    xx = np.linspace(resistance_values[0], resistance_values[-1], 2000)
-    yy = output(xx, *popt)
-
-    scatters.append(go.Marker(
-                      x=resistance_values,
-                      y=output_power_values,
-                      mode='markers',
-                      marker=go.Marker(color=list(range(10)), colorscale="Viridis"),
-                      name="",
-                      showlegend=False
-                      ))
-
-    trendlines.append(go.Marker(
-                      x=xx,
-                      y=yy,
-                      mode='lines',
-                      marker=go.Marker(color="#ff007f"),
-                      name='${{\dot{W}}}$',
-                      ))
-
-    annotations.append(go.layout.Annotation(
-                      x=resistance_values[-4],
-                      y=output(resistance_values[-4], *popt),
-                      text='$\dot{{W}} = \\frac{{{}^2 R_L}}{{{{{} + R_L}}^2}}$'.format(round(popt[1], 3), round(popt[0], 3)),
-                      showarrow=True,
-                      ax=70,
-                      ay=-5,
-                      ))
-
-    annotations.append(go.layout.Annotation(
-                      x=resistance_values[-1],
-                      y=output(resistance_values[-1], *popt),
-                      text='$R^2 = {}$'.format(r_squared),
-                      showarrow=True,
-                      ax=70,
-                      ay=-5,
-                      ))
-
-    layout = go.Layout(
-                title='$\\textbf{Output Power Graph}$',
-                plot_bgcolor='rgb(229, 229, 229)',
-                xaxis=go.layout.XAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$R_L (\Omega)$'),
-                yaxis=go.layout.YAxis(zerolinecolor='rgb(255,255,255)', gridcolor='rgb(255,255,255)', title='$\dot{W} (W)$'),
-                annotations=annotations,
-                )
-
-    scatter_trendline_pair = list(zip(scatters, trendlines))
-    final_data = [item for sublist in scatter_trendline_pair for item in sublist]
-    fig = go.Figure(data=final_data, layout=layout)
-
-    # Show plot on localhost
-    pio.show(fig)
-
-# Function to graph (Ẇ, T_c and T_h) vs time
-def plot_data_graphs(resistance):
-    # Read xlsx file as dataframe
-    df = pd.read_excel("dataset.xlsx", sheet_name=f"{resistance}ohm", usecols = "A:D", na_values=" ", skiprows=range(0, 67)).dropna()
-    
-    stage_1 = pd.read_excel("dataset.xlsx", sheet_name=f"{resistance}ohm", usecols = "A:D", na_values=" ", skiprows=[0], nrows=64).dropna()
-    # Voltage is 3.0 V
-    stage_1_power = stage_1[" Current(A)"].to_numpy() * 3
-    stage_1_time = stage_1["time(s)"].to_numpy()
-
-    # Get room temperature based on average of last two data (to minimise systematic error)
-    room_temp = (df["CH1(C)"].to_numpy()[-1] + df["CH2(C)"].to_numpy()[-1]) / 2
-
-    # Parse columns as numpy.ndarray
-    time = df["time(s)"].to_numpy()
-    t_c = df["CH1(C)"].to_numpy() - room_temp
-    t_h = df["CH2(C)"].to_numpy() - room_temp
-    power = df[" Voltage(V)"].to_numpy() ** 2 / resistance
-
-    time = time.astype(np.float64)
-    t_c = t_c.astype(np.float64)
-    t_h = t_h.astype(np.float64)
-    power = power.astype(np.float64)
-
-    # Answer question (4)
-    tc_coeffs, th_coeffs, w_coeffs = plot_exponential_graph_set(time, t_c, t_h, power, resistance, room_temp)
-    # Answer question (5)
-    plot_qc_qh_and_w(time, tc_coeffs, th_coeffs, resistance, room_temp)
-    # Answer question (7)
-    input = stage_1_sum(stage_1_time, stage_1_power, resistance)
-    print(f"Energy input: {input} J")
-    # Answer question (8)
-    output = stage_2_sum(time, w_coeffs, resistance)
-    print(f"Energy output: {output} J")
-    # Answer question (9)
-    efficiency = output/input * 100
-    print(f"Efficiency: {efficiency}%")
-
-if __name__ == '__main__':
-    for resistance in [1]:
-        plot_data_graphs(resistance)
-    
-    # Manually define values for ΔT = 16.2 °C
-    resistance_values = np.array([1, 1.2, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6]).astype(np.float64)
-    voltage_values = np.array([0.176193, 0.193620666666667, 0.247881, 0.3116985, 0.317917, 0.3412875, 0.361429, 0.384375, 0.40167875]).astype(np.float64)
-    output_power_values = voltage_values ** 2 / resistance_values
-    output_power_values = output_power_values.astype(np.float64)
-    
-    # Evidence/justification to remove anomalies
-    detect_anomalies(output_power_values)
-    
-    # Remove anomaly
-    actual_resistances = np.delete(resistance_values, [3])
-    actual_voltages = np.delete(voltage_values, [3])
-    
-    actual_output_power = actual_voltages ** 2 / actual_resistances
-    actual_output_power = actual_output_power.astype(np.float64)
-    
-    plot_actual_boxplot(actual_output_power)
-    
-    # Answer question (11)
-    output_power_graph(actual_resistances, actual_output_power)
diff --git a/dataset.xlsx b/hoa1/dataset.xlsx
similarity index 100%
rename from dataset.xlsx
rename to hoa1/dataset.xlsx
diff --git a/experiment_1.xlsx b/hoa1/experiment_1.xlsx
similarity index 100%
rename from experiment_1.xlsx
rename to hoa1/experiment_1.xlsx
diff --git a/hoa1/hoa1.py b/hoa1/hoa1.py
new file mode 100644
index 0000000..8a1a2ba
--- /dev/null
+++ b/hoa1/hoa1.py
@@ -0,0 +1,505 @@
+# 10.008 Hands-on Activity 1 Plotting Script
+# This will spawn application window and browser tab popups to display the graphs
+# NOTE: Ignore the pandas.util.testing and plotly.graph_objs.Marker deprecation warnings in the console when running this script
+# Created by James Raphael Tiovalen (2020)
+
+# Setup plotly
+import chart_studio
+
+chart_studio.tools.set_credentials_file(username="", api_key="")
+chart_studio.tools.set_config_file(world_readable=False, sharing="private")
+import chart_studio.plotly as py
+import plotly.graph_objects as go
+import plotly.io as pio
+
+# Import scientific libraries
+import pandas as pd
+import numpy as np
+from scipy.optimize import curve_fit
+from sklearn import metrics
+from sympy import symbols, diff, integrate, lambdify, exp
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+# Define exponential function to fit to
+def exponential(x, a, b):
+    return a * np.exp(-b * (x - 61.9))
+
+
+# Define output power function
+def output(rl, rin, vse):
+    return ((vse ** 2) * rl) / ((rin + rl) ** 2)
+
+
+# Plot exponential graphs
+def plot_exponential_graph_set(time, t_c, t_h, power, resistance, room_temp):
+    data_dict = {
+        "T_c": (t_c, "#17becf", 40, 20),
+        "T_h": (t_h, "#de1738", 40, -20),
+        "\dot{W}": (power, "#7851a9", 40, -20),
+    }
+
+    scatters = []
+    trendlines = []
+    annotations = []
+
+    for dataname, datapair in data_dict.items():
+        # Use dogleg algorithm
+        popt, pcov = curve_fit(exponential, time, datapair[0], method="dogbox")
+
+        # Prep constants for return values
+        if dataname == "T_c":
+            tc_coeffs = popt
+        elif dataname == "T_h":
+            th_coeffs = popt
+        elif dataname == "\dot{W}":
+            w_coeffs = popt
+
+        # Compute R^2 value
+        r_squared = round(metrics.r2_score(exponential(time, *popt), datapair[0]), 3)
+
+        xx = np.linspace(time[0], time[-1], 2000)
+        yy = exponential(xx, *popt)
+
+        if dataname == "T_c" or dataname == "T_h":
+            # Creating the dataset, and generating the plot
+            scatters.append(
+                go.Marker(
+                    x=time,
+                    y=datapair[0] + room_temp,
+                    mode="markers",
+                    marker=go.Marker(color=list(range(140)), colorscale="Viridis"),
+                    name="",
+                    showlegend=False,
+                )
+            )
+
+            trendlines.append(
+                go.Marker(
+                    x=xx,
+                    y=yy + room_temp,
+                    mode="lines",
+                    marker=go.Marker(color=datapair[1]),
+                    name="${{{}}}$".format(dataname),
+                )
+            )
+
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[-20],
+                    y=exponential(time[-20], *popt) + room_temp,
+                    text="${} = {}e^{{-{}(t - 61.9)}} + {{{}}}$".format(
+                        dataname, round(popt[0], 3), round(popt[1], 3), room_temp
+                    ),
+                    showarrow=True,
+                    ax=datapair[2],
+                    ay=datapair[3],
+                )
+            )
+
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[-1],
+                    y=exponential(time[-1], *popt) + room_temp,
+                    text="$R^2 = {}$".format(r_squared),
+                    showarrow=True,
+                    ax=datapair[2],
+                    ay=datapair[3],
+                )
+            )
+
+        elif dataname == "\dot{W}":
+            scatters.append(
+                go.Marker(
+                    x=time,
+                    y=datapair[0],
+                    mode="markers",
+                    marker=go.Marker(color=list(range(140)), colorscale="Viridis"),
+                    name="",
+                    showlegend=False,
+                )
+            )
+
+            trendlines.append(
+                go.Marker(
+                    x=xx,
+                    y=yy,
+                    mode="lines",
+                    marker=go.Marker(color=datapair[1]),
+                    name="${{{}}}$".format(dataname),
+                )
+            )
+
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[-20],
+                    y=exponential(time[-20], *popt),
+                    text="${} = {}e^{{-{}(t - 61.9)}}$".format(
+                        dataname, round(popt[0], 3), round(popt[1], 3)
+                    ),
+                    showarrow=True,
+                    ax=datapair[2],
+                    ay=datapair[3],
+                )
+            )
+
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[-1],
+                    y=exponential(time[-1], *popt),
+                    text="$R^2 = {}$".format(r_squared),
+                    showarrow=True,
+                    ax=datapair[2],
+                    ay=datapair[3],
+                )
+            )
+
+    layout = go.Layout(
+        title="$\\textbf{{Graph for }} {} \\textbf{{ }} \Omega$".format(resistance),
+        plot_bgcolor="rgb(229, 229, 229)",
+        xaxis=go.layout.XAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$t (s)$",
+        ),
+        yaxis=go.layout.YAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$T (^\circ C)$",
+        ),
+        annotations=annotations,
+    )
+
+    scatter_trendline_pair = list(zip(scatters, trendlines))
+    final_data = [item for sublist in scatter_trendline_pair for item in sublist]
+    fig = go.Figure(data=final_data, layout=layout)
+
+    # Show plot on localhost
+    pio.show(fig)
+    # pio.write_html(fig, file='T_c vs Time', auto_open=True)
+    return tc_coeffs, th_coeffs, w_coeffs
+
+
+def plot_qc_qh_and_w(time, tc_coeffs, th_coeffs, resistance, room_temp):
+    # Define symbols for differentiation
+    t, Tc, Th = symbols("t Tc Th")
+    Tc = tc_coeffs[0] * exp(-tc_coeffs[1] * (t - 61.9)) + room_temp
+    Th = th_coeffs[0] * exp(-th_coeffs[1] * (t - 61.9)) + room_temp
+
+    dTcdt = diff(Tc, t)
+    dThdt = diff(Th, t)
+
+    # mc = 21.87 J/°C
+    dQcdt = 21.87 * dTcdt
+    dQhdt = 21.87 * dThdt
+    dWdt = dQhdt + dQcdt
+    actual_Qc = lambdify(t, dQcdt)
+    actual_Qh = lambdify(t, dQhdt)
+    actual_W = lambdify(t, dWdt)
+
+    # Start plotting
+    data_dict = {
+        "\dot{Q_c}": (tc_coeffs[1], actual_Qc, "#17becf", dQcdt, -40, 20),
+        "\dot{Q_h}": (th_coeffs[1], actual_Qh, "#de1738", dQhdt, 40, -20),
+        "\dot{W}": (th_coeffs[1] + tc_coeffs[1], actual_W, "#7851a9", dWdt, 40, 20),
+    }
+    plots = []
+    annotations = []
+
+    for dataname, datapair in data_dict.items():
+        xx = np.linspace(time[0], time[-1], 2000)
+        yy = abs(datapair[1](xx))
+
+        plots.append(
+            go.Marker(
+                x=xx,
+                y=yy,
+                mode="lines",
+                marker=go.Marker(color=datapair[2]),
+                name="${{{}}}$".format(dataname),
+            )
+        )
+
+        if dataname == "\dot{Q_c}" or dataname == "\dot{Q_h}":
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[20],
+                    y=abs(datapair[1](time[20])),
+                    text="${} = {}e^{{-{}t}}$".format(
+                        dataname,
+                        abs(round(datapair[3].args[0], 3)),
+                        round(datapair[0], 3),
+                    ),
+                    showarrow=True,
+                    ax=datapair[4],
+                    ay=datapair[5],
+                )
+            )
+
+        elif dataname == "\dot{W}":
+            annotations.append(
+                go.layout.Annotation(
+                    x=time[-20],
+                    y=datapair[1](time[-20]),
+                    text="${} = {}e^{{-{}t}} - {}e^{{-{}t}}$".format(
+                        dataname,
+                        abs(round(dQhdt.args[0], 3)),
+                        round(th_coeffs[1], 3),
+                        abs(round(dQcdt.args[0], 3)),
+                        round(tc_coeffs[1], 3),
+                    ),
+                    showarrow=True,
+                    ax=datapair[4],
+                    ay=datapair[5],
+                )
+            )
+
+    layout = go.Layout(
+        title="$\\textbf{{Graph for }} {} \\textbf{{ }} \Omega$".format(resistance),
+        plot_bgcolor="rgb(229, 229, 229)",
+        xaxis=go.layout.XAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$t (s)$",
+        ),
+        yaxis=go.layout.YAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$\dot{Q} (J/s)$",
+        ),
+        annotations=annotations,
+    )
+
+    fig = go.Figure(data=plots, layout=layout)
+
+    # Show plot on localhost
+    pio.show(fig)
+
+
+# Trapezoidal rule method
+def stage_1_sum(stage_1_time, stage_1_power, resistance):
+    total_energy = 0
+    for i in range(len(stage_1_time) - 1):
+        total_energy += ((stage_1_power[i] + stage_1_power[i + 1]) / 2) * (
+            stage_1_time[i + 1] - stage_1_time[i]
+        )
+
+    return total_energy
+
+
+# Integration method
+def stage_2_sum(time, w_coeffs, resistance):
+    t, W = symbols("t W")
+    W = w_coeffs[0] * exp(-w_coeffs[1] * (t - 61.9))
+
+    IWdt = integrate(W, (t, time[0], time[-1]))
+
+    return IWdt
+
+
+# Plot boxplot to detect anomalies
+def detect_anomalies(output_power_values):
+    sns.set(style="whitegrid")
+    ax = sns.boxplot(
+        x=output_power_values,
+        orient="h",
+        palette="Set3",
+        flierprops=dict(markerfacecolor="r", markersize=5, linestyle="none"),
+    )
+    ax.set(xlabel="Output power, $\dot{W}$")
+    plt.show()
+
+
+# Plot actual boxplot
+def plot_actual_boxplot(actual_output_power):
+    sns.set(style="whitegrid")
+    ax = sns.boxplot(
+        x=actual_output_power,
+        orient="h",
+        palette="Set3",
+        flierprops=dict(markerfacecolor="r", markersize=5, linestyle="none"),
+    )
+    ax.set(xlabel="Actual output power, $\dot{W}$")
+    plt.show()
+
+
+# Manual plotting
+def output_power_graph(resistance_values, output_power_values):
+    scatters = []
+    trendlines = []
+    annotations = []
+
+    # Use dogleg algorithm
+    popt, pcov = curve_fit(
+        output, resistance_values, output_power_values, method="dogbox"
+    )
+
+    # Compute R^2 value
+    r_squared = round(
+        metrics.r2_score(output(resistance_values, *popt), output_power_values), 3
+    )
+
+    xx = np.linspace(resistance_values[0], resistance_values[-1], 2000)
+    yy = output(xx, *popt)
+
+    scatters.append(
+        go.Marker(
+            x=resistance_values,
+            y=output_power_values,
+            mode="markers",
+            marker=go.Marker(color=list(range(10)), colorscale="Viridis"),
+            name="",
+            showlegend=False,
+        )
+    )
+
+    trendlines.append(
+        go.Marker(
+            x=xx,
+            y=yy,
+            mode="lines",
+            marker=go.Marker(color="#ff007f"),
+            name="${{\dot{W}}}$",
+        )
+    )
+
+    annotations.append(
+        go.layout.Annotation(
+            x=resistance_values[-4],
+            y=output(resistance_values[-4], *popt),
+            text="$\dot{{W}} = \\frac{{{}^2 R_L}}{{{{{} + R_L}}^2}}$".format(
+                round(popt[1], 3), round(popt[0], 3)
+            ),
+            showarrow=True,
+            ax=70,
+            ay=-5,
+        )
+    )
+
+    annotations.append(
+        go.layout.Annotation(
+            x=resistance_values[-1],
+            y=output(resistance_values[-1], *popt),
+            text="$R^2 = {}$".format(r_squared),
+            showarrow=True,
+            ax=70,
+            ay=-5,
+        )
+    )
+
+    layout = go.Layout(
+        title="$\\textbf{Output Power Graph}$",
+        plot_bgcolor="rgb(229, 229, 229)",
+        xaxis=go.layout.XAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$R_L (\Omega)$",
+        ),
+        yaxis=go.layout.YAxis(
+            zerolinecolor="rgb(255,255,255)",
+            gridcolor="rgb(255,255,255)",
+            title="$\dot{W} (W)$",
+        ),
+        annotations=annotations,
+    )
+
+    scatter_trendline_pair = list(zip(scatters, trendlines))
+    final_data = [item for sublist in scatter_trendline_pair for item in sublist]
+    fig = go.Figure(data=final_data, layout=layout)
+
+    # Show plot on localhost
+    pio.show(fig)
+
+
+# Function to graph (Ẇ, T_c and T_h) vs time
+def plot_data_graphs(resistance):
+    # Read xlsx file as dataframe
+    df = pd.read_excel(
+        "dataset.xlsx",
+        sheet_name=f"{resistance}ohm",
+        usecols="A:D",
+        na_values=" ",
+        skiprows=range(0, 67),
+    ).dropna()
+
+    stage_1 = pd.read_excel(
+        "dataset.xlsx",
+        sheet_name=f"{resistance}ohm",
+        usecols="A:D",
+        na_values=" ",
+        skiprows=[0],
+        nrows=64,
+    ).dropna()
+    # Voltage is 3.0 V
+    stage_1_power = stage_1[" Current(A)"].to_numpy() * 3
+    stage_1_time = stage_1["time(s)"].to_numpy()
+
+    # Get room temperature based on average of last two data (to minimise systematic error)
+    room_temp = (df["CH1(C)"].to_numpy()[-1] + df["CH2(C)"].to_numpy()[-1]) / 2
+
+    # Parse columns as numpy.ndarray
+    time = df["time(s)"].to_numpy()
+    t_c = df["CH1(C)"].to_numpy() - room_temp
+    t_h = df["CH2(C)"].to_numpy() - room_temp
+    power = df[" Voltage(V)"].to_numpy() ** 2 / resistance
+
+    time = time.astype(np.float64)
+    t_c = t_c.astype(np.float64)
+    t_h = t_h.astype(np.float64)
+    power = power.astype(np.float64)
+
+    # Answer question (4)
+    tc_coeffs, th_coeffs, w_coeffs = plot_exponential_graph_set(
+        time, t_c, t_h, power, resistance, room_temp
+    )
+    # Answer question (5)
+    plot_qc_qh_and_w(time, tc_coeffs, th_coeffs, resistance, room_temp)
+    # Answer question (7)
+    input = stage_1_sum(stage_1_time, stage_1_power, resistance)
+    print(f"Energy input: {input} J")
+    # Answer question (8)
+    output = stage_2_sum(time, w_coeffs, resistance)
+    print(f"Energy output: {output} J")
+    # Answer question (9)
+    efficiency = output / input * 100
+    print(f"Efficiency: {efficiency}%")
+
+
+if __name__ == "__main__":
+    for resistance in [1]:
+        plot_data_graphs(resistance)
+
+    # Manually define values for ΔT = 16.2 °C
+    resistance_values = np.array([1, 1.2, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6]).astype(
+        np.float64
+    )
+    voltage_values = np.array(
+        [
+            0.176193,
+            0.193620666666667,
+            0.247881,
+            0.3116985,
+            0.317917,
+            0.3412875,
+            0.361429,
+            0.384375,
+            0.40167875,
+        ]
+    ).astype(np.float64)
+    output_power_values = voltage_values ** 2 / resistance_values
+    output_power_values = output_power_values.astype(np.float64)
+
+    # Evidence/justification to remove anomalies
+    detect_anomalies(output_power_values)
+
+    # Remove anomaly
+    actual_resistances = np.delete(resistance_values, [3])
+    actual_voltages = np.delete(voltage_values, [3])
+
+    actual_output_power = actual_voltages ** 2 / actual_resistances
+    actual_output_power = actual_output_power.astype(np.float64)
+
+    plot_actual_boxplot(actual_output_power)
+
+    # Answer question (11)
+    output_power_graph(actual_resistances, actual_output_power)
diff --git a/images/Output.png b/hoa1/images/Output.png
similarity index 100%
rename from images/Output.png
rename to hoa1/images/Output.png
diff --git a/images/QhQc.png b/hoa1/images/QhQc.png
similarity index 100%
rename from images/QhQc.png
rename to hoa1/images/QhQc.png
diff --git a/images/ThTc.png b/hoa1/images/ThTc.png
similarity index 100%
rename from images/ThTc.png
rename to hoa1/images/ThTc.png
diff --git a/images/W.png b/hoa1/images/W.png
similarity index 100%
rename from images/W.png
rename to hoa1/images/W.png
diff --git a/images/actual_boxplot.png b/hoa1/images/actual_boxplot.png
similarity index 100%
rename from images/actual_boxplot.png
rename to hoa1/images/actual_boxplot.png
diff --git a/images/anomaly_boxplot.png b/hoa1/images/anomaly_boxplot.png
similarity index 100%
rename from images/anomaly_boxplot.png
rename to hoa1/images/anomaly_boxplot.png
diff --git a/instructions.pdf b/hoa1/instructions.pdf
similarity index 100%
rename from instructions.pdf
rename to hoa1/instructions.pdf
diff --git a/worksheet.pdf b/hoa1/worksheet.pdf
similarity index 100%
rename from worksheet.pdf
rename to hoa1/worksheet.pdf
diff --git a/hoa3/hoa3.py b/hoa3/hoa3.py
new file mode 100644
index 0000000..ab6747c
--- /dev/null
+++ b/hoa3/hoa3.py
@@ -0,0 +1,35 @@
+# 10.008 Hands-on Activity 3 Verifier Script
+# Check validity of relationship between intensity and height
+# Created by James Raphael Tiovalen (2020)
+
+# Import scientific libraries
+import numpy as np
+from scipy.optimize import curve_fit
+from sklearn import metrics
+
+
+def inverse_square(h, k):
+    return k / (h ** 2)
+
+
+if __name__ == "__main__":
+    heights = np.array([6.5, 9.0, 12.0, 17.0]).astype(np.float64)
+
+    red_intensities = np.array(
+        [0.8319219036, 0.4896690263, 0.249988924, 0.1288602701]
+    ).astype(np.float64)
+    green_intensities = np.array(
+        [0.410783097, 0.2675019527, 0.1692426359, 0.0913581608]
+    ).astype(np.float64)
+    blue_intensities = np.array(
+        [1.002523336, 0.5135877745, 0.398544113, 0.1931090032]
+    ).astype(np.float64)
+
+    for intensity in [red_intensities, green_intensities, blue_intensities]:
+        popt, pcov = curve_fit(inverse_square, heights, intensity, method="lm")
+        r_squared = round(
+            metrics.r2_score(inverse_square(heights, *popt), intensity), 3
+        )
+
+        print(f"Equation: I = {popt[0]} / h^2")
+        print(f"R^2 = {r_squared}")
diff --git a/hoa3/instructions.pdf b/hoa3/instructions.pdf
new file mode 100644
index 0000000..82e8a0b
Binary files /dev/null and b/hoa3/instructions.pdf differ
diff --git a/hoa3/worksheet.pdf b/hoa3/worksheet.pdf
new file mode 100644
index 0000000..ad1c7de
Binary files /dev/null and b/hoa3/worksheet.pdf differ