diff --git a/README.md b/README.md
index ae71ee394944e598675482f9d5ff0df7a6d308e9..50a44cef3dbd530fc9129c2cde5596bf3cb362f5 100644
--- a/README.md
+++ b/README.md
@@ -41,3 +41,5 @@ TODO
 ---------------------
 - ?Put the header and src files in a dedicated directory
 - ?Have the compiled files be written into build
+- Use sphinx for auto generation of docs
+- Use queueing system/asynchronous task queue
diff --git a/test_david.py b/test_david.py
new file mode 100644
index 0000000000000000000000000000000000000000..59c22675547d5ea27854be53e79981dffdaccede
--- /dev/null
+++ b/test_david.py
@@ -0,0 +1,100 @@
+#!/usr/bin/python3
+import os
+
+import binary_c
+
+
+############################################################
+# Test script to run a binary using the binary_c Python
+# module.
+############################################################
+
+def run_test_binary():
+    m1 = 15.0 # Msun
+    m2 = 14.0 # Msun
+    separation = 0 # 0 = ignored, use period
+    orbital_period = 4530.0 # days
+    eccentricity = 0.0
+    metallicity = 0.02
+    max_evolution_time = 15000
+    buffer = ""
+    # argstring = "binary_c M_1 {0:g} M_2 {1:g} separation {2:g} orbital_period {3:g} eccentricity {4:g} metallicity {5:g} max_evolution_time {6:g}  ".format(m1,m2,separation,orbital_period,eccentricity,metallicity,max_evolution_time)
+    
+    argstring = f"binary_c M_1 {m1} M_2 {m2} separation {separation} orbital_period {orbital_period} eccentricity {eccentricity} metallicity {metallicity} max_evolution_time {max_evolution_time}"
+
+    output = binary_c.run_binary(argstring)
+
+    # print ("Binary_c output:\n\n")
+    print (output)
+
+
+def run_simple_loop_binary():
+    # Some simple function with a loop
+
+    m2 = 14.0 # Msun
+    separation = 0 # 0 = ignored, use period
+    orbital_period = 4530.0 # days
+    eccentricity = 0.0
+    metallicity = 0.02
+    max_evolution_time = 15000
+    buffer = ""
+
+    print ("Binary_c grid output:")
+
+    mass_1 = []
+    time = []
+    initial_m1 = []
+
+    # Set up for loop. Beware: this is only a proof of concept example. big grid should not be run in this way. rather one should use a queueing/threading system. 
+    for m1 in range(15, 20): 
+        argstring = f"binary_c M_1 {m1} M_2 {m2} separation {separation} orbital_period {orbital_period} eccentricity {eccentricity} metallicity {metallicity} max_evolution_time {max_evolution_time}"
+        
+        # Get the output string
+        output = binary_c.run_binary(argstring)
+
+        # split output on newlines
+        for line in output.split('\n'):
+            # Example of splitting off certain lines of output
+            if line.startswith('TESTLOG_BINARY_POP'):
+                value_array = line.split()[1:]
+                # print(value_array)
+
+                # Either save the results in arrays or lists
+                mass_1.append(value_array[1])
+                time.append(value_array[0])
+                initial_m1.append(value_array[2])
+
+                # write headers. can be cleaner
+                if not os.path.isfile("test_output.txt"):
+                    with open("test_output.txt", "a") as f:
+                        f.write('time M1 M1_zams\n')
+
+                # Or write them to a file with a file:
+                with open("test_output.txt", "a") as myfile:
+                    myfile.write(' '.join(value_array[:3])+'\n')
+
+    print(mass_1)
+    print(initial_m1)
+    print(time)
+
+
+# print("Current binary_c object class functions: ")
+# print(dir(binary_c))
+
+# print("Current binary_c.new_system object class functions: ")
+# print(dir(binary_c.new_system()))
+## WHat is the use of new_system?
+
+# print("Current binary_c.run_binary object class functions: ")
+# print(dir(binary_c.run_binary()))
+
+# binary_star=binary_c.new_system()
+# print(binary_star)
+# print(dir(binary_star))
+
+# Test single system
+# run_test_binary()
+
+# Test grid-like
+run_simple_loop_binary()
+