[SOLVED] What is wrong with my MVP matrix operations?

I'm trying to program a 3D graphic engine (almost game engine) using Vulkan, following https://vulkan-tutorial.com/. I've gotten to end of the Uniform Buffer section but instead of nicely angled rotating square, my square is pivoting through the void.

I'm updating the Uniform Buffers like this:

void updateUniformBuffer(ref Application app, uint currentImage)
{
    MonoTime currentTime = MonoTime.currTime;
    float time = 0.0f;

    time = (currentTime - app.startTime).total!"msecs"() / 50.0f;
    // writeln("time: ", time);
    vec3 rotation_data = [time, 0.0f , 1.0f];
    mat4 modelMatrix;
    modelMatrix = mat4.identity().rotate_new(90.0f, rotation_data); 

    UniformBufferObject ubo = {
        model: modelMatrix,
        view: lookAt(vec3([2.0f, 02.0f, 2.0f])
            , vec3([0.0f, 0.0f, 0.0])
            , vec3([0.0f, -1.0f, -1.0f])
        ),
        proj: perspective(45.0f, app.width/cast(float)app.height, 0.1f, 10.0f),
    };

//  For debugging
//    ubo.model = mat4.identity();
//    ubo.view = mat4.identity();
//    ubo.proj = mat4.identity();

    //ubo.proj[2] = 1;
    writeln("Uniform buffer Model: \n", ubo.model);
    // writeln("Uniform buffer View: \n", ubo.view);

    ubo.proj[5] = -1 * ubo.proj[5];
    ubo.proj[15] = 1;
    // writeln("Uniform buffer Projection: \n", ubo.proj);

    void* data;
    vkMapMemory(app.device, app.uniformBuffersMemory[currentImage], 0, ubo.sizeof, 0, &data);
    memcpy(data, &ubo, ubo.sizeof);
    vkUnmapMemory(app.device, app.uniformBuffersMemory[currentImage]);
}

My Matrix Operations code looks like this:

module matrix;

import std.math;

import dlib;

vec3 normalise(ref vec3 v)
{
    float len = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
    if (len == 0.0f) 
        return v;
    v[0] /= len;
    v[1] /= len;
    v[2] /= len;
    return v;
}

vec3 crossProduct(vec3 v1, vec3 v2)
{
    vec3 res = [ 0.0f, 0.0f, 0.0f ];
    res[0] = v1[1]*v2[2] - v1[2]*v2[1];
    res[1] = v1[2]*v2[0] - v1[0]*v2[2];
    res[2] = v1[0]*v2[1] - v1[1]*v2[0];
    return res;
}

// https://www.3dgep.com/understanding-the-view-matrix/
/***
    * 
    * @param camera 
    * @param target 
    * @param up 
    * @returns 
    *
    * 1. Calculate the camera's forward vector.
    * 2. Calculate the camera's right vector.
    * 3. Calculate the camera's up vector.
    * 4. Calculate the orientation matrix.
    * 5. Calculate the translation matrix.
    * 6. Return the view matrix.
*/
mat4 lookAt(vec3 camera, vec3 target, vec3 up) 
{
    mat4 res;

    auto a = target - camera;
    normalise(a);

    auto b = crossProduct(a, up);
    normalise(b);

    auto c = crossProduct(b, a);

    res[0] = b[0]; 
    res[1] = c[0]; 
    res[2] = -a[0]; 
    res[3] = 0.0f; 
    res[4] = b[1]; 
    res[5] = c[1]; 
    res[6] = -a[1]; 
    res[7] = 0.0f; 
    res[8] = b[2]; 
    res[9] = c[2]; 
    res[10] = -a[2]; 
    res[11] = 0.0f; 
    res[12] = 0.0f;
    res[13] = 0.0f;
    res[14] = 0.0f;
    res[15] = 1.0f;    
    mat4 orientation = [
        b[0] , c[0] , a[0] , 0.0f,
        b[1] , c[1] , a[1] , 0.0f,
        b[2] , c[2] , a[2] , 0.0f,
        0.0f , 0.0f , 0.0f , 1.0f
    ];

    mat4 translation = [
        1.0f      , 0.0f      , 0.0f      , 0.0f, 
        0.0f      , 1.0f      , 0.0f      , 0.0f,  
        0.0f      , 0.0f      , 1.0f      , 0.0f, 
        -camera[0], -camera[1], -camera[2], 1.0f
    ];
    return res;
//    return orientation.multiply(translation);
}

mat4 perspective(float degrees, float ratio, float near, float far)
{
    // https://vincent-p.github.io/posts/vulkan_perspective_matrix/
    mat4 returnme = mat4.zero();
    float y = (1.0f / cos(degrees * PI / 180.0f));
    // float y = 1.0f / tan(degrees * (PI / 180.0f) * 0.5f);                        // cos(radians * PI / 180.0f);
    float x = y / ratio;
    float visible_length = far - near;

    returnme[0] = x;
    returnme[5] = -y;
    returnme[10] = -((far + near) / visible_length);  // -((far + near) / visible_length);
    returnme[11] = -1.0f; // (near * far) / (visible_length) ;
    returnme[14] = -((2.0f * near * far) / visible_length); // 1.0f;

    return returnme;
}

