khora_core/math/
matrix.rs

1// Copyright 2025 eraflo
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7//     http://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15//! Defines the `Mat3` and `Mat4` types and associated operations.
16
17use super::{Quaternion, Vec2, Vec3, Vec4, EPSILON};
18use std::ops::{Index, IndexMut, Mul};
19
20// --- Mat3 ---
21
22/// A 3x3 column-major matrix, typically used for 2D affine transformations (scale, rotation).
23///
24/// While it can represent any 3x3 matrix, its primary role in a 3D engine is often
25/// as the upper-left rotation and scale part of a `Mat4`.
26#[derive(Debug, Clone, Copy, PartialEq)]
27#[repr(C)]
28pub struct Mat3 {
29    /// The columns of the matrix. `cols[0]` is the first column, and so on.
30    pub cols: [Vec3; 3],
31}
32
33impl Mat3 {
34    /// The 3x3 identity matrix.
35    pub const IDENTITY: Self = Self {
36        cols: [Vec3::X, Vec3::Y, Vec3::Z],
37    };
38
39    /// A 3x3 matrix with all elements set to 0.
40    pub const ZERO: Self = Self {
41        cols: [Vec3::ZERO; 3],
42    };
43
44    /// Creates a new matrix from three column vectors.
45    #[inline]
46    pub fn from_cols(c0: Vec3, c1: Vec3, c2: Vec3) -> Self {
47        Self { cols: [c0, c1, c2] }
48    }
49
50    /// Returns a row of the matrix as a `Vec3`.
51    #[allow(dead_code)]
52    #[inline]
53    fn get_row(&self, index: usize) -> Vec3 {
54        Vec3 {
55            x: self.cols[0].get(index),
56            y: self.cols[1].get(index),
57            z: self.cols[2].get(index),
58        }
59    }
60
61    /// Creates a 2D scaling matrix.
62    ///
63    /// The Z-axis scale is set to 1.0, making it a no-op in that dimension.
64    #[inline]
65    pub fn from_scale_vec2(scale: Vec2) -> Self {
66        Self::from_scale(Vec3::new(scale.x, scale.y, 1.0))
67    }
68
69    /// Creates a 3D scaling matrix.
70    #[inline]
71    pub fn from_scale(scale: Vec3) -> Self {
72        Self {
73            cols: [
74                Vec3::new(scale.x, 0.0, 0.0),
75                Vec3::new(0.0, scale.y, 0.0),
76                Vec3::new(0.0, 0.0, scale.z),
77            ],
78        }
79    }
80
81    /// Creates a matrix for a rotation around the X-axis.
82    ///
83    /// # Arguments
84    ///
85    /// * `angle_radians`: The angle of rotation in radians.
86    #[inline]
87    pub fn from_rotation_x(angle_radians: f32) -> Self {
88        let (s, c) = angle_radians.sin_cos();
89        Self {
90            cols: [
91                Vec3::new(1.0, 0.0, 0.0),
92                Vec3::new(0.0, c, s),
93                Vec3::new(0.0, -s, c),
94            ],
95        }
96    }
97
98    /// Creates a matrix for a right-handed rotation around the Y-axis.
99    ///
100    /// # Arguments
101    ///
102    /// * `angle_radians`: The angle of rotation in radians.
103    #[inline]
104    pub fn from_rotation_y(angle_radians: f32) -> Self {
105        let (s, c) = angle_radians.sin_cos();
106        Self {
107            cols: [
108                Vec3::new(c, 0.0, -s),
109                Vec3::new(0.0, 1.0, 0.0),
110                Vec3::new(s, 0.0, c),
111            ],
112        }
113    }
114
115    /// Creates a matrix for a rotation around the Z-axis.
116    ///
117    /// # Arguments
118    ///
119    /// * `angle_radians`: The angle of rotation in radians.
120    #[inline]
121    pub fn from_rotation_z(angle_radians: f32) -> Self {
122        let (s, c) = angle_radians.sin_cos();
123        Self {
124            cols: [
125                Vec3::new(c, s, 0.0),
126                Vec3::new(-s, c, 0.0),
127                Vec3::new(0.0, 0.0, 1.0),
128            ],
129        }
130    }
131
132    /// Creates a rotation matrix from a normalized axis and an angle.
133    ///
134    /// # Arguments
135    ///
136    /// * `axis`: The axis of rotation. Must be a unit vector.
137    /// * `angle_radians`: The angle of rotation in radians.
138    #[inline]
139    pub fn from_axis_angle(axis: Vec3, angle_radians: f32) -> Self {
140        let (s, c) = angle_radians.sin_cos();
141        let t = 1.0 - c;
142        let x = axis.x;
143        let y = axis.y;
144        let z = axis.z;
145        Self {
146            cols: [
147                Vec3::new(t * x * x + c, t * x * y + s * z, t * x * z - s * y),
148                Vec3::new(t * y * x - s * z, t * y * y + c, t * y * z + s * x),
149                Vec3::new(t * z * x + s * y, t * z * y - s * x, t * z * z + c),
150            ],
151        }
152    }
153
154    /// Creates a rotation matrix from a quaternion.
155    /// The quaternion is normalized before conversion to ensure a valid rotation matrix.
156    #[inline]
157    pub fn from_quat(q: Quaternion) -> Self {
158        let q = q.normalize();
159        let x = q.x;
160        let y = q.y;
161        let z = q.z;
162        let w = q.w;
163        let x2 = x + x;
164        let y2 = y + y;
165        let z2 = z + z;
166        let xx = x * x2;
167        let xy = x * y2;
168        let xz = x * z2;
169        let yy = y * y2;
170        let yz = y * z2;
171        let zz = z * z2;
172        let wx = w * x2;
173        let wy = w * y2;
174        let wz = w * z2;
175
176        Self::from_cols(
177            Vec3::new(1.0 - (yy + zz), xy + wz, xz - wy),
178            Vec3::new(xy - wz, 1.0 - (xx + zz), yz + wx),
179            Vec3::new(xz + wy, yz - wx, 1.0 - (xx + yy)),
180        )
181    }
182
183    /// Creates a `Mat3` from the upper-left 3x3 corner of a [`Mat4`].
184    /// This effectively extracts the rotation and scale components, discarding translation.
185    #[inline]
186    pub fn from_mat4(m4: &Mat4) -> Self {
187        Self::from_cols(
188            m4.cols[0].truncate(),
189            m4.cols[1].truncate(),
190            m4.cols[2].truncate(),
191        )
192    }
193
194    /// Computes the determinant of the matrix.
195    ///
196    /// The determinant is a scalar value indicating the volume scaling factor of the
197    /// linear transformation. A determinant of 0 means the matrix is not invertible.
198    #[inline]
199    pub fn determinant(&self) -> f32 {
200        let c0 = self.cols[0];
201        let c1 = self.cols[1];
202        let c2 = self.cols[2];
203        c0.x * (c1.y * c2.z - c2.y * c1.z) - c1.x * (c0.y * c2.z - c2.y * c0.z)
204            + c2.x * (c0.y * c1.z - c1.y * c0.z)
205    }
206
207    /// Returns the transpose of the matrix, where rows and columns are swapped.
208    #[inline]
209    pub fn transpose(&self) -> Self {
210        Self::from_cols(
211            Vec3::new(self.cols[0].x, self.cols[1].x, self.cols[2].x),
212            Vec3::new(self.cols[0].y, self.cols[1].y, self.cols[2].y),
213            Vec3::new(self.cols[0].z, self.cols[1].z, self.cols[2].z),
214        )
215    }
216
217    /// Computes the inverse of the matrix.
218    ///
219    /// If the matrix is not invertible (i.e., its determinant is close to zero),
220    /// this method returns `None`.
221    pub fn inverse(&self) -> Option<Self> {
222        let c0 = self.cols[0];
223        let c1 = self.cols[1];
224        let c2 = self.cols[2];
225        let m00 = c1.y * c2.z - c2.y * c1.z;
226        let m10 = c2.y * c0.z - c0.y * c2.z;
227        let m20 = c0.y * c1.z - c1.y * c0.z;
228        let det = c0.x * m00 + c1.x * m10 + c2.x * m20;
229
230        if det.abs() < EPSILON {
231            return None;
232        }
233
234        let inv_det = 1.0 / det;
235        let m01 = c2.x * c1.z - c1.x * c2.z;
236        let m11 = c0.x * c2.z - c2.x * c0.z;
237        let m21 = c1.x * c0.z - c0.x * c1.z;
238        let m02 = c1.x * c2.y - c2.x * c1.y;
239        let m12 = c2.x * c0.y - c0.x * c2.y;
240        let m22 = c0.x * c1.y - c1.x * c0.y;
241
242        Some(Self::from_cols(
243            Vec3::new(m00, m10, m20) * inv_det,
244            Vec3::new(m01, m11, m21) * inv_det,
245            Vec3::new(m02, m12, m22) * inv_det,
246        ))
247    }
248
249    /// Converts this `Mat3` into a [`Mat4`], preserving its values in the upper-left corner.
250    /// The new fourth column and row are set to `(0, 0, 0, 1)`.
251    #[inline]
252    pub fn to_mat4(&self) -> Mat4 {
253        Mat4::from_cols(
254            Vec4::from_vec3(self.cols[0], 0.0),
255            Vec4::from_vec3(self.cols[1], 0.0),
256            Vec4::from_vec3(self.cols[2], 0.0),
257            Vec4::W,
258        )
259    }
260}
261
262// --- Operator Overloads ---
263
264impl Default for Mat3 {
265    /// Returns the 3x3 identity matrix.
266    #[inline]
267    fn default() -> Self {
268        Self::IDENTITY
269    }
270}
271
272impl Mul<Mat3> for Mat3 {
273    type Output = Self;
274    /// Multiplies this matrix by another `Mat3`.
275    #[inline]
276    fn mul(self, rhs: Mat3) -> Self::Output {
277        Self::from_cols(self * rhs.cols[0], self * rhs.cols[1], self * rhs.cols[2])
278    }
279}
280
281impl Mul<Vec3> for Mat3 {
282    type Output = Vec3;
283    /// Transforms a `Vec3` by this matrix.
284    #[inline]
285    fn mul(self, v: Vec3) -> Self::Output {
286        self.cols[0] * v.x + self.cols[1] * v.y + self.cols[2] * v.z
287    }
288}
289
290impl Index<usize> for Mat3 {
291    type Output = Vec3;
292    /// Allows accessing a matrix column by index.
293    #[inline]
294    fn index(&self, index: usize) -> &Self::Output {
295        &self.cols[index]
296    }
297}
298
299impl IndexMut<usize> for Mat3 {
300    /// Allows mutably accessing a matrix column by index.
301    #[inline]
302    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
303        &mut self.cols[index]
304    }
305}
306
307// --- Mat4 ---
308
309/// A 4x4 column-major matrix, used for 3D affine transformations.
310///
311/// This is the primary type for representing transformations (translation, rotation,
312/// scale) in 3D space. It is also used for camera view and projection matrices.
313/// The memory layout is column-major, which is compatible with modern graphics APIs
314/// like Vulkan, Metal, and DirectX.
315#[derive(Debug, Clone, Copy, PartialEq)]
316#[repr(C)]
317pub struct Mat4 {
318    /// The columns of the matrix. `cols[0]` is the first column, and so on.
319    pub cols: [Vec4; 4],
320}
321
322impl Mat4 {
323    /// The 4x4 identity matrix.
324    pub const IDENTITY: Self = Self {
325        cols: [Vec4::X, Vec4::Y, Vec4::Z, Vec4::W],
326    };
327
328    /// A 4x4 matrix with all elements set to 0.
329    pub const ZERO: Self = Self {
330        cols: [Vec4::ZERO; 4],
331    };
332
333    /// Creates a new matrix from four column vectors.
334    #[inline]
335    pub fn from_cols(c0: Vec4, c1: Vec4, c2: Vec4, c3: Vec4) -> Self {
336        Self {
337            cols: [c0, c1, c2, c3],
338        }
339    }
340
341    /// Returns a row of the matrix as a `Vec4`.
342    #[inline]
343    fn get_row(&self, index: usize) -> Vec4 {
344        Vec4 {
345            x: self.cols[0].get(index),
346            y: self.cols[1].get(index),
347            z: self.cols[2].get(index),
348            w: self.cols[3].get(index),
349        }
350    }
351
352    /// Creates a translation matrix.
353    ///
354    /// # Arguments
355    ///
356    /// * `v`: The translation vector to apply.
357    #[inline]
358    pub fn from_translation(v: Vec3) -> Self {
359        Self {
360            cols: [
361                Vec4::new(1.0, 0.0, 0.0, 0.0),
362                Vec4::new(0.0, 1.0, 0.0, 0.0),
363                Vec4::new(0.0, 0.0, 1.0, 0.0),
364                Vec4::new(v.x, v.y, v.z, 1.0),
365            ],
366        }
367    }
368
369    /// Creates a non-uniform scaling matrix.
370    #[inline]
371    pub fn from_scale(scale: Vec3) -> Self {
372        Self {
373            cols: [
374                Vec4::new(scale.x, 0.0, 0.0, 0.0),
375                Vec4::new(0.0, scale.y, 0.0, 0.0),
376                Vec4::new(0.0, 0.0, scale.z, 0.0),
377                Vec4::new(0.0, 0.0, 0.0, 1.0),
378            ],
379        }
380    }
381
382    /// Creates a matrix for a rotation around the X-axis.
383    ///
384    /// # Arguments
385    ///
386    /// * `angle`: The angle of rotation in radians.
387    #[inline]
388    pub fn from_rotation_x(angle: f32) -> Self {
389        let c = angle.cos();
390        let s = angle.sin();
391        Self {
392            cols: [
393                Vec4::new(1.0, 0.0, 0.0, 0.0),
394                Vec4::new(0.0, c, s, 0.0),
395                Vec4::new(0.0, -s, c, 0.0),
396                Vec4::new(0.0, 0.0, 0.0, 1.0),
397            ],
398        }
399    }
400
401    /// Creates a matrix for a right-handed rotation around the Y-axis.
402    ///
403    /// # Arguments
404    ///
405    /// * `angle`: The angle of rotation in radians.
406    #[inline]
407    pub fn from_rotation_y(angle: f32) -> Self {
408        let c = angle.cos();
409        let s = angle.sin();
410        Self {
411            cols: [
412                Vec4::new(c, 0.0, -s, 0.0),
413                Vec4::new(0.0, 1.0, 0.0, 0.0),
414                Vec4::new(s, 0.0, c, 0.0),
415                Vec4::new(0.0, 0.0, 0.0, 1.0),
416            ],
417        }
418    }
419
420    /// Creates a matrix for a rotation around the Z-axis.
421    ///
422    /// # Arguments
423    ///
424    /// * `angle`: The angle of rotation in radians.
425    #[inline]
426    pub fn from_rotation_z(angle: f32) -> Self {
427        let c = angle.cos();
428        let s = angle.sin();
429        Self {
430            cols: [
431                Vec4::new(c, s, 0.0, 0.0),
432                Vec4::new(-s, c, 0.0, 0.0),
433                Vec4::new(0.0, 0.0, 1.0, 0.0),
434                Vec4::new(0.0, 0.0, 0.0, 1.0),
435            ],
436        }
437    }
438
439    /// Creates a rotation matrix from a normalized axis and an angle.
440    ///
441    /// # Arguments
442    ///
443    /// * `axis`: The axis of rotation. Must be a unit vector.
444    /// * `angle`: The angle of rotation in radians.
445    #[inline]
446    pub fn from_axis_angle(axis: Vec3, angle: f32) -> Self {
447        let c = angle.cos();
448        let s = angle.sin();
449        let t = 1.0 - c;
450        let x = axis.x;
451        let y = axis.y;
452        let z = axis.z;
453
454        Self {
455            cols: [
456                Vec4::new(t * x * x + c, t * x * y - s * z, t * x * z + s * y, 0.0),
457                Vec4::new(t * y * x + s * z, t * y * y + c, t * y * z - s * x, 0.0),
458                Vec4::new(t * z * x - s * y, t * z * y + s * x, t * z * z + c, 0.0),
459                Vec4::new(0.0, 0.0, 0.0, 1.0),
460            ],
461        }
462    }
463
464    /// Creates a rotation matrix from a quaternion.
465    #[inline]
466    pub fn from_quat(q: Quaternion) -> Self {
467        let x = q.x;
468        let y = q.y;
469        let z = q.z;
470        let w = q.w;
471        let x2 = x + x;
472        let y2 = y + y;
473        let z2 = z + z;
474        let xx = x * x2;
475        let xy = x * y2;
476        let xz = x * z2;
477        let yy = y * y2;
478        let yz = y * z2;
479        let zz = z * z2;
480        let wx = w * x2;
481        let wy = w * y2;
482        let wz = w * z2;
483
484        Self::from_cols(
485            Vec4::new(1.0 - (yy + zz), xy + wz, xz - wy, 0.0),
486            Vec4::new(xy - wz, 1.0 - (xx + zz), yz + wx, 0.0),
487            Vec4::new(xz + wy, yz - wx, 1.0 - (xx + yy), 0.0),
488            Vec4::W,
489        )
490    }
491
492    /// Creates a right-handed perspective projection matrix with a [0, 1] depth range (ZO).
493    ///
494    /// # Arguments
495    ///
496    /// * `fov_y_radians`: Vertical field of view in radians.
497    /// * `aspect_ratio`: Width divided by height of the viewport.
498    /// * `z_near`: Distance to the near clipping plane (must be positive).
499    /// * `z_far`: Distance to the far clipping plane (must be positive and > `z_near`).
500    #[inline]
501    pub fn perspective_rh_zo(
502        fov_y_radians: f32,
503        aspect_ratio: f32,
504        z_near: f32,
505        z_far: f32,
506    ) -> Self {
507        assert!(z_near > 0.0 && z_far > z_near);
508        let tan_half_fovy = (fov_y_radians / 2.0).tan();
509        let f = 1.0 / tan_half_fovy;
510        let aa = f / aspect_ratio;
511        let bb = f;
512        let cc = z_far / (z_near - z_far);
513        let dd = (z_near * z_far) / (z_near - z_far);
514
515        Self::from_cols(
516            Vec4::new(aa, 0.0, 0.0, 0.0),
517            Vec4::new(0.0, bb, 0.0, 0.0),
518            Vec4::new(0.0, 0.0, cc, -1.0),
519            Vec4::new(0.0, 0.0, dd, 0.0),
520        )
521    }
522
523    /// Creates a right-handed orthographic projection matrix with a [0, 1] depth range (ZO).
524    #[inline]
525    pub fn orthographic_rh_zo(
526        left: f32,
527        right: f32,
528        bottom: f32,
529        top: f32,
530        z_near: f32,
531        z_far: f32,
532    ) -> Self {
533        let rml = right - left;
534        let rpl = right + left;
535        let tmb = top - bottom;
536        let tpb = top + bottom;
537        let fmn = z_far - z_near;
538        let aa = 2.0 / rml;
539        let bb = 2.0 / tmb;
540        let cc = -1.0 / fmn;
541        let dd = -rpl / rml;
542        let ee = -tpb / tmb;
543        let ff = -z_near / fmn;
544
545        Self::from_cols(
546            Vec4::new(aa, 0.0, 0.0, 0.0),
547            Vec4::new(0.0, bb, 0.0, 0.0),
548            Vec4::new(0.0, 0.0, cc, 0.0),
549            Vec4::new(dd, ee, ff, 1.0),
550        )
551    }
552
553    /// Creates a right-handed view matrix for a camera looking from `eye` towards `target`.
554    ///
555    /// # Arguments
556    ///
557    /// * `eye`: The position of the camera in world space.
558    /// * `target`: The point in world space that the camera is looking at.
559    /// * `up`: A vector indicating the "up" direction of the world (commonly `Vec3::Y`).
560    ///
561    /// # Returns
562    ///
563    /// Returns `Some(Mat4)` if a valid view matrix can be constructed, or `None` if
564    /// `eye` and `target` are too close, or if `up` is parallel to the view direction.
565    #[inline]
566    pub fn look_at_rh(eye: Vec3, target: Vec3, up: Vec3) -> Option<Self> {
567        let forward = target - eye;
568        if forward.length_squared() < crate::math::EPSILON * crate::math::EPSILON {
569            return None;
570        }
571        let f = forward.normalize();
572        let s = f.cross(up);
573        if s.length_squared() < crate::math::EPSILON * crate::math::EPSILON {
574            return None;
575        }
576        let s = s.normalize();
577        let u = s.cross(f);
578
579        Some(Self::from_cols(
580            Vec4::new(s.x, u.x, -f.x, 0.0),
581            Vec4::new(s.y, u.y, -f.y, 0.0),
582            Vec4::new(s.z, u.z, -f.z, 0.0),
583            Vec4::new(-eye.dot(s), -eye.dot(u), eye.dot(f), 1.0),
584        ))
585    }
586
587    /// Returns the transpose of the matrix, where rows and columns are swapped.
588    #[inline]
589    pub fn transpose(&self) -> Self {
590        Self::from_cols(
591            Vec4::new(
592                self.cols[0].x,
593                self.cols[1].x,
594                self.cols[2].x,
595                self.cols[3].x,
596            ),
597            Vec4::new(
598                self.cols[0].y,
599                self.cols[1].y,
600                self.cols[2].y,
601                self.cols[3].y,
602            ),
603            Vec4::new(
604                self.cols[0].z,
605                self.cols[1].z,
606                self.cols[2].z,
607                self.cols[3].z,
608            ),
609            Vec4::new(
610                self.cols[0].w,
611                self.cols[1].w,
612                self.cols[2].w,
613                self.cols[3].w,
614            ),
615        )
616    }
617
618    /// Computes the determinant of the matrix.
619    pub fn determinant(&self) -> f32 {
620        let c0 = self.cols[0];
621        let c1 = self.cols[1];
622        let c2 = self.cols[2];
623        let c3 = self.cols[3];
624
625        let m00 = c1.y * (c2.z * c3.w - c3.z * c2.w) - c2.y * (c1.z * c3.w - c3.z * c1.w)
626            + c3.y * (c1.z * c2.w - c2.z * c1.w);
627        let m01 = c0.y * (c2.z * c3.w - c3.z * c2.w) - c2.y * (c0.z * c3.w - c3.z * c0.w)
628            + c3.y * (c0.z * c2.w - c2.z * c0.w);
629        let m02 = c0.y * (c1.z * c3.w - c3.z * c1.w) - c1.y * (c0.z * c3.w - c3.z * c0.w)
630            + c3.y * (c0.z * c1.w - c1.z * c0.w);
631        let m03 = c0.y * (c1.z * c2.w - c2.z * c1.w) - c1.y * (c0.z * c2.w - c2.z * c0.w)
632            + c2.y * (c0.z * c1.w - c1.z * c0.w);
633
634        c0.x * m00 - c1.x * m01 + c2.x * m02 - c3.x * m03
635    }
636
637    /// Computes the inverse of the matrix.
638    /// Returns `None` if the matrix is not invertible.
639    pub fn inverse(&self) -> Option<Self> {
640        let c0 = self.cols[0];
641        let c1 = self.cols[1];
642        let c2 = self.cols[2];
643        let c3 = self.cols[3];
644
645        let a00 = c1.y * (c2.z * c3.w - c3.z * c2.w) - c2.y * (c1.z * c3.w - c3.z * c1.w)
646            + c3.y * (c1.z * c2.w - c2.z * c1.w);
647        let a01 = -(c1.x * (c2.z * c3.w - c3.z * c2.w) - c2.x * (c1.z * c3.w - c3.z * c1.w)
648            + c3.x * (c1.z * c2.w - c2.z * c1.w));
649        let a02 = c1.x * (c2.y * c3.w - c3.y * c2.w) - c2.x * (c1.y * c3.w - c3.y * c1.w)
650            + c3.x * (c1.y * c2.w - c2.y * c1.w);
651        let a03 = -(c1.x * (c2.y * c3.z - c3.y * c2.z) - c2.x * (c1.y * c3.z - c3.y * c1.z)
652            + c3.x * (c1.y * c2.z - c2.y * c1.z));
653
654        let a10 = -(c0.y * (c2.z * c3.w - c3.z * c2.w) - c2.y * (c0.z * c3.w - c3.z * c0.w)
655            + c3.y * (c0.z * c2.w - c2.z * c0.w));
656        let a11 = c0.x * (c2.z * c3.w - c3.z * c2.w) - c2.x * (c0.z * c3.w - c3.z * c0.w)
657            + c3.x * (c0.z * c2.w - c2.z * c0.w);
658        let a12 = -(c0.x * (c2.y * c3.w - c3.y * c2.w) - c2.x * (c0.y * c3.w - c3.y * c0.w)
659            + c3.x * (c0.y * c2.w - c2.y * c0.w));
660        let a13 = c0.x * (c2.y * c3.z - c3.y * c2.z) - c2.x * (c0.y * c3.z - c3.y * c0.z)
661            + c3.x * (c0.y * c2.z - c2.y * c0.z);
662
663        let a20 = c0.y * (c1.z * c3.w - c3.z * c1.w) - c1.y * (c0.z * c3.w - c3.z * c0.w)
664            + c3.y * (c0.z * c1.w - c1.z * c0.w);
665        let a21 = -(c0.x * (c1.z * c3.w - c3.z * c1.w) - c1.x * (c0.z * c3.w - c3.z * c0.w)
666            + c3.x * (c0.z * c1.w - c1.z * c0.w));
667        let a22 = c0.x * (c1.y * c3.w - c3.y * c1.w) - c1.x * (c0.y * c3.w - c3.y * c0.w)
668            + c3.x * (c0.y * c1.w - c1.y * c0.w);
669        let a23 = -(c0.x * (c1.y * c3.z - c3.y * c1.z) - c1.x * (c0.y * c3.z - c3.y * c0.z)
670            + c3.x * (c0.y * c1.z - c1.y * c0.z));
671
672        let a30 = -(c0.y * (c1.z * c2.w - c2.z * c1.w) - c1.y * (c0.z * c2.w - c2.z * c0.w)
673            + c2.y * (c0.z * c1.w - c1.z * c0.w));
674        let a31 = c0.x * (c1.z * c2.w - c2.z * c1.w) - c1.x * (c0.z * c2.w - c2.z * c0.w)
675            + c2.x * (c0.z * c1.w - c1.z * c0.w);
676        let a32 = -(c0.x * (c1.y * c2.w - c2.y * c1.w) - c1.x * (c0.y * c2.w - c2.y * c0.w)
677            + c2.x * (c0.y * c1.w - c1.y * c0.w));
678        let a33 = c0.x * (c1.y * c2.z - c2.y * c1.z) - c1.x * (c0.y * c2.z - c2.y * c0.z)
679            + c2.x * (c0.y * c1.z - c1.y * c0.z);
680
681        let det = c0.x * a00 + c1.x * a10 + c2.x * a20 + c3.x * a30;
682        if det.abs() < crate::math::EPSILON {
683            return None;
684        }
685        let inv_det = 1.0 / det;
686
687        Some(Self::from_cols(
688            Vec4::new(a00 * inv_det, a10 * inv_det, a20 * inv_det, a30 * inv_det),
689            Vec4::new(a01 * inv_det, a11 * inv_det, a21 * inv_det, a31 * inv_det),
690            Vec4::new(a02 * inv_det, a12 * inv_det, a22 * inv_det, a32 * inv_det),
691            Vec4::new(a03 * inv_det, a13 * inv_det, a23 * inv_det, a33 * inv_det),
692        ))
693    }
694
695    /// Computes the inverse of an affine transformation matrix more efficiently
696    /// and with better numerical stability than the general `inverse` method.
697    ///
698    /// An affine matrix is one composed of only translation, rotation, and scale.
699    ///
700    /// # Returns
701    ///
702    /// `None` if the matrix is not affine or is not invertible.
703    #[inline]
704    pub fn affine_inverse(&self) -> Option<Self> {
705        let c0 = self.cols[0].truncate();
706        let c1 = self.cols[1].truncate();
707        let c2 = self.cols[2].truncate();
708        let translation = self.cols[3].truncate();
709        let det3x3 = c0.x * (c1.y * c2.z - c2.y * c1.z) - c1.x * (c0.y * c2.z - c2.y * c0.z)
710            + c2.x * (c0.y * c1.z - c1.y * c0.z);
711
712        if det3x3.abs() < crate::math::EPSILON {
713            return None;
714        }
715
716        let inv_det3x3 = 1.0 / det3x3;
717        let inv00 = (c1.y * c2.z - c2.y * c1.z) * inv_det3x3;
718        let inv10 = -(c2.y * c0.z - c0.y * c2.z) * inv_det3x3;
719        let inv20 = (c0.y * c1.z - c1.y * c0.z) * inv_det3x3;
720        let inv01 = -(c2.x * c1.z - c1.x * c2.z) * inv_det3x3;
721        let inv11 = (c0.x * c2.z - c2.x * c0.z) * inv_det3x3;
722        let inv21 = -(c1.x * c0.z - c0.x * c1.z) * inv_det3x3;
723        let inv02 = (c1.x * c2.y - c2.x * c1.y) * inv_det3x3;
724        let inv12 = -(c2.x * c0.y - c0.x * c2.y) * inv_det3x3;
725        let inv22 = (c0.x * c1.y - c1.x * c0.y) * inv_det3x3;
726        let inv_tx = -(inv00 * translation.x + inv01 * translation.y + inv02 * translation.z);
727        let inv_ty = -(inv10 * translation.x + inv11 * translation.y + inv12 * translation.z);
728        let inv_tz = -(inv20 * translation.x + inv21 * translation.y + inv22 * translation.z);
729
730        Some(Self::from_cols(
731            Vec4::new(inv00, inv10, inv20, 0.0),
732            Vec4::new(inv01, inv11, inv21, 0.0),
733            Vec4::new(inv02, inv12, inv22, 0.0),
734            Vec4::new(inv_tx, inv_ty, inv_tz, 1.0),
735        ))
736    }
737}
738
739// --- Operators Overloading ---
740
741impl Default for Mat4 {
742    /// Returns the 4x4 identity matrix.
743    #[inline]
744    fn default() -> Self {
745        Self::IDENTITY
746    }
747}
748
749impl Mul<Mat4> for Mat4 {
750    type Output = Self;
751    /// Multiplies this matrix by another `Mat4`. Note that matrix multiplication is not commutative.
752    #[inline]
753    fn mul(self, rhs: Mat4) -> Self::Output {
754        let mut result_cols = [Vec4 {
755            x: 0.0,
756            y: 0.0,
757            z: 0.0,
758            w: 0.0,
759        }; 4];
760        for (c_idx, target_col_ref_mut) in result_cols.iter_mut().enumerate().take(4) {
761            let col_from_rhs = rhs.cols[c_idx];
762            *target_col_ref_mut = Vec4 {
763                x: self.get_row(0).dot(col_from_rhs),
764                y: self.get_row(1).dot(col_from_rhs),
765                z: self.get_row(2).dot(col_from_rhs),
766                w: self.get_row(3).dot(col_from_rhs),
767            };
768        }
769        Mat4 { cols: result_cols }
770    }
771}
772
773impl Mul<Vec4> for Mat4 {
774    type Output = Vec4;
775    /// Transforms a `Vec4` by this matrix.
776    #[inline]
777    fn mul(self, rhs: Vec4) -> Self::Output {
778        self.cols[0] * rhs.x + self.cols[1] * rhs.y + self.cols[2] * rhs.z + self.cols[3] * rhs.w
779    }
780}
781
782// --- Tests ---
783
784#[cfg(test)]
785mod tests {
786    use super::*;
787    use crate::math::{approx_eq, matrix::Mat4, quaternion::Quaternion, vector::Vec3, PI};
788
789    fn vec3_approx_eq(a: Vec3, b: Vec3) -> bool {
790        approx_eq(a.x, b.x) && approx_eq(a.y, b.y) && approx_eq(a.z, b.z)
791    }
792
793    fn mat3_approx_eq(a: Mat3, b: Mat3) -> bool {
794        vec3_approx_eq(a.cols[0], b.cols[0])
795            && vec3_approx_eq(a.cols[1], b.cols[1])
796            && vec3_approx_eq(a.cols[2], b.cols[2])
797    }
798
799    fn vec4_approx_eq(a: Vec4, b: Vec4) -> bool {
800        approx_eq(a.x, b.x) && approx_eq(a.y, b.y) && approx_eq(a.z, b.z) && approx_eq(a.w, b.w)
801    }
802
803    fn mat4_approx_eq(a: Mat4, b: Mat4) -> bool {
804        vec4_approx_eq(a.cols[0], b.cols[0])
805            && vec4_approx_eq(a.cols[1], b.cols[1])
806            && vec4_approx_eq(a.cols[2], b.cols[2])
807            && vec4_approx_eq(a.cols[3], b.cols[3])
808    }
809
810    // --- Tests for Mat3 ---
811
812    #[test]
813    fn test_mat3_identity_default() {
814        assert_eq!(Mat3::default(), Mat3::IDENTITY);
815
816        let m = Mat3::from_scale(Vec3::new(1.0, 2.0, 3.0));
817        assert!(mat3_approx_eq(m * Mat3::IDENTITY, m));
818        assert!(mat3_approx_eq(Mat3::IDENTITY * m, m));
819    }
820
821    #[test]
822    fn test_mat3_from_scale() {
823        let s = Vec3::new(2.0, -3.0, 0.5);
824        let m = Mat3::from_scale(s);
825        let v = Vec3::new(1.0, 1.0, 1.0);
826        assert!(vec3_approx_eq(m * v, s)); // Scaling (1,1,1) should yield the scale vector
827    }
828
829    #[test]
830    fn test_mat3_rotations() {
831        let angle = PI / 6.0; // 30 degrees
832        let mx = Mat3::from_rotation_x(angle);
833        let my = Mat3::from_rotation_y(angle);
834        let mz = Mat3::from_rotation_z(angle);
835
836        let p = Vec3::Y; // Point on Y axis
837        let expected_px = Vec3::new(0.0, angle.cos(), angle.sin());
838        assert!(vec3_approx_eq(mx * p, expected_px));
839
840        let p = Vec3::X; // Point on X axis
841        let expected_py = Vec3::new(angle.cos(), 0.0, -angle.sin()); // RH
842        assert!(vec3_approx_eq(my * p, expected_py));
843
844        let p = Vec3::X; // Point on X axis
845        let expected_pz = Vec3::new(angle.cos(), angle.sin(), 0.0);
846        assert!(vec3_approx_eq(mz * p, expected_pz));
847    }
848
849    #[test]
850    fn test_mat3_from_axis_angle() {
851        let axis = Vec3::new(1.0, 1.0, 1.0).normalize();
852        let angle = 1.2 * PI;
853        let m = Mat3::from_axis_angle(axis, angle);
854
855        // Rotation around (1,1,1) axis should permute basis vectors
856        let v = Vec3::X;
857        let v_rotated = m * v;
858
859        // Check if length is preserved
860        assert!(approx_eq(v_rotated.length(), v.length()));
861
862        // Specific check is hard without known values, but ensure it's not identity or zero
863        assert!(v_rotated.distance_squared(v) > EPSILON);
864    }
865
866    #[test]
867    fn test_mat3_from_quat() {
868        let axis = Vec3::new(1.0, -2.0, 3.0).normalize();
869        let angle = PI / 7.0;
870        let q = Quaternion::from_axis_angle(axis, angle);
871        let m_from_q = Mat3::from_quat(q);
872
873        let v = Vec3::new(0.5, 1.0, -0.2);
874        let v_rotated_q = q * v;
875        let v_rotated_m = m_from_q * v;
876
877        assert!(vec3_approx_eq(v_rotated_q, v_rotated_m));
878    }
879
880    #[test]
881    fn test_mat3_determinant() {
882        assert!(approx_eq(Mat3::IDENTITY.determinant(), 1.0));
883        assert!(approx_eq(Mat3::ZERO.determinant(), 0.0));
884
885        let m_scale = Mat3::from_scale(Vec3::new(2.0, 3.0, 4.0));
886        assert!(approx_eq(m_scale.determinant(), 24.0));
887
888        let m_rot = Mat3::from_rotation_y(PI / 5.0);
889        assert!(approx_eq(m_rot.determinant(), 1.0)); // Rotations preserve volume
890    }
891
892    #[test]
893    fn test_mat3_transpose() {
894        let m = Mat3::from_cols(
895            Vec3::new(1.0, 2.0, 3.0),
896            Vec3::new(4.0, 5.0, 6.0),
897            Vec3::new(7.0, 8.0, 9.0),
898        );
899        let mt = m.transpose();
900        let expected_mt = Mat3::from_cols(
901            Vec3::new(1.0, 4.0, 7.0),
902            Vec3::new(2.0, 5.0, 8.0),
903            Vec3::new(3.0, 6.0, 9.0),
904        );
905
906        assert!(mat3_approx_eq(mt, expected_mt));
907        assert!(mat3_approx_eq(m.transpose().transpose(), m)); // Double transpose
908    }
909
910    #[test]
911    fn test_mat3_inverse() {
912        let m = Mat3::from_rotation_z(PI / 3.0) * Mat3::from_scale(Vec3::new(1.0, 2.0, 0.5));
913        let inv_m = m.inverse().expect("Matrix should be invertible");
914        let identity = m * inv_m;
915        assert!(
916            mat3_approx_eq(identity, Mat3::IDENTITY),
917            "M * inv(M) should be Identity"
918        );
919
920        let singular = Mat3::from_scale(Vec3::new(1.0, 0.0, 1.0));
921        assert!(
922            singular.inverse().is_none(),
923            "Singular matrix inverse should be None"
924        );
925    }
926
927    #[test]
928    fn test_mat3_mul_vec3() {
929        let m = Mat3::from_rotation_z(PI / 2.0); // Rotate 90 deg around Z
930        let v = Vec3::X; // (1, 0, 0)
931        let expected_v = Vec3::Y; // (0, 1, 0)
932        assert!(vec3_approx_eq(m * v, expected_v));
933    }
934
935    #[test]
936    fn test_mat3_mul_mat3() {
937        let rot90z = Mat3::from_rotation_z(PI / 2.0);
938        let rot180z = rot90z * rot90z;
939        let expected_rot180z = Mat3::from_rotation_z(PI);
940        assert!(mat3_approx_eq(rot180z, expected_rot180z));
941    }
942
943    #[test]
944    fn test_mat3_conversions() {
945        let m4 = Mat4::from_translation(Vec3::new(10., 20., 30.)) * Mat4::from_rotation_x(PI / 4.0);
946        let m3 = Mat3::from_mat4(&m4);
947        let m4_again = m3.to_mat4();
948
949        // Check if rotation part was extracted correctly
950        let v = Vec3::Y;
951        let v_rot_m3 = m3 * v;
952        let v_rot_m4 = Mat4::from_rotation_x(PI / 4.0) * Vec4::from_vec3(v, 0.0); // Rotate as vector
953        assert!(vec3_approx_eq(v_rot_m3, v_rot_m4.truncate()));
954
955        // Check if embedding back into Mat4 worked (translation should be zero)
956        let origin = Vec4::new(0.0, 0.0, 0.0, 1.0);
957        let transformed_origin = m4_again * origin;
958        assert!(approx_eq(transformed_origin.x, 0.0));
959        assert!(approx_eq(transformed_origin.y, 0.0));
960        assert!(approx_eq(transformed_origin.z, 0.0));
961        assert!(approx_eq(transformed_origin.w, 1.0));
962    }
963
964    #[test]
965    fn test_mat3_index() {
966        let mut m = Mat3::from_cols(Vec3::X, Vec3::Y, Vec3::Z);
967        assert_eq!(m[0], Vec3::X);
968        assert_eq!(m[1], Vec3::Y);
969        assert_eq!(m[2], Vec3::Z);
970        m[0] = Vec3::ONE;
971        assert_eq!(m.cols[0], Vec3::ONE);
972    }
973
974    #[test]
975    #[should_panic]
976    fn test_mat3_index_out_of_bounds() {
977        let m = Mat3::IDENTITY;
978        let _ = m[3]; // Should panic
979    }
980
981    // --- End of Mat3 Tests ---
982
983    // --- Tests for Mat4 ---
984
985    #[test]
986    fn test_identity() {
987        assert_eq!(Mat4::default(), Mat4::IDENTITY);
988        let m = Mat4::from_translation(Vec3::new(1.0, 2.0, 3.0));
989        assert!(mat4_approx_eq(m * Mat4::IDENTITY, m));
990        assert!(mat4_approx_eq(Mat4::IDENTITY * m, m));
991    }
992
993    #[test]
994    fn test_from_quat() {
995        let axis = Vec3::new(1.0, 2.0, 3.0).normalize();
996        let angle = PI / 5.0;
997        let q = Quaternion::from_axis_angle(axis, angle);
998        let m_from_q = Mat4::from_quat(q);
999
1000        let v = Vec3::new(5.0, -1.0, 2.0);
1001
1002        let v_rotated_q = q * v; // Rotate using quaternion directly
1003                                 // Rotate using matrix: convert v to Vec4(point), multiply, convert back
1004        let v4 = Vec4::from_vec3(v, 1.0);
1005        let v_rotated_m4 = m_from_q * v4;
1006        let v_rotated_m = v_rotated_m4.truncate();
1007
1008        // Compare results
1009        assert!(approx_eq(v_rotated_q.x, v_rotated_m.x));
1010        assert!(approx_eq(v_rotated_q.y, v_rotated_m.y));
1011        assert!(approx_eq(v_rotated_q.z, v_rotated_m.z));
1012    }
1013
1014    #[test]
1015    fn test_translation() {
1016        let t = Vec3::new(1.0, 2.0, 3.0);
1017        let m = Mat4::from_translation(t);
1018        let p = Vec4::new(1.0, 1.0, 1.0, 1.0);
1019        let expected_p = Vec4::new(2.0, 3.0, 4.0, 1.0);
1020
1021        assert!(vec4_approx_eq(m * p, expected_p));
1022    }
1023
1024    #[test]
1025    fn test_scale() {
1026        let s = Vec3::new(2.0, 3.0, 4.0);
1027        let m = Mat4::from_scale(s);
1028        let p = Vec4::new(1.0, 1.0, 1.0, 1.0);
1029        let expected_p = Vec4::new(2.0, 3.0, 4.0, 1.0);
1030        assert!(vec4_approx_eq(m * p, expected_p));
1031    }
1032
1033    #[test]
1034    fn test_rotation_x() {
1035        let angle = PI / 2.0; // 90 degrees
1036        let m = Mat4::from_rotation_x(angle);
1037        let p = Vec4::new(0.0, 1.0, 0.0, 1.0); // Point on Y axis
1038        let expected_p = Vec4::new(0.0, 0.0, 1.0, 1.0); // Should rotate to Z axis
1039        assert!(vec4_approx_eq(m * p, expected_p));
1040    }
1041
1042    #[test]
1043    fn test_rotation_y() {
1044        let angle = PI / 2.0; // 90 degrees
1045        let m = Mat4::from_rotation_y(angle);
1046        let p = Vec4::new(1.0, 0.0, 0.0, 1.0); // Point on X axis
1047        let expected_p = Vec4::new(0.0, 0.0, -1.0, 1.0); // Should rotate to -Z axis
1048
1049        assert!(vec4_approx_eq(m * p, expected_p));
1050    }
1051
1052    #[test]
1053    fn test_rotation_z() {
1054        let angle = PI / 2.0; // 90 degrees
1055        let m = Mat4::from_rotation_z(angle);
1056        let p = Vec4::new(1.0, 0.0, 0.0, 1.0); // Point on X axis
1057        let expected_p = Vec4::new(0.0, 1.0, 0.0, 1.0); // Should rotate to Y axis
1058        assert!(vec4_approx_eq(m * p, expected_p));
1059    }
1060
1061    #[test]
1062    fn test_transpose() {
1063        let m = Mat4::from_cols(
1064            Vec4::new(1., 2., 3., 4.),
1065            Vec4::new(5., 6., 7., 8.),
1066            Vec4::new(9., 10., 11., 12.),
1067            Vec4::new(13., 14., 15., 16.),
1068        );
1069        let mt = m.transpose();
1070        let expected_mt = Mat4::from_cols(
1071            Vec4::new(1., 5., 9., 13.),
1072            Vec4::new(2., 6., 10., 14.),
1073            Vec4::new(3., 7., 11., 15.),
1074            Vec4::new(4., 8., 12., 16.),
1075        );
1076        assert_eq!(mt.cols[0], expected_mt.cols[0]); // Compare columns after transpose
1077        assert_eq!(mt.cols[1], expected_mt.cols[1]);
1078        assert_eq!(mt.cols[2], expected_mt.cols[2]);
1079        assert_eq!(mt.cols[3], expected_mt.cols[3]);
1080
1081        // Test double transpose
1082        assert!(mat4_approx_eq(m.transpose().transpose(), m));
1083    }
1084
1085    #[test]
1086    fn test_mul_mat4() {
1087        let t = Mat4::from_translation(Vec3::new(1.0, 0.0, 0.0));
1088        let r = Mat4::from_rotation_z(PI / 2.0);
1089
1090        // Order matters: Translate then Rotate
1091        let tr = r * t;
1092        let p = Vec4::new(1.0, 0.0, 0.0, 1.0); // Point at (1,0,0)
1093                                               // 1. Translate: p becomes (2, 0, 0, 1)
1094                                               // 2. Rotate Z 90: (2, 0, 0) becomes (0, 2, 0)
1095        let expected_tr = Vec4::new(0.0, 2.0, 0.0, 1.0);
1096        assert!(vec4_approx_eq(tr * p, expected_tr));
1097
1098        // Order matters: Rotate then Translate
1099        let rt = t * r;
1100        // 1. Rotate Z 90: p becomes (0, 1, 0, 1)
1101        // 2. Translate: (0, 1, 0) becomes (1, 1, 0)
1102        let expected_rt = Vec4::new(1.0, 1.0, 0.0, 1.0);
1103        assert!(vec4_approx_eq(rt * p, expected_rt));
1104    }
1105
1106    #[test]
1107    fn test_inverse() {
1108        let m = Mat4::from_translation(Vec3::new(1., 2., 3.))
1109            * Mat4::from_rotation_y(PI / 4.0)
1110            * Mat4::from_scale(Vec3::new(1., 2., 1.));
1111
1112        let inv_m = m.inverse().expect("Matrix should be invertible");
1113        let identity = m * inv_m;
1114
1115        // Check if M * M^-1 is close to identity
1116        assert!(
1117            mat4_approx_eq(identity, Mat4::IDENTITY),
1118            "M * inv(M) should be Identity"
1119        );
1120
1121        // Check singular matrix (e.g., scale with zero)
1122        let singular = Mat4::from_scale(Vec3::new(1.0, 0.0, 1.0));
1123        assert!(
1124            singular.inverse().is_none(),
1125            "Singular matrix inverse should be None"
1126        );
1127    }
1128
1129    #[test]
1130    fn test_affine_inverse() {
1131        let t = Mat4::from_translation(Vec3::new(1., 2., 3.));
1132        let r = Mat4::from_rotation_y(PI / 3.0);
1133        let s = Mat4::from_scale(Vec3::new(1., 2., 0.5));
1134        let m = t * r * s; // Combined affine transform
1135
1136        let inv_m = m.inverse().expect("Matrix should be invertible");
1137        let affine_inv_m = m
1138            .affine_inverse()
1139            .expect("Matrix should be affine invertible");
1140
1141        // Check if affine inverse matches general inverse for this case
1142        assert!(
1143            mat4_approx_eq(inv_m, affine_inv_m),
1144            "Affine inverse should match general inverse"
1145        );
1146
1147        // Check M * inv(M) == Identity using affine inverse
1148        let identity = m * affine_inv_m;
1149        assert!(
1150            mat4_approx_eq(identity, Mat4::IDENTITY),
1151            "M * affine_inv(M) should be Identity"
1152        );
1153
1154        // Test singular affine matrix
1155        let singular_s = Mat4::from_scale(Vec3::new(1.0, 0.0, 1.0));
1156        let singular_m = t * singular_s;
1157        assert!(
1158            singular_m.affine_inverse().is_none(),
1159            "Singular affine matrix inverse should be None"
1160        );
1161    }
1162
1163    #[test]
1164    fn test_perspective_rh_zo() {
1165        let fov = PI / 4.0; // 45 degrees
1166        let aspect = 16.0 / 9.0;
1167        let near = 0.1;
1168        let far = 100.0;
1169
1170        let m = Mat4::perspective_rh_zo(fov, aspect, near, far);
1171        assert!(approx_eq(m.cols[0].x, 1.0 / (aspect * (fov / 2.0).tan())));
1172        assert!(approx_eq(m.cols[1].y, 1.0 / ((fov / 2.0).tan())));
1173        assert!(approx_eq(m.cols[2].z, -far / (far - near)));
1174        assert!(approx_eq(m.cols[3].z, -(far * near) / (far - near)));
1175    }
1176
1177    #[test]
1178    fn test_orthographic_rh_zo() {
1179        let left = -1.0;
1180        let right = 1.0;
1181        let bottom = -1.0;
1182        let top = 1.0;
1183        let near = 0.1;
1184        let far = 100.0;
1185        let m = Mat4::orthographic_rh_zo(left, right, bottom, top, near, far);
1186
1187        // Check scale factors
1188        assert!(approx_eq(m.cols[0].x, 2.0 / (right - left)));
1189        assert!(approx_eq(m.cols[1].y, 2.0 / (top - bottom)));
1190        assert!(approx_eq(m.cols[2].z, -1.0 / (far - near)));
1191
1192        // Check translation factors
1193        assert!(approx_eq(m.cols[3].x, -(right + left) / (right - left)));
1194        assert!(approx_eq(m.cols[3].y, -(top + bottom) / (top - bottom)));
1195        assert!(approx_eq(m.cols[3].z, -near / (far - near))); // -near and not -(far + near)
1196    }
1197
1198    #[test]
1199    fn test_look_at_rh() {
1200        let eye = Vec3::new(0.0, 0.0, 5.0);
1201        let target = Vec3::new(0.0, 0.0, 0.0);
1202        let up = Vec3::new(0.0, 1.0, 0.0);
1203
1204        let m = Mat4::look_at_rh(eye, target, up).expect("look_at_rh should return Some(Mat4)");
1205
1206        // Forward direction (third column, third row): should be +1.0 for a right-handed system
1207        assert!(approx_eq(m.cols[2].z, 1.0));
1208
1209        // Translation part (fourth column, third row): should be -eye ยท forward = -5.0
1210        assert!(approx_eq(m.cols[3].z, -5.0));
1211    }
1212
1213    #[test]
1214    fn test_look_at_rh_invalid() {
1215        let eye = Vec3::new(0.0, 0.0, 5.0);
1216        let target = Vec3::new(0.0, 0.0, 5.0); // Same as eye
1217        let up = Vec3::new(0.0, 1.0, 0.0);
1218
1219        // This should panic or return None (depending on implementation)
1220        assert!(Mat4::look_at_rh(eye, target, up).is_none());
1221    }
1222
1223    // --- End of Tests For Mat4 ---
1224}