Are quaternions useful as coordinate systems for two-dimensional manifolds?

This is complementary information about my question at MSE "Are quaternions useful as coordinate systems for two-dimensional manifolds?" that I cannot add there due to the weight of the images and the length of the question. 

 This is the evolution of the first 500 steps of a test:

Same with some motion blur (basically is just a visualization test):

This is the Python code, please use and modify it freely:

 def quat():
    from sympy import mobius, factorint, totient
    from gmpy2 import is_prime, is_square
    import matplotlib.pyplot as plt
    import csv
    from random import randint
    from math import sqrt, log, cos , sin , tan, pi
  
    test_n_limit = 7000
    test_loop_limit = 500
    loa = []
    loi = []
    loj = []
    lok = []
    A=2
    I=100
    J=100
    K=4
      
    # Generating first positions of the elements
    for n in range (0,test_n_limit):
        loa.append(randint(0,A-1))
        loi.append(randint(0,I-1))
        loj.append(randint(0,J-1))
        lok.append(randint(0,K-1))
  
    for i in range(0,test_loop_limit):
        nloa = []
        nloi = []
        nloj = []
        nlok = []
        removed=[]
        for n in range (0,len(loa)):
            if n in removed:
                continue
            collided=False
            for m in range (n+1,len(loa)):
                if m in removed:
                    continue
                if loa[n]==loa[m] and loi[n]==loi[m] and loj[n]==loj[m] and lok[n]==lok[m]:
                    collided = True
                    removed.append(m)
                    # The collision produces a new point generated by a quaternion sum and the former points dissapear
                    nloi.append((loi[n]+loi[m])%I)
                    if (loi[n]+loi[m] >= I):
                        # The point goes out of the surface from the right, so it goes to the other side of the Mobius band
                        nloa.append((loa[n]+1)%A)
                        nloj.append(J-1-((loj[n]+loj[m])%J))
                    else:
                        nloa.append(loa[n])
                        nloj.append((loj[n]+loj[m])%J)
                    nlok.append((lok[n]+lok[m])%K)
                  
                elif loi[n]==loi[m] and loj[n]==loj[m] and lok[n]==lok[m]:
                    collided = True
                    removed.append(m)
                    # The collision produces a new point generated by a Hamiltonian product and the former points dissapear
                    nloa.append(((loa[n]*loa[m])-(loi[n]*loi[m])-(loj[n]*loj[m])-(lok[n]*lok[m]))%A)
                    nloi.append(((loa[n]*loi[m])+(loi[n]*loa[m])+(loj[n]*lok[m])-(lok[n]*loj[m]))%I)
                    nloj.append(((loa[n]*loj[m])-(loi[n]*lok[m])+(loj[n]*loa[m])+(lok[n]*loi[m]))%J)
                    nlok.append(((loa[n]*lok[m])+(loi[n]*loj[m])-(loj[n]*loi[m])+(lok[n]*loa[m]))%K)
                  
            if collided==False:
                # A collision did not happen and the point moves in a random walk
                incx = randint(-1,1)
                incy = randint(-1,1)
              
                nloi.append((loi[n]+incx)%I)
                if (loi[n]+incx >= I) or (loi[n]+incx < 0):
                    # The point goes out of the surface from the right or left, so it goes to the other side of the Mobius band
                    nloa.append((loa[n]+1)%A)
                    nloj.append(J-1-((loj[n]+incy)%J))
                else:
                    nloa.append(loa[n])
                    nloj.append((loj[n]+incy)%J)
                nlok.append(lok[n])
              
        print_up_i_0=[]
        print_up_j_0=[]
        print_down_i_0=[]
        print_down_j_0=[]
        print_up_i_1=[]
        print_up_j_1=[]
        print_down_i_1=[]
        print_down_j_1=[]
        print_up_i_2=[]
        print_up_j_2=[]
        print_down_i_2=[]
        print_down_j_2=[]
        print_up_i_3=[]
        print_up_j_3=[]
        print_down_i_3=[]
        print_down_j_3=[]
        for n in range (0,len(loa)):
            if loa[n]==0:
                if lok[n]==0:
                    print_up_i_0.append(loi[n])
                    print_up_j_0.append(loj[n])
                elif lok[n]==1:
                    print_up_i_1.append(loi[n])
                    print_up_j_1.append(loj[n])
                elif lok[n]==2:
                    print_up_i_2.append(loi[n])
                    print_up_j_2.append(loj[n])
                elif lok[n]==3:
                    print_up_i_3.append(loi[n])
                    print_up_j_3.append(loj[n])
              
            else:
                if lok[n]==0:
                    print_down_i_0.append(loi[n])
                    print_down_j_0.append(loj[n])
                elif lok[n]==1:
                    print_down_i_1.append(loi[n])
                    print_down_j_1.append(loj[n])
                elif lok[n]==2:
                    print_down_i_2.append(loi[n])
                    print_down_j_2.append(loj[n])
                elif lok[n]==3:
                    print_down_i_3.append(loi[n])
                    print_down_j_3.append(loj[n])
      
        # Draw current status
        print_up_i_0.append(0)
        print_up_j_0.append(0)
        print_up_i_0.append(I)
        print_up_j_0.append(J)
        plt.plot(print_up_i_0,print_up_j_0,"b,")
        plt.plot(print_down_i_0,print_down_j_0,"r,")
        plt.plot(print_up_i_1,print_up_j_1,"b.")
        plt.plot(print_down_i_1,print_down_j_1,"r.")
        plt.plot(print_up_i_2,print_up_j_2,"b*")
        plt.plot(print_down_i_2,print_down_j_2,"r*")
        plt.plot(print_up_i_3,print_up_j_3,"b+")
        plt.plot(print_down_i_3,print_down_j_3,"r+")
        #plt.show()
        #figure.set_size_inches(18, 16)
        plt.savefig("Topology_MobiusBand_quaternion_"+str(i)+".png")
        plt.clf()
      
        loa=nloa
        loi=nloi
        loj=nloj
        lok=nlok

quat()