diff --git a/kiln-tuner.py b/kiln-tuner.py
new file mode 100644
index 0000000..5e69f26
--- /dev/null
+++ b/kiln-tuner.py
@@ -0,0 +1,110 @@
+import csv
+import argparse
+
+
+# Using the method from https://www.ias.ac.in/article/fulltext/reso/025/10/1385-1397 or https://www.youtube.com/watch?v=nvAQHSe-Ax4
+
+
+def line(a, b, x):
+    return a * x + b
+
+
+def invline(a, b, y):
+    return (y - b) / a
+
+
+def plot(xdata, ydata,
+         tangent_min, tangent_max, tangent_slope, tangent_offset,
+         lower_crossing_x, upper_crossing_x):
+    from matplotlib import pyplot
+
+    minx = min(xdata)
+    maxx = max(xdata)
+    miny = min(ydata)
+    maxy = max(ydata)
+
+    pyplot.scatter(xdata, ydata)
+
+    pyplot.plot([minx, maxx], [miny, miny], '--', color='purple')
+    pyplot.plot([minx, maxx], [maxy, maxy], '--', color='purple')
+
+    pyplot.plot(tangent_min[0], tangent_min[1], 'v', color='red')
+    pyplot.plot(tangent_max[0], tangent_max[1], 'v', color='red')
+    pyplot.plot([minx, maxx], [line(tangent_slope, tangent_offset, minx), line(tangent_slope, tangent_offset, maxx)], '--', color='red')
+
+    pyplot.plot([lower_crossing_x, lower_crossing_x], [miny, maxy], '--', color='black')
+    pyplot.plot([upper_crossing_x, upper_crossing_x], [miny, maxy], '--', color='black')
+
+    pyplot.show()
+
+
+def calculate(filename, tangentdivisor, showplot):
+    # parse the csv file
+    xdata = []
+    ydata = []
+    filemintime = None
+    with open(filename) as f:
+        for row in csv.DictReader(f):
+            try:
+                time = float(row['pid_time'])
+                temp = float(row['pid_ispoint'])
+                if filemintime is None:
+                    filemintime = time
+
+                xdata.append(time - filemintime)
+                ydata.append(temp)
+            except ValueError:
+                continue  # just ignore bad values!
+
+    # gather points for tangent line
+    miny = min(ydata)
+    maxy = max(ydata)
+    midy = (maxy + miny) / 2
+    yoffset = int((maxy - miny) / tangentdivisor)
+    tangent_min = tangent_max = None
+    for i in range(0, len(xdata)):
+        rowx = xdata[i]
+        rowy = ydata[i]
+
+        if rowy >= (midy - yoffset) and tangent_min is None:
+            tangent_min = (rowx, rowy)
+        elif rowy >= (midy + yoffset) and tangent_max is None:
+            tangent_max = (rowx, rowy)
+
+    # calculate tangent line to the main temperature curve
+    tangent_slope = (tangent_max[1] - tangent_min[1]) / (tangent_max[0] - tangent_min[0])
+    tangent_offset = tangent_min[1] - line(tangent_slope, 0, tangent_min[0])
+
+    # determine the point at which the tangent line crosses the min/max temperaturess
+    lower_crossing_x = invline(tangent_slope, tangent_offset, miny)
+    upper_crossing_x = invline(tangent_slope, tangent_offset, maxy)
+
+    # compute parameters
+    L = lower_crossing_x - min(xdata)
+    T = upper_crossing_x - lower_crossing_x
+    Kp = 1.2 * (T / L)
+    Ti = 2 * L
+    Td = 0.5 * L
+    Ki = Kp / Ti
+    Kd = Kp * Td
+
+    # outut to the user
+    print(Kp, 1 / Ki, Kd)
+
+    if showplot:
+        plot(xdata, ydata,
+             tangent_min, tangent_max, tangent_slope, tangent_offset,
+             lower_crossing_x, upper_crossing_x)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Perform Ziegler-Nichols PID tuning')
+    parser.add_argument('csvfile', type=str, help="The CSV file to read from. Must contain two columns called pid_time (time in seconds) and pid_ispoint (observed temperature)")
+    parser.add_argument('--showplot', action='store_true', help="If set, also plot results (requires pyplot to be pip installed)")
+    parser.add_argument('--tangentdivisor', type=float, default=4, help="Adjust the tangent calculation to fit better. Must be >= 2.")
+    args = parser.parse_args()
+
+    if args.tangentdivisor < 2:
+        raise ValueError("tangentdivisor must be >= 2")
+
+    calculate(args.csvfile, args.tangentdivisor, args.showplot)