I'm learning modern OpenGL after dabbling with NeHes tutorials a few years ago. I started following the tutorials here, but I ended up with a problem when I got to the fifth tutorial in the "basics" arc. It required the C++ library GLM, which doesn't have a binding for D, which is the language I'm using for this. Having written up the code now for a simple set of matrix and vector routines, I'm testing and getting some weird results. First off, my rotating cube seems to be horrifically skewed and stretched, but only when I introduce my camera and frustum matrices. Secondly, my cube seems to start clipping when it reaches closer than 0.4 in Z screen space. I have confirmed both my matrix multiplication code, and my frustum matrix having the correct values. I'm not sure what it could be that is affecting this, so this is the code for my "idle" function which is called whenever the "update" function has returned;

auto cam = new camera;cam.translate( move_x, 0, move_z+2 );
auto model = new matrix;

// frustum matches glFrustum, see that function for arguments
auto frust = frustum(-0.5, 0.5, -0.5, 0.5, 0.1, 100 );
auto axis = (new vector( 1, 1, 1 )).normalize;
model.translate( move, 0, 0 );
model.rotate( angle, axis );
model.scale( 0.3 );

matrix transform = model * cam * frust;

glUseProgram( program );
if( trans_m != -1 )
  glUniformMatrix4fv( trans_m, 1, GL_FALSE, transform.cols.ptr );

While the drawing code just takes the vertices and colours and renders them, and the shader uses the matrix created above to transform each vertex. The cube stretches off into the distance, clips before reaching the camera, and skews when viewed from even slightly a different angle to head-on. Needless to say, I'm getting frustrated. I know there probably isn't enough here to go on, but to avoid an info-dump, if you could let me know which parts you would like to look at to work this out, I'll post them up ASAP.

Thanks in advance for any help.

EDIT: OK, so plenty of views, but no answers? Is it just that people don't have enough experience in D to help, or would it help if I add some more code? I'll include my math module below, sans any functions not being used at the moment, if that helps:

module math3d;

/+

OpenGL GLSL matrices;

Matrices are in column order, so:
    [ 0, 4,  8, 12 ]
    [ 1, 5,  9, 13 ]
    [ 2, 6, 10, 14 ]
    [ 3, 7, 11, 15 ]

Translation matrix is:
    [ 1, 0, 0, x ]
    [ 0, 1, 0, y ]
    [ 0, 0, 1, z ]
    [ 0, 0, 0, 1 ]

Rotation( angle, axis ) matrix is:
    [ (x^2)*(1-cos(angle))+cos(angle), x*y*(1-cos(angle))-z*sin(angle), x*z*(1-cos(angle))+y*sin(angle), 0 ]
    [ x*y*(1-cos(angle))+z*sin(angle), (y^2)*(1-cos(angle))+cos(angle), y*z*(1-cos(angle))-x*sin(angle), 0 ]
    [ x*z*(1-cos(angle))-y*sin(angle), y*z*(1-cos(angle))+x*sin(angle), (z^2)*(1-cos(angle))+cos(angle), 0 ]
    [                 0,                                    0,                                0,                 1 ]

Scale matrix is:
    [ x, 0, 0, 0 ]
    [ 0, y, 0, 0 ]
    [ 0, 0, z, 0 ]
    [ 0, 0, 0, 1 ]

+/

public import std.math;

import std.stdio;

class vector{
    float[4] array;
    
    this( float a, float b, float c, float d = 1.0f ){
        array = [ a, b, c, d ];
    }
    
    this(){
        this( 0.0, 0.0, 0.0 );
    }
    
    this( float[4] r ){
        array = r;
    }
    
    this( float[] r ){
        if( r.length >=4 )
            array = r;
        else{
            array[0..r.length] = r;
            array[r.length..$] = 0;
            array[3] = 1;
        }
    }
    
    @property float x(){ return array[0]; }
    @property float y(){ return array[1]; }
    @property float z(){ return array[2]; }
    @property float w(){ return array[3]; }
    
    @property float x( float c ){ return array[0] = c; }
    @property float y( float c ){ return array[1] = c; }
    @property float z( float c ){ return array[2] = c; }
    @property float w( float c ){ return array[3] = c; }
    
    vector opBinary( string op )( vector v ){
        return mixin( "new vector( array"~op~"v.array )" );
    }
    
    vector opBinary( string op )( matrix m ){
        vector ret = new vector;
        static if( op == "*" ){
            ret.x = x*m.cols[0][0] + x*m.cols[0][1] + x*m.cols[0][2] + x*m.cols[0][3];
            ret.y = y*m.cols[0][0] + y*m.cols[0][1] + y*m.cols[0][2] + y*m.cols[0][3];
            ret.z = z*m.cols[0][0] + z*m.cols[0][1] + z*m.cols[0][2] + z*m.cols[0][3];
            ret.w = w*m.cols[0][0] + w*m.cols[0][1] + w*m.cols[0][2] + w*m.cols[0][3];
        }
    }
    
