This is the code mentioned in the appendix to my thesis. Please note that the runtime at the moment is long. If you have any suggestions as to how to optimise this, feel free to let me know.