mat4 multiply(const mat4 m1, const mat4 m2)
{
    mat4 res;
    res[0] = m1[0]*m2[0] + m1[1]*m2[4] + m1[2]*m2[8]  + m1[3]*m2[12];
    res[1] = m1[0]*m2[1] + m1[1]*m2[5] + m1[2]*m2[9]  + m1[3]*m2[13];
    res[2] = m1[0]*m2[2] + m1[1]*m2[6] + m1[2]*m2[10] + m1[3]*m2[14];
    res[3] = m1[0]*m2[3] + m1[1]*m2[7] + m1[2]*m2[11] + m1[3]*m2[15];
    res[4] = m1[4]*m2[0] + m1[5]*m2[4] + m1[6]*m2[8]  + m1[7]*m2[12];
    res[5] = m1[4]*m2[1] + m1[5]*m2[5] + m1[6]*m2[9]  + m1[7]*m2[13];
    res[6] = m1[4]*m2[2] + m1[5]*m2[6] + m1[6]*m2[10] + m1[7]*m2[14];
    res[7] = m1[4]*m2[3] + m1[5]*m2[7] + m1[6]*m2[11] + m1[7]*m2[15];
    res[8] = m1[8]*m2[0] + m1[9]*m2[4] + m1[10]*m2[8]  + m1[11]*m2[12];
    res[9] = m1[8]*m2[1] + m1[9]*m2[5] + m1[10]*m2[9]  + m1[11]*m2[13];
    res[10] = m1[8]*m2[2] + m1[9]*m2[6] + m1[10]*m2[10] + m1[11]*m2[14];
    res[11] = m1[8]*m2[3] + m1[9]*m2[7] + m1[10]*m2[11] + m1[11]*m2[15];
    res[12] = m1[12]*m2[0] + m1[13]*m2[4] + m1[14]*m2[8]  + m1[15]*m2[12];
    res[13] = m1[12]*m2[1] + m1[13]*m2[5] + m1[14]*m2[9]  + m1[15]*m2[13];
    res[14] = m1[12]*m2[2] + m1[13]*m2[6] + m1[14]*m2[10] + m1[15]*m2[14];
    res[15] = m1[12]*m2[3] + m1[13]*m2[7] + m1[14]*m2[11] + m1[15]*m2[15];

    return res;
}

vec4 multiply(const mat4 m1, const vec3 v)
{
    vec4 res;
    res[0] = m1[0]*v[0] + m1[1]*v[1] + m1[2]*v[2];
    res[1] = m1[4]*v[1] + m1[5]*v[1] + m1[6]*v[2];
    res[2] = m1[8]*v[2] + m1[9]*v[1] + m1[10]*v[2];
    res[3] = 1.0f;

    return res;
}

mat4 multiply(const vec3 v, const mat4 m2)
{
    mat4 m1 = [ v[0], v[1], v[2], 1.0f , 
        0.0f, 0.0f, 0.0f, 0.0f,
        0.0f, 0.0f, 0.0f, 0.0f,
        0.0f, 0.0f, 0.0f, 0.0f,
    ];

    return multiply(m1, m2);
}

/**
Rotate a matrix by an angle around a vector
Parameters
    m       Input matrix multiplied by this rotation matrix.
    angle   Rotation angle expressed in radians.
    v       Rotation axis, recommended to be normalized. 
*/
mat4 rotate_new(const mat4 m, float angle, const vec3 v) 
{
  float α = (v[0] * PI) / 180.0f;
  float β = (v[1] * PI) / 180.0f; //radian(v[1]); 
  float γ = (v[2] * PI) / 180.0f; //radian(v[2]);
  // rotate around x, then y, then z
  // http://web.archive.org/web/20140515121518/http://inside.mines.edu:80/~gmurray/ArbitraryAxisRotation/ArbitraryAxisRotation.html
  mat4 mt = [
        cos(α)*cos(β), (cos(α)*sin(β)*sin(γ)) - (sin(α)*cos(γ)), (cos(α)*sin(β)*cos(γ)) + (sin(β)*sin(γ)), 0.0f,
        sin(α)*cos(β), (sin(α)*sin(β)*sin(γ)) + (cos(α)*cos(γ)), (sin(α)*sin(β)*cos(γ)) - (cos(β)*sin(γ)), 0.0f,
        -sin(β)      ,  cos(β)*sin(γ)                          ,  cos(β)*cos(γ)                          , 0.0f,
        0.0f        ,  0.0f                                   ,  0.0f,                                    1.0f 
    ]
  ;
  auto r = mt.multiply(m);
  return(r);
}

unittest {
    mat4 m = [  1, 2, 3, 4 ,
                5, 6, 7, 8 ,
                9, 10, 11, 12 ,
                13, 14, 15, 16 ];
    assert(m.multiply(m) == [ 90.0f, 100.0f, 110.0f, 120.0f, 
        202.0f, 228.0f, 254.0f, 280.0f, 
        314.0f, 356.0f, 398.0f, 440.0f, 
        426.0f, 484.0f, 542.0f, 600.0f]);

    // vec3 v = [0, 0, 0];
    // writeln(m.rotate_new(0, v));
    // assert
}

Edit 1: After altering the program as per @Thomas suggestion

    modelMatrix = mat4.identity().rotate_new(time, rotation_data); 

    UniformBufferObject ubo = {
        model: modelMatrix,
        view: lookAt(vec3([2.0f, 2.0f, 2.0f])
            , vec3([0.0f, 0.0f, 0.0])
            , vec3([0.0f, -1.0f, -1.0f])
        ),
        proj: perspective(45.0f, app.width/cast(float)app.height, 0.1f, 10.0f),
    };

the image remains but doesn't rotate. It is still at an odd angle, and clipped diagonally.