void ident_matrix(matrix *mat)
{
     mat->m[0][0] = 1.0f; mat->m[0][1] = 0.0f;
     mat->m[0][2] = 0.0f; mat->m[0][3] = 0.0f;

     mat->m[1][0] = 0.0f; mat->m[1][1] = 1.0f;
     mat->m[1][2] = 0.0f; mat->m[1][3] = 0.0f;

     mat->m[2][0] = 0.0f; mat->m[2][1] = 0.0f;
     mat->m[2][2] = 1.0f; mat->m[2][3] = 0.0f;

     mat->m[3][0] = 0.0f; mat->m[3][1] = 0.0f;
     mat->m[3][2] = 0.0f; mat->m[3][3] = 1.0f;
}

void rotate_matrix(matrix *m, double xtheta, double ytheta, double ztheta)
{
     float xscale, yscale, zscale;
     float sinx, cosx;
     float siny, cosy;
     float sinz, cosz;

     xscale = m->m[0][0];
     yscale = m->m[1][1];
     zscale = m->m[2][2];
     sinx = (float) sin(xtheta);
     cosx = (float) cos(xtheta);
     siny = (float) sin(ytheta);
     cosy = (float) cos(ytheta);
     sinz = (float) sin(ztheta);
     cosz = (float) cos(ztheta);

     m->m[0][0] = (float) ((cosz * cosy) * xscale);
     m->m[0][1] = (float) ((cosz * -siny * -sinx + sinz * cosx) * yscale);
     m->m[0][2] = (float) ((cosz * -siny * cosx + sinz * sinx) * zscale);

     m->m[1][0] = (float) (-sinz * cosy * xscale);
     m->m[1][1] = (float) ((-sinz * -siny * -sinx + cosz * cosx) * yscale);
     m->m[1][2] = (float) ((-sinz * -siny * cosx + cosz * sinx) * zscale);

     m->m[2][0] = (float) (siny * xscale);
     m->m[2][1] = (float) (cosy * -sinx * yscale);
     m->m[2][2] = (float) (cosy * cosx * zscale);
}

void transform_point(point3d *ptout, point3d *ptin, matrix *mat)
{
     double x, y, z;

     x = (ptin->x * mat->m[0][0]) + (ptin->y * mat->m[0][1]) +
         (ptin->z * mat->m[0][2]) + mat->m[0][3];

     y = (ptin->x * mat->m[1][0]) + (ptin->y * mat->m[1][1]) +
         (ptin->z * mat->m[1][2]) + mat->m[1][3];

     z = (ptin->x * mat->m[2][0]) + (ptin->y * mat->m[2][1])  +
   (ptin->z * mat->m[2][2]) + mat->m[2][3];

     ptout->x = (float) x;
     ptout->y = (float) y;
     ptout->z = (float) z;
}