diff --git a/hw4-pybind/python_functions.py b/hw4-pybind/python_functions.py
index 6f11d704..813316a0 100644
--- a/hw4-pybind/python_functions.py
+++ b/hw4-pybind/python_functions.py
@@ -1,34 +1,104 @@
 import numpy as np
+import subprocess
 import os
+from scipy import optimize
 
-# readPositions reads the trajectory of a given planet and makes a numpy array out of it
-# returns: 365 × 3 matrix, the columns being the three components of the planet position.
+# reads the trajectory of a given planet and makes a numpy array out of it
+#   returns: 365 × 3 matrix, the columns being the three components of the planet position.
 def readPositions(planet_name, directory):
-    out = np.zeros((365, 3))
-    for s in range(0, 365):
+    out = np.zeros((365, 3)) # initialize out matrix
+    for s in range(0, 365):  # for each time step
         filename = directory + "/step-" + str(s).zfill(5) + ".csv"
         #print(filename)
         with open(filename) as file:
             for line in file:
                 line = line.split()
                 for word in line:
                     #print(word)
                     if word == planet_name:
                         out[s, 0] = line[0]
                         out[s, 1] = line[1]
                         out[s, 2] = line[2]
                         break;
     return out
 
 # computes L2 norm of the error in position
 def computeError(positions, positions_ref):
     return np.sqrt( np.sum( (positions-positions_ref)**2) )
 
+# scales velocity of given planet in given input file
+def generateInput(scale, planet_name, input_filename, output_filename):
+    v = np.zeros(3)
 
-# Testing functions
+    # find and scale input velocity
+    content = []
+    with open(input_filename) as file:
+        for line in file:
+            content.append(line)
+            line = line.split()
+            for word in line:
+                if word == planet_name:
+                    v[0] = scale * float(line[3])
+                    v[1] = scale * float(line[4])
+                    v[2] = scale * float(line[5])
+                    break; # break for loop of word search if planet found. Go to next line.
+
+    # create new output file content
+    new_content = []
+    for line in content:
+        line = line.split()
+        if planet_name in line:
+            new_content.append( line[0] + " " + line[1] + " " + line[2] + " " + str(v[0]) + " " + str(v[1]) + " " + str(v[2]) + " " + " ".join(line[6:]) + "\n" )
+        else:
+            new_content.append( " ".join(line) + "\n" )
+
+    # write content to file
+    with open(output_filename, 'w') as new_file:
+        new_file.writelines(new_content)
+
+
+# launches the particle code on an provided input
+def launchParticles(input, nb_steps, freq):
+    ### Calling c++ executable:
+    res = subprocess.call(['./particles', str(nb_steps), str(freq), input, "planet", "1"])
+    ### Calling python function:
+    #res = subprocess.call(['python', "main.py", str(nb_steps), str(freq), input, "planet", "1"])
+
+
+# gives the error for a given scaling velocity factor of a given planet
+def runAndComputeError(scale, planet_name, input, nb_steps, freq):
+    scaled_input = input[:-4] + "_scaled.csv"
+    #print(scaled_input)
+    generateInput(scale, planet_name, input, scaled_input)
+    launchParticles(scaled_input, nb_steps, freq)
+    positions     = readPositions(planet_name, "dumps")
+    positions_ref = readPositions(planet_name, "../trajectories")
+    error = computeError(positions, positions_ref)
+    print("runAndComputeError: The L2 error for planet ", planet, " is ", error)
+    return error
+
+
+####### Testing functions #########
+
+### Note:
+### relative file paths are given here such that the script can be run from the build folder!
+
+### Exercise 5 - Question 1
 planet = "mercury"
-positions     = readPositions(planet, "dumps")
-positions_ref = readPositions(planet, "../trajectories")
+# positions     = readPositions(planet, "dumps")
+# positions_ref = readPositions(planet, "../trajectories")
+#
+### Exercise 5 - Question 2
+# error = computeError(positions, positions_ref)
+# print("The L2 error for planet ", planet, " is ", error)
+
+### Exercise 6 - Question 1
+#generateInput(2.0, planet, "../init.csv", "../init_scaled.csv")
+
+### Exercise 6 - Question 2
+#launchParticles("../init.csv", 365, 1)
 
-error = computeError(positions, positions_ref)
-print(error)
+### Exercise 6 - Question 3    
+runAndComputeError(1.0, planet, "../init.csv", 365, 1)
+approximate_scaling = optimize.fmin(runAndComputeError, 1.0, args=(planet, "../init.csv", 365, 1), xtol=0.0001, ftol=0.0001, maxiter=100000)
+print("The approximate scaling of velocity is ", approximate_scaling[0])