    vector normalize(){
        import std.math;
        
        auto sqr = (x*x)+(y*y)+(z*z);
        if( sqr != 1.0f ){
            auto srt = sqrt( sqr );
            array[0] /= srt;
            array[1] /= srt;
            array[2] /= srt;
            array[3]  = 1;
        }
        
        return this;
    }
}

/// simple matrix class. Matrix values are stored in the array in OpenGL( i.e.; column-major ) order.
class matrix{
    float[16] cols;
    private matrix previous;
    
    /// allows us to treat the matrix as though it were the array of 16 floats
    alias cols this;
    
    ////////////////
    // CONSTRUCTORS
    ////////////////
    this(){
        cols =
            [ 1.0, 0.0, 0.0, 0.0,
              0.0, 1.0, 0.0, 0.0,
              0.0, 0.0, 1.0, 0.0,
              0.0, 0.0, 0.0, 1.0 ];
        
    }
    
    this( float[16] vals ){
        cols = vals.dup;
    }
    
    this( float[4][4] vals ){
        cols = vals[0]~vals[1]~vals[2]~vals[3];
    }
    
    this( matrix m ){
        cols = m.cols.dup;
    }
    
    /// returns an identity matrix.
    /// Not difficult, as the identity matrix is the default for the matrix class
    static matrix identity(){
        return new matrix();
    }
    
    /// operator overloading. Currently only overloads the '*' operator.
    matrix opBinary( string op )( matrix m ){
        matrix r = new matrix;
        
        /++
        
        At the moment, if the matrices are layed out in the following way;
        
                               [ t00 t04 t08 t12 ]
                               [ t01 t05 t09 t13 ]
                               [ t02 t06 t10 t14 ]
                               [ t03 t07 t11 t15 ]
        
        [ m00 m04 m08 m12 ]    [ r00 r04 r08 r12 ]
        [ m01 m05 m09 m13 ]    [ r01 r05 r09 r13 ]
        [ m02 m06 m10 m14 ]    [ r02 r06 r10 r14 ]
        [ m03 m07 m11 m15 ]    [ r03 r07 r11 r15 ]
        
        where rXX is the resulting value, tXX is the "this" matrix, and mXX is
        the passed-in matrix. If the matrices are messing up, then I may need to
        switch the m and t matrices.
        
        ++/
        static if( op == "*" ){
            auto t = cols;
            r[0]  = m[0]*t[0]  + m[4]*t[1]  + m[8]* t[2]  + m[12]*t[3];
            r[1]  = m[1]*t[0]  + m[5]*t[1]  + m[9]* t[2]  + m[13]*t[3];
            r[2]  = m[2]*t[0]  + m[6]*t[1]  + m[10]*t[2]  + m[14]*t[3];
            r[3]  = m[3]*t[0]  + m[7]*t[1]  + m[11]*t[2]  + m[15]*t[3];
            
            r[4]  = m[0]*t[4]  + m[4]*t[5]  + m[8]* t[6]  + m[12]*t[7];
            r[5]  = m[1]*t[4]  + m[5]*t[5]  + m[9]* t[6]  + m[13]*t[7];
            r[6]  = m[2]*t[4]  + m[6]*t[5]  + m[10]*t[6]  + m[14]*t[7];
            r[7]  = m[3]*t[4]  + m[7]*t[5]  + m[11]*t[6]  + m[15]*t[7];
            
            r[8]  = m[0]*t[8]  + m[4]*t[9]  + m[8]* t[10] + m[12]*t[11];
            r[9]  = m[1]*t[8]  + m[5]*t[9]  + m[9]* t[10] + m[13]*t[11];
            r[10] = m[2]*t[8]  + m[6]*t[9]  + m[10]*t[10] + m[14]*t[11];
            r[11] = m[3]*t[8]  + m[7]*t[9]  + m[11]*t[10] + m[15]*t[11];
            
            r[12] = m[0]*t[12] + m[4]*t[13] + m[8]* t[14] + m[12]*t[15];
            r[13] = m[1]*t[12] + m[5]*t[13] + m[9]* t[14] + m[13]*t[15];
            r[14] = m[2]*t[12] + m[6]*t[13] + m[10]*t[14] + m[14]*t[15];
            r[15] = m[3]*t[12] + m[7]*t[13] + m[11]*t[14] + m[15]*t[15];
            
        }
        else static assert( 0, "Unsupported binary op in math3d.matrix: "~op );
        
        return r;
    }
    
    /// utility method to perform a translation on a matrix
    matrix translate( vector v ){
        return translate( v.x, v.y, v.z );
    }
    
