matrix rotation code

COP

                    for mid, cen, med, nor in zip(midledge, centrs, medians, normals):
                        MatrixLoc = Matrix.Translation(cen)
                        
                        MatrixRot_z = cen.rotation_difference(nor).to_matrix().to_4x4()
                        x,y,z = nor[:]
                        quota = sqrt( x**2 + y**2 )
                        if quota > 1e-12:
                            Nx = Vector(( -y/quota,   x/quota,  0 ))
                            Ny = Vector(( -x*z/quota,-y*z/quota,quota ))
                            Nz = Vector((  x,         y,        z ))
                        else:
                            Nx = Vector((1,0,0))
                            Ny = Vector((0,1,0))
                            Nz = Vector((0,0,1))
                        Y = Vector((0,1,0))
                        anglf = Y.rotation_difference(mid).angle
                        Mx = Nx*cos(anglf)-Ny*sin(anglf)
                        My = Nx*sin(anglf)+Ny*cos(anglf)
                        Mz = nor
                        MatrixRot_y = Matrix(((Mx[:]),(My[:]),(Mz[:]))).to_4x4()
                        Matrix_final = MatrixLoc * MatrixRot_y * MatrixRot_z
                        lM = [ j[:] for j in Matrix_final ]
                        '''
                        # old solution
                        loc = Matrix.Translation(cen)
                        vecz = Vector((0, 1e-6, 1))
                        q_rot0 = vecz.rotation_difference(nor).to_matrix().to_4x4()
                        q_rot2 = nor.rotation_difference(vecz).to_matrix().to_4x4()
                        vecy = Vector((1e-6, 1, 0)) * q_rot2
                        q_rot1 = vecy.rotation_difference(med).to_matrix().to_4x4()
                        M = loc*q_rot1*q_rot0
                        lM = [ j[:] for j in M ]
                        '''
                        mat_collect_.append(lM)
                    mat_collect.extend(mat_collect_)