    /// utility method to perform a translation on a matrix
    matrix translate( float x, float y, float z ){
        
        cols[0] += cols[3]*x;    cols[4] += cols[7]*x;    cols[8] += cols[11]*x;    cols[12] += cols[15]*x;
        cols[1] += cols[3]*y;    cols[5] += cols[7]*y;    cols[9] += cols[11]*y;    cols[13] += cols[15]*y;
        cols[2] += cols[3]*z;    cols[6] += cols[7]*z;    cols[10] += cols[11]*z;    cols[14] += cols[15]*z;
        
        version(test_matrix) print(this);
        return this;
    }
    
    /// utility method to perform a rotation on a matrix
    matrix rotate( float angle, vector axis ){
        version(test_matrix) writefln("entered rotate(float, vector(%s))", axis);
        return rotate( angle, axis.x, axis.y, axis.z );
        version(test_matrix) writeln("leaving rotate(float, vector)");
    }
    
    /// utility method to perform a rotation on a matrix
    matrix rotate( float angle, float x, float y, float z ){
        version(test_matrix) writefln( "entered rotate(%s, %s, %s, %s)", angle, x, y, z );
        matrix m = new matrix;
        version(test_matrix) writeln( "\nmatrix m before construction of rotation" );
        
        float c = cos(angle*(PI/180.0));
        float s = sin(angle*(PI/180.0));
        version(test_matrix){
            writefln( "c: %s\ns: %s\n", c, s );
            writefln( "cos(angle*(PI/180.0)): %s\nsin(angle*(PI/180.0)): %s\n", cos(angle*(PI/180.0)), sin(angle*(PI/180.0)) );
        }
        
        auto xx = x*x,            xy = x*y,
             xz = x*z,            yy = y*y,
             yz = y*z,            zz = z*z;
        
        m.cols[0] = xx * (1 - c) + c;
        m.cols[1] = xy * (1 - c) + z * s;
        m.cols[2] = xz * (1 - c) - y * s;
        m.cols[3]= 0;
        m.cols[4] = xy * (1 - c) - z * s;
        m.cols[5] = yy * (1 - c) + c;
        m.cols[6] = yz * (1 - c) + x * s;
        m.cols[7]= 0;
        m.cols[8] = xz * (1 - c) + y * s;
        m.cols[9] = yz * (1 - c) - x * s;
        m.cols[10]= zz * (1 - c) + c;
        m.cols[11]= 0;
        m.cols[12] = 0;
        m.cols[13] = 0;
        m.cols[14]= 0;
        m.cols[15]= 1;
        
        version(test_matrix){
            writeln( "\nmatrix m after rotation created:" );
            print(m);
        }
        
        cols = ( m*this ).cols;
        
        version(test_matrix){
            writeln( "\nfinal matrix after rotation applied:" );
            print(this);
        }
        
        version(test_matrix) writeln("leaving rotate(float, float, float, float)");
        return this;
    }
    
    matrix rotateX( float angle ){
        version(test_matrix) writeln("entered rotateX");
        return rotate( angle, 1.0, 0.0, 0.0 );
    }
    matrix rotateY( float angle ){
        version(test_matrix) writeln("entered rotateY");
        return rotate( angle, 0.0, 1.0, 0.0 );
    }
    matrix rotateZ( float angle ){
        version(test_matrix) writeln("entered rotateX");
        return rotate( angle, 0.0, 0.0, 1.0 );
    }
    
    /// utility method to scale the matrix
    matrix scale( vector v ){
        return scale( v.x, v.y, v.z );
    }
    
    matrix scale( float f ){
        return scale( f, f, f );
    }
    
    matrix scale( float x, float y, float z ){
        cols[0] = cols[0]*x;   cols[1] = cols[1]*x;   cols[2] = cols[2]*x;   cols[3] = cols[3]*x;
        cols[4] = cols[4]*y;   cols[5] = cols[5]*y;   cols[6] = cols[6]*y;   cols[7] = cols[7]*y;
        cols[8] = cols[8]*z;   cols[9] = cols[9]*z;   cols[10]= cols[10]*z;  cols[11]= cols[11]*z;
        
        return this;
    }
}

/// replacement for gluPerspective. Matches the parameters and results as accurately as the floats allow
matrix perspective( float fovInDeg, float aspectRatio, float znear, float zfar ){
    float ymax, xmax;
    ymax = znear * tan( fovInDeg * PI / 360.0 );
    xmax = ymax * aspectRatio;
    return frustum( -xmax, xmax, -ymax, ymax, znear, zfar );
}

/// replacement for glFrustum. Matches the parameters and results as accurately as the floats allow
matrix frustum( float left, float right, float bottom, float top, float znear, float zfar ){
    float t1, t2, t3, t4;
    t1 = 2.0 * znear;
    t2 = right - left;
    t3 = top - bottom;
    t4 = zfar - znear;
    
    return new matrix(
        [
            t1/t2, 0, 0, 0,
            0, t1/t3, 0, 0,
            (right + left)/t2, (top + bottom)/t3, -(zfar + znear)/t4, -1,
            0, 0, -(t1*zfar)/t4, 0
        ]
    );
}

I hope this helps, and I'd appreciate any help you guys can give me