docs: update ray sphere intersection brief

clarify some rotation rules

CREDITS

@@ -82,3 +82,17 @@ http://github.com/microsoft/DirectXMath
 17. Pick Matrix
 
 glu project -> project.c
+
+18. Ray sphere intersection
+
+RAY TRACING GEMS
+HIGH-QUALITY AND REAL-TIME RENDERING WITH DXR AND OTHER APIS
+
+CHAPTER 7
+Precision Improvements for Ray/Sphere Intersection
+Eric Haines (1), Johannes Günther (2), and Tomas Akenine-Möller (1)
+  (1) NVIDIA
+  (2) Intel
+
+Wyman, C., and Haines, E. Getting Started with RTX Ray Tracing.
+https://github.com/NVIDIAGameWorks/GettingStartedWithRTXRayTracing

Makefile.am

@@ -216,7 +216,8 @@ cglm_struct_HEADERS = include/cglm/struct/mat4.h \
                       include/cglm/struct/project.h \
                       include/cglm/struct/sphere.h \
                       include/cglm/struct/color.h \
-                      include/cglm/struct/curve.h
+                      include/cglm/struct/curve.h \
+                      include/cglm/struct/ray.h
 
 cglm_struct_clipspacedir=$(includedir)/cglm/struct/clipspace
 cglm_struct_clipspace_HEADERS = include/cglm/struct/clipspace/persp_lh_no.h \

docs/source/affine-pre.rst

@@ -110,7 +110,24 @@ Functions documentation
 
 .. c:function:: void  glm_rotate(mat4 m, float angle, vec3 axis)
 
-    rotate existing transform matrix around Z axis by angle and axis
+    rotate existing transform matrix around given axis by angle at ORIGIN (0,0,0)
+
+    **❗️IMPORTANT ❗️**
+    
+    If you need to rotate object around itself e.g. center of object or at
+    some point [of object] then `glm_rotate_at()` would be better choice to do so.
+    
+    Even if object's model transform is identiy, rotation may not be around
+    center of object if object does not lay out at ORIGIN perfectly.
+    
+    Using `glm_rotate_at()` with center of bounding shape ( AABB, Sphere ... )
+    would be an easy option to rotate around object if object is not at origin.
+    
+    One another option to rotate around itself at any point is `glm_spin()`
+    which is perfect if only rotating around model position is desired e.g. not
+    specific point on model for instance center of geometry or center of mass,
+    again if geometry is not perfectly centered at origin at identity transform,
+    rotation may not be around geometry.
 
     Parameters:
       | *[in, out]* **m**     affine transform

docs/source/mat2.rst

@@ -180,7 +180,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat2_make(float * __restrict src, mat2 dest)
+.. c:function:: void glm_mat2_make(const float * __restrict src, mat2 dest)
 
     Create mat2 matrix from pointer
 

docs/source/mat2x3.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat2x3_make(float * __restrict src, mat2x3 dest)
+.. c:function:: void glm_mat2x3_make(const float * __restrict src, mat2x3 dest)
 
     Create mat2x3 matrix from pointer
 

docs/source/mat2x4.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat2x4_make(float * __restrict src, mat2x4 dest)
+.. c:function:: void glm_mat2x4_make(const float * __restrict src, mat2x4 dest)
 
     Create mat2x4 matrix from pointer
 

docs/source/mat3.rst

@@ -190,7 +190,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat3_make(float * __restrict src, mat3 dest)
+.. c:function:: void glm_mat3_make(const float * __restrict src, mat3 dest)
 
     Create mat3 matrix from pointer
 

docs/source/mat3x2.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat3x2_make(float * __restrict src, mat3x2 dest)
+.. c:function:: void glm_mat3x2_make(const float * __restrict src, mat3x2 dest)
 
     Create mat3x2 matrix from pointer
 

docs/source/mat3x4.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat3x4_make(float * __restrict src, mat3x4 dest)
+.. c:function:: void glm_mat3x4_make(const float * __restrict src, mat3x4 dest)
 
     Create mat3x4 matrix from pointer
 

docs/source/mat4.rst

@@ -304,7 +304,7 @@ Functions documentation
     Returns:
         scalar value e.g. Matrix1x1
 
-.. c:function:: void glm_mat4_make(float * __restrict src, mat4 dest)
+.. c:function:: void glm_mat4_make(const float * __restrict src, mat4 dest)
 
     Create mat4 matrix from pointer
 

docs/source/mat4x2.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat4x2_make(float * __restrict src, mat4x2 dest)
+.. c:function:: void glm_mat4x2_make(const float * __restrict src, mat4x2 dest)
 
     Create mat4x2 matrix from pointer
 

docs/source/mat4x3.rst

@@ -41,7 +41,7 @@ Functions documentation
     Parameters:
       | *[in,out]* **mat**  matrix
 
-.. c:function:: void glm_mat4x3_make(float * __restrict src, mat4x3 dest)
+.. c:function:: void glm_mat4x3_make(const float * __restrict src, mat4x3 dest)
 
     Create mat4x3 matrix from pointer
 

docs/source/quat.rst

@@ -422,7 +422,7 @@ Functions documentation
       | *[in]*       **q**      quaternion
       | *[in]*       **pivot**  pivot
 
-.. c:function:: void glm_quat_make(float * __restrict src, versor dest)
+.. c:function:: void glm_quat_make(const float * __restrict src, versor dest)
 
     Create quaternion from pointer
 

docs/source/ray.rst

@@ -13,6 +13,8 @@ Table of contents (click to go):
 Functions:
 
 1. :c:func:`glm_ray_triangle`
+#. :c:func:`glm_ray_sphere`
+#. :c:func:`glm_ray_at`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -29,3 +31,37 @@ Functions documentation
       | *[in]*       **v2**            third vertex of triangle
       | *[in, out]*  **d**             float pointer to save distance to intersection
       | *[out]*      **intersection**  whether there is intersection
+
+.. c:function:: bool glm_ray_sphere(vec3 origin, vec3 dir, vec4 s, float * __restrict t1, float * __restrict t2)
+
+    ray sphere intersection
+
+    - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+    - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+    - t1 < 0, t2 < 0: no intersection ahead of the ray
+    - the caller can check if the intersection points (t1 and t2) fall within a
+      specific range (for example, tmin < t1, t2 < tmax) to determine if the
+      intersections are within a desired segment of the ray
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     normalized ray direction
+      | *[in]*   **s**       sphere  [center.x, center.y, center.z, radii]
+      | *[out]*  **t1**      near point1 (closer to origin)
+      | *[out]*  **t2**      far point2 (farther from origin)
+
+    Return:
+      | whether there is intersection
+
+.. c:function:: bool  glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point)
+
+    point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     ray direction
+      | *[out]*  **t**       parameter
+      | *[out]*  **point**   point at t
+
+    Return:
+      | point at t

docs/source/vec2.rst

@@ -53,6 +53,8 @@ Functions:
 #. :c:func:`glm_vec2_clamp`
 #. :c:func:`glm_vec2_lerp`
 #. :c:func:`glm_vec2_make`
+#. :c:func:`glm_vec2_reflect`
+#. :c:func:`glm_vec2_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -385,7 +387,7 @@ Functions documentation
       | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec2_make(float * __restrict src, vec2 dest)
+.. c:function:: void glm_vec2_make(const float * __restrict src, vec2 dest)
 
     Create two dimensional vector from pointer
 
@@ -394,3 +396,25 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec2_reflect(vec2 I, vec2 N, vec2 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

docs/source/vec3.rst

@@ -80,6 +80,9 @@ Functions:
 #. :c:func:`glm_vec3_clamp`
 #. :c:func:`glm_vec3_lerp`
 #. :c:func:`glm_vec3_make`
+#. :c:func:`glm_vec3_faceforward`
+#. :c:func:`glm_vec3_reflect`
+#. :c:func:`glm_vec3_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -503,7 +506,7 @@ Functions documentation
       | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec3_make(float * __restrict src, vec3 dest)
+.. c:function:: void glm_vec3_make(const float * __restrict src, vec3 dest)
 
     Create three dimensional vector from pointer
 
@@ -512,3 +515,35 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest)
+
+    A vector pointing in the same direction as another
+
+    Parameters:
+      | *[in]*  **N**     vector to orient
+      | *[in]*  **I**     incident vector
+      | *[in]*  **Nref**  reference vector
+      | *[out]* **dest**  destination: oriented vector, pointing away from the surface.
+
+.. c:function:: void glm_vec3_reflect(vec3 I, vec3 N, vec3 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

docs/source/vec4.rst

@@ -60,6 +60,8 @@ Functions:
 #. :c:func:`glm_vec4_lerp`
 #. :c:func:`glm_vec4_cubic`
 #. :c:func:`glm_vec4_make`
+#. :c:func:`glm_vec4_reflect`
+#. :c:func:`glm_vec4_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -415,7 +417,7 @@ Functions documentation
       | *[in]*  **s**      parameter
       | *[out]* **dest**   destination
 
-.. c:function:: void glm_vec4_make(float * __restrict src, vec4 dest)
+.. c:function:: void glm_vec4_make(const float * __restrict src, vec4 dest)
 
     Create four dimensional vector from pointer
 
@@ -424,3 +426,29 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec4_reflect(vec4 I, vec4 N, vec4 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    this implementation does not explicitly preserve the 'w' component of the
+    incident vector 'I' in the output 'dest', users requiring the preservation of
+    the 'w' component should manually adjust 'dest' after calling this function.
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

include/cglm/affine-pre.h

@@ -207,7 +207,25 @@ glm_rotate_z(mat4 m, float angle, mat4 dest) {
 }
 
 /*!
- * @brief rotate existing transform matrix around given axis by angle
+ * @brief rotate existing transform matrix 
+ *        around given axis by angle at ORIGIN (0,0,0)
+ *
+ *   **❗️IMPORTANT ❗️**
+ *
+ *   If you need to rotate object around itself e.g. center of object or at
+ *   some point [of object] then `glm_rotate_at()` would be better choice to do so.
+ *
+ *   Even if object's model transform is identiy, rotation may not be around
+ *   center of object if object does not lay out at ORIGIN perfectly.
+ *
+ *   Using `glm_rotate_at()` with center of bounding shape ( AABB, Sphere ... )
+ *   would be an easy option to rotate around object if object is not at origin.
+ *
+ *   One another option to rotate around itself at any point is `glm_spin()`
+ *   which is perfect if only rotating around model position is desired e.g. not
+ *   specific point on model for instance center of geometry or center of mass,
+ *   again if geometry is not perfectly centered at origin at identity transform,
+ *   rotation may not be around geometry.
  *
  * @param[in, out]  m      affine transform
  * @param[in]       angle  angle (radians)
@@ -268,7 +286,8 @@ glm_rotate_atm(mat4 m, vec3 pivot, float angle, vec3 axis) {
 }
 
 /*!
- * @brief rotate existing transform matrix around given axis by angle around self (doesn't affected by position)
+ * @brief rotate existing transform matrix 
+ *        around given axis by angle around self (doesn't affected by position)
  *
  * @param[in, out]  m      affine transform
  * @param[in]       angle  angle (radians)

include/cglm/call/mat2.h

@@ -75,7 +75,7 @@ glmc_mat2_rmc(vec2 r, mat2 m, vec2 c);
 
 CGLM_EXPORT
 void
-glmc_mat2_make(float * __restrict src, mat2 dest);
+glmc_mat2_make(const float * __restrict src, mat2 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat2x3.h

@@ -23,7 +23,7 @@ glmc_mat2x3_zero(mat2x3 mat);
 
 CGLM_EXPORT
 void
-glmc_mat2x3_make(float * __restrict src, mat2x3 dest);
+glmc_mat2x3_make(const float * __restrict src, mat2x3 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat2x4.h

@@ -23,7 +23,7 @@ glmc_mat2x4_zero(mat2x4 mat);
 
 CGLM_EXPORT
 void
-glmc_mat2x4_make(float * __restrict src, mat2x4 dest);
+glmc_mat2x4_make(const float * __restrict src, mat2x4 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat3.h

@@ -82,7 +82,7 @@ glmc_mat3_rmc(vec3 r, mat3 m, vec3 c);
 
 CGLM_EXPORT
 void
-glmc_mat3_make(float * __restrict src, mat3 dest);
+glmc_mat3_make(const float * __restrict src, mat3 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat3x2.h

@@ -23,7 +23,7 @@ glmc_mat3x2_zero(mat3x2 mat);
 
 CGLM_EXPORT
 void
-glmc_mat3x2_make(float * __restrict src, mat3x2 dest);
+glmc_mat3x2_make(const float * __restrict src, mat3x2 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat3x4.h

@@ -23,7 +23,7 @@ glmc_mat3x4_zero(mat3x4 mat);
 
 CGLM_EXPORT
 void
-glmc_mat3x4_make(float * __restrict src, mat3x4 dest);
+glmc_mat3x4_make(const float * __restrict src, mat3x4 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat4.h

@@ -123,7 +123,7 @@ glmc_mat4_rmc(vec4 r, mat4 m, vec4 c);
 
 CGLM_EXPORT
 void
-glmc_mat4_make(float * __restrict src, mat4 dest);
+glmc_mat4_make(const float * __restrict src, mat4 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/mat4x2.h

@@ -23,7 +23,7 @@ glmc_mat4x2_zero(mat4x2 mat);
 
 CGLM_EXPORT
 void
-glmc_mat4x2_make(float * __restrict src, mat4x2 dest);
+glmc_mat4x2_make(const float * __restrict src, mat4x2 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/mat4x3.h

@@ -23,7 +23,7 @@ glmc_mat4x3_zero(mat4x3 mat);
 
 CGLM_EXPORT
 void
-glmc_mat4x3_make(float * __restrict src, mat4x3 dest);
+glmc_mat4x3_make(const float * __restrict src, mat4x3 dest);
 
 CGLM_EXPORT
 void

include/cglm/call/quat.h

@@ -163,7 +163,7 @@ glmc_quat_rotate_atm(mat4 m, versor q, vec3 pivot);
 
 CGLM_EXPORT
 void
-glmc_quat_make(float * __restrict src, versor dest);
+glmc_quat_make(const float * __restrict src, versor dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/ray.h

@@ -20,7 +20,19 @@ glmc_ray_triangle(vec3   origin,
                   vec3   v1,
                   vec3   v2,
                   float *d);
-    
+
+CGLM_EXPORT
+bool
+glmc_ray_sphere(vec3 origin,
+                vec3 dir,
+                vec4 s,
+                float * __restrict t1,
+                float * __restrict t2);
+
+CGLM_EXPORT
+void
+glmc_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec2.h

@@ -195,7 +195,15 @@ glmc_vec2_complex_conjugate(vec2 a, vec2 dest);
 
 CGLM_EXPORT
 void
-glmc_vec2_make(float * __restrict src, vec2 dest);
+glmc_vec2_make(const float * __restrict src, vec2 dest);
+
+CGLM_EXPORT
+void
+glmc_vec2_reflect(vec2 I, vec2 N, vec2 dest);
+
+CGLM_EXPORT
+void
+glmc_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/vec3.h

@@ -332,7 +332,19 @@ glmc_vec3_sqrt(vec3 v, vec3 dest);
 
 CGLM_EXPORT
 void
-glmc_vec3_make(float * __restrict src, vec3 dest);
+glmc_vec3_make(const float * __restrict src, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_reflect(vec3 I, vec3 N, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest);
 
 #ifdef __cplusplus
 }

include/cglm/call/vec4.h

@@ -309,7 +309,15 @@ glmc_vec4_sqrt(vec4 v, vec4 dest);
 
 CGLM_EXPORT
 void
-glmc_vec4_make(float * __restrict src, vec4 dest);
+glmc_vec4_make(const float * __restrict src, vec4 dest);
+
+CGLM_EXPORT
+void
+glmc_vec4_reflect(vec4 I, vec4 N, vec4 dest);
+
+CGLM_EXPORT
+void
+glmc_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest);
 
 #ifdef __cplusplus
 }

include/cglm/mat2.h

@@ -354,7 +354,7 @@ glm_mat2_rmc(vec2 r, mat2 m, vec2 c) {
  */
 CGLM_INLINE
 void
-glm_mat2_make(float * __restrict src, mat2 dest) {
+glm_mat2_make(const float * __restrict src, mat2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[1][0] = src[2];

include/cglm/mat2x3.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat2x3_copy(mat2x3 mat, mat2x3 dest);
    CGLM_INLINE void glm_mat2x3_zero(mat2x3 mat);
-   CGLM_INLINE void glm_mat2x3_make(float * __restrict src, mat2x3 dest);
+   CGLM_INLINE void glm_mat2x3_make(const float * __restrict src, mat2x3 dest);
    CGLM_INLINE void glm_mat2x3_mul(mat2x3 m1, mat3x2 m2, mat2 dest);
    CGLM_INLINE void glm_mat2x3_mulv(mat2x3 m, vec3 v, vec2 dest);
    CGLM_INLINE void glm_mat2x3_transpose(mat2x3 m, mat3x2 dest);
@@ -68,7 +68,7 @@ glm_mat2x3_zero(mat2x3 mat) {
  */
 CGLM_INLINE
 void
-glm_mat2x3_make(float * __restrict src, mat2x3 dest) {
+glm_mat2x3_make(const float * __restrict src, mat2x3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat2x4.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat2x4_copy(mat2x4 mat, mat2x4 dest);
    CGLM_INLINE void glm_mat2x4_zero(mat2x4 mat);
-   CGLM_INLINE void glm_mat2x4_make(float * __restrict src, mat2x4 dest);
+   CGLM_INLINE void glm_mat2x4_make(const float * __restrict src, mat2x4 dest);
    CGLM_INLINE void glm_mat2x4_mul(mat2x4 m1, mat4x2 m2, mat2 dest);
    CGLM_INLINE void glm_mat2x4_mulv(mat2x4 m, vec4 v, vec2 dest);
    CGLM_INLINE void glm_mat2x4_transpose(mat2x4 m, mat4x2 dest);
@@ -64,7 +64,7 @@ glm_mat2x4_zero(mat2x4 mat) {
  */
 CGLM_INLINE
 void
-glm_mat2x4_make(float * __restrict src, mat2x4 dest) {
+glm_mat2x4_make(const float * __restrict src, mat2x4 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat3.h

@@ -436,7 +436,7 @@ glm_mat3_rmc(vec3 r, mat3 m, vec3 c) {
  */
 CGLM_INLINE
 void
-glm_mat3_make(float * __restrict src, mat3 dest) {
+glm_mat3_make(const float * __restrict src, mat3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat3x2.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat3x2_copy(mat3x2 mat, mat3x2 dest);
    CGLM_INLINE void glm_mat3x2_zero(mat3x2 mat);
-   CGLM_INLINE void glm_mat3x2_make(float * __restrict src, mat3x2 dest);
+   CGLM_INLINE void glm_mat3x2_make(const float * __restrict src, mat3x2 dest);
    CGLM_INLINE void glm_mat3x2_mul(mat3x2 m1, mat2x3 m2, mat3 dest);
    CGLM_INLINE void glm_mat3x2_mulv(mat3x2 m, vec2 v, vec3 dest);
    CGLM_INLINE void glm_mat3x2_transpose(mat3x2 m, mat2x3 dest);
@@ -69,7 +69,7 @@ glm_mat3x2_zero(mat3x2 mat) {
  */
 CGLM_INLINE
 void
-glm_mat3x2_make(float * __restrict src, mat3x2 dest) {
+glm_mat3x2_make(const float * __restrict src, mat3x2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
 

include/cglm/mat3x4.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat3x4_copy(mat3x4 mat, mat3x4 dest);
    CGLM_INLINE void glm_mat3x4_zero(mat3x4 mat);
-   CGLM_INLINE void glm_mat3x4_make(float * __restrict src, mat3x4 dest);
+   CGLM_INLINE void glm_mat3x4_make(const float * __restrict src, mat3x4 dest);
    CGLM_INLINE void glm_mat3x4_mul(mat3x4 m1, mat4x3 m2, mat3 dest);
    CGLM_INLINE void glm_mat3x4_mulv(mat3x4 m, vec4 v, vec3 dest);
    CGLM_INLINE void glm_mat3x4_transpose(mat3x4 m, mat4x3 dest);
@@ -66,7 +66,7 @@ glm_mat3x4_zero(mat3x4 mat) {
  */
 CGLM_INLINE
 void
-glm_mat3x4_make(float * __restrict src, mat3x4 dest) {
+glm_mat3x4_make(const float * __restrict src, mat3x4 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/mat4.h

@@ -790,7 +790,7 @@ glm_mat4_rmc(vec4 r, mat4 m, vec4 c) {
  */
 CGLM_INLINE
 void
-glm_mat4_make(float * __restrict src, mat4 dest) {
+glm_mat4_make(const float * __restrict src, mat4 dest) {
   dest[0][0] = src[0];   dest[1][0] = src[4];
   dest[0][1] = src[1];   dest[1][1] = src[5];
   dest[0][2] = src[2];   dest[1][2] = src[6];

include/cglm/mat4x2.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat4x2_copy(mat4x2 mat, mat4x2 dest);
    CGLM_INLINE void glm_mat4x2_zero(mat4x2 mat);
-   CGLM_INLINE void glm_mat4x2_make(float * __restrict src, mat4x2 dest);
+   CGLM_INLINE void glm_mat4x2_make(const float * __restrict src, mat4x2 dest);
    CGLM_INLINE void glm_mat4x2_mul(mat4x2 m1, mat2x4 m2, mat4 dest);
    CGLM_INLINE void glm_mat4x2_mulv(mat4x2 m, vec2 v, vec4 dest);
    CGLM_INLINE void glm_mat4x2_transpose(mat4x2 m, mat2x4 dest);
@@ -72,7 +72,7 @@ glm_mat4x2_zero(mat4x2 mat) {
  */
 CGLM_INLINE
 void
-glm_mat4x2_make(float * __restrict src, mat4x2 dest) {
+glm_mat4x2_make(const float * __restrict src, mat4x2 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
 

include/cglm/mat4x3.h

@@ -13,7 +13,7 @@
  Functions:
    CGLM_INLINE void glm_mat4x3_copy(mat4x3 mat, mat4x3 dest);
    CGLM_INLINE void glm_mat4x3_zero(mat4x3 mat);
-   CGLM_INLINE void glm_mat4x3_make(float * __restrict src, mat4x3 dest);
+   CGLM_INLINE void glm_mat4x3_make(const float * __restrict src, mat4x3 dest);
    CGLM_INLINE void glm_mat4x3_mul(mat4x3 m1, mat3x4 m2, mat4 dest);
    CGLM_INLINE void glm_mat4x3_mulv(mat4x3 m, vec3 v, vec4 dest);
    CGLM_INLINE void glm_mat4x3_transpose(mat4x3 m, mat3x4 dest);
@@ -77,7 +77,7 @@ glm_mat4x3_zero(mat4x3 mat) {
  */
 CGLM_INLINE
 void
-glm_mat4x3_make(float * __restrict src, mat4x3 dest) {
+glm_mat4x3_make(const float * __restrict src, mat4x3 dest) {
   dest[0][0] = src[0];
   dest[0][1] = src[1];
   dest[0][2] = src[2];

include/cglm/quat.h

@@ -894,7 +894,7 @@ glm_quat_rotate_atm(mat4 m, versor q, vec3 pivot) {
  */
 CGLM_INLINE
 void
-glm_quat_make(float * __restrict src, versor dest) {
+glm_quat_make(const float * __restrict src, versor dest) {
   dest[0] = src[0]; dest[1] = src[1];
   dest[2] = src[2]; dest[3] = src[3];
 }

include/cglm/ray.h

@@ -7,12 +7,18 @@
 
 /*
  Functions:
-   CGLM_INLINE bool glm_line_triangle_intersect(vec3   origin,
+   CGLM_INLINE bool glm_ray_triangle(vec3   origin,
-                                                vec3   direction,
+                                     vec3   direction,
-                                                vec3   v0,
+                                     vec3   v0,
-                                                vec3   v1,
+                                     vec3   v1,
-                                                vec3   v2,
+                                     vec3   v2,
-                                                float *d);
+                                     float *d);
+ CGLM_INLINE bool glm_ray_sphere(vec3 origin,
+                                 vec3 dir,
+                                 vec4 s,
+                                 float * __restrict t1,
+                                 float * __restrict t2)
+ CGLM_INLINE void glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
 */
 
 #ifndef cglm_ray_h
@@ -31,7 +37,6 @@
  * @param[in, out] d         distance to intersection
  * @return whether there is intersection
  */
-
 CGLM_INLINE
 bool
 glm_ray_triangle(vec3   origin,
@@ -74,4 +79,93 @@ glm_ray_triangle(vec3   origin,
   return dist > epsilon;
 }
 
+/*!
+ * @brief ray sphere intersection
+ *
+ * - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+ * - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+ * - t1 < 0, t2 < 0: no intersection ahead of the ray
+ * - the caller can check if the intersection points (t1 and t2) fall within a
+ *   specific range (for example, tmin < t1, t2 < tmax) to determine if the
+ *   intersections are within a desired segment of the ray
+ *
+ * @param[in]  origin ray origin
+ * @param[out] dir    normalized ray direction
+ * @param[in]  s      sphere  [center.x, center.y, center.z, radii]
+ * @param[in]  t1     near point1 (closer to origin)
+ * @param[in]  t2     far point2 (farther from origin)
+ */
+CGLM_INLINE
+bool 
+glm_ray_sphere(vec3 origin,
+               vec3 dir,
+               vec4 s,
+               float * __restrict t1,
+               float * __restrict t2) {
+  vec3  dp;
+  float r2, ddp, dpp, dscr, q, tmp, _t1, _t2;
+
+  /* ensure dir is normalized */
+  glm_vec3_sub(s, origin, dp);
+
+  ddp = glm_vec3_dot(dir, dp);
+  dpp = glm_vec3_norm2(dp);
+
+  /* compute the remedy term for numerical stability */
+  glm_vec3_mulsubs(dir, ddp, dp); /* dp: remedy term */
+
+  r2   = s[3] * s[3];
+  dscr = r2 - glm_vec3_norm2(dp);
+
+  if (dscr < 0.0f) {
+    /* no intersection */
+    return false;
+  }
+
+  dscr = sqrtf(dscr);
+  q    = (ddp >= 0.0f) ? (ddp + dscr) : (ddp - dscr);
+
+  /*
+     include Press, William H., Saul A. Teukolsky,
+     William T. Vetterling, and Brian P. Flannery,
+     "Numerical Recipes in C," Cambridge University Press, 1992.
+   */
+  _t1 = q;
+  _t2 = (dpp - r2) / q;
+
+  /* adjust t1 and t2 to ensure t1 is the closer intersection */
+  if (_t1 > _t2) {
+    tmp = _t1;
+    _t1 = _t2;
+    _t2 = tmp;
+  }
+
+  *t1 = _t1;
+  *t2 = _t2;
+
+  /* check if the closest intersection (t1) is behind the ray's origin */
+  if (_t1 < 0.0f && _t2 < 0.0f) {
+    /* both intersections are behind the ray, no visible intersection */
+    return false;
+  }
+
+  return true;
+}
+
+/*!
+ * @brief point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+ *
+ * @param[in]  orig  origin of ray
+ * @param[in]  dir   direction of ray
+ * @param[in]  t     parameter
+ * @param[out] point point at t
+ */
+CGLM_INLINE
+void
+glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point) {
+  vec3 dst;
+  glm_vec3_scale(dir, t, dst);
+  glm_vec3_add(orig, dst, point);
+}
+
 #endif

include/cglm/struct.h

@@ -41,6 +41,7 @@ extern "C" {
 #include "struct/sphere.h"
 #include "struct/curve.h"
 #include "struct/affine2d.h"
+#include "struct/ray.h"
 
 #ifdef __cplusplus
 }

include/cglm/struct/mat2.h

@@ -27,7 +27,7 @@
    CGLM_INLINE void  glms_mat2_swap_col(mat2 mat, int col1, int col2)
    CGLM_INLINE void  glms_mat2_swap_row(mat2 mat, int row1, int row2)
    CGLM_INLINE float glms_mat2_rmc(vec2 r, mat2 m, vec2 c)
-   CGLM_INLINE mat2s glms_mat2_make(float * __restrict src);
+   CGLM_INLINE mat2s glms_mat2_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat2_h
@@ -267,7 +267,7 @@ glms_mat2_(rmc)(vec2s r, mat2s m, vec2s c) {
  */
 CGLM_INLINE
 mat2s
-glms_mat2_(make)(float * __restrict src) {
+glms_mat2_(make)(const float * __restrict src) {
   mat2s r;
   glm_mat2_make(src, r.raw);
   return r;

include/cglm/struct/mat2x3.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat2x3s glms_mat2x3_zero(void);
-   CGLM_INLINE mat2x3s glms_mat2x3_make(float * __restrict src);
+   CGLM_INLINE mat2x3s glms_mat2x3_make(const float * __restrict src);
    CGLM_INLINE mat2s   glms_mat2x3_mul(mat2x3s m1, mat3x2s m2);
    CGLM_INLINE vec2s   glms_mat2x3_mulv(mat2x3s m, vec3s v);
    CGLM_INLINE mat3x2s glms_mat2x3_transpose(mat2x3s m);
@@ -55,7 +55,7 @@ glms_mat2x3_(zero)(void) {
  */
 CGLM_INLINE
 mat2x3s
-glms_mat2x3_(make)(float * __restrict src) {
+glms_mat2x3_(make)(const float * __restrict src) {
   mat2x3s r;
   glm_mat2x3_make(src, r.raw);
   return r;

include/cglm/struct/mat2x4.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat2x4s glms_mat2x4_zero(void);
-   CGLM_INLINE mat2x4s glms_mat2x4_make(float * __restrict src);
+   CGLM_INLINE mat2x4s glms_mat2x4_make(const float * __restrict src);
    CGLM_INLINE mat2s   glms_mat2x4_mul(mat2x4s m1, mat4x2s m2);
    CGLM_INLINE vec2s   glms_mat2x4_mulv(mat2x4s m, vec4s v);
    CGLM_INLINE mat4x2s glms_mat2x4_transpose(mat2x4s m);
@@ -55,7 +55,7 @@ glms_mat2x4_(zero)(void) {
  */
 CGLM_INLINE
 mat2x4s
-glms_mat2x4_(make)(float * __restrict src) {
+glms_mat2x4_(make)(const float * __restrict src) {
   mat2x4s r;
   glm_mat2x4_make(src, r.raw);
   return r;

include/cglm/struct/mat3.h

@@ -28,7 +28,7 @@
    CGLM_INLINE mat3s  glms_mat3_swap_col(mat3s mat, int col1, int col2);
    CGLM_INLINE mat3s  glms_mat3_swap_row(mat3s mat, int row1, int row2);
    CGLM_INLINE float  glms_mat3_rmc(vec3s r, mat3s m, vec3s c);
-   CGLM_INLINE mat3s  glms_mat3_make(float * __restrict src);
+   CGLM_INLINE mat3s  glms_mat3_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat3s_h
@@ -294,7 +294,7 @@ glms_mat3_(rmc)(vec3s r, mat3s m, vec3s c) {
  */
 CGLM_INLINE
 mat3s
-glms_mat3_(make)(float * __restrict src) {
+glms_mat3_(make)(const float * __restrict src) {
   mat3s r;
   glm_mat3_make(src, r.raw);
   return r;

include/cglm/struct/mat3x2.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat3x2s glms_mat3x2_zero(void);
-   CGLM_INLINE mat3x2s glms_mat3x2_make(float * __restrict src);
+   CGLM_INLINE mat3x2s glms_mat3x2_make(const float * __restrict src);
    CGLM_INLINE mat3s   glms_mat3x2_mul(mat3x2s m1, mat2x3s m2);
    CGLM_INLINE vec3s   glms_mat3x2_mulv(mat3x2s m, vec2s v);
    CGLM_INLINE mat2x3s glms_mat3x2_transpose(mat3x2s m);
@@ -55,7 +55,7 @@ glms_mat3x2_(zero)(void) {
  */
 CGLM_INLINE
 mat3x2s
-glms_mat3x2_(make)(float * __restrict src) {
+glms_mat3x2_(make)(const float * __restrict src) {
   mat3x2s r;
   glm_mat3x2_make(src, r.raw);
   return r;

include/cglm/struct/mat3x4.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat3x4s glms_mat3x4_zero(void);
-   CGLM_INLINE mat3x4s glms_mat3x4_make(float * __restrict src);
+   CGLM_INLINE mat3x4s glms_mat3x4_make(const float * __restrict src);
    CGLM_INLINE mat3s   glms_mat3x4_mul(mat3x4s m1, mat4x3s m2);
    CGLM_INLINE vec3s   glms_mat3x4_mulv(mat3x4s m, vec4s v);
    CGLM_INLINE mat4x3s glms_mat3x4_transpose(mat3x4s m);
@@ -55,7 +55,7 @@ glms_mat3x4_(zero)(void) {
  */
 CGLM_INLINE
 mat3x4s
-glms_mat3x4_(make)(float * __restrict src) {
+glms_mat3x4_(make)(const float * __restrict src) {
   mat3x4s r;
   glm_mat3x4_make(src, r.raw);
   return r;

include/cglm/struct/mat4.h

@@ -42,7 +42,7 @@
    CGLM_INLINE mat4s   glms_mat4_swap_col(mat4s mat, int col1, int col2);
    CGLM_INLINE mat4s   glms_mat4_swap_row(mat4s mat, int row1, int row2);
    CGLM_INLINE float   glms_mat4_rmc(vec4s r, mat4s m, vec4s c);
-   CGLM_INLINE mat4s   glms_mat4_make(float * __restrict src);
+   CGLM_INLINE mat4s   glms_mat4_make(const float * __restrict src);
  */
 
 #ifndef cglms_mat4s_h
@@ -468,7 +468,7 @@ glms_mat4_(rmc)(vec4s r, mat4s m, vec4s c) {
  */
 CGLM_INLINE
 mat4s
-glms_mat4_(make)(float * __restrict src) {
+glms_mat4_(make)(const float * __restrict src) {
   mat4s r;
   glm_mat4_make(src, r.raw);
   return r;

include/cglm/struct/mat4x2.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat4x2s glms_mat4x2_zero(void);
-   CGLM_INLINE mat4x2s glms_mat4x2_make(float * __restrict src);
+   CGLM_INLINE mat4x2s glms_mat4x2_make(const float * __restrict src);
    CGLM_INLINE mat4s   glms_mat4x2_mul(mat4x2s m1, mat2x4s m2);
    CGLM_INLINE vec4s   glms_mat4x2_mulv(mat4x2s m, vec2s v);
    CGLM_INLINE mat2x4s glms_mat4x2_transpose(mat4x2s m);
@@ -56,7 +56,7 @@ glms_mat4x2_(zero)(void) {
  */
 CGLM_INLINE
 mat4x2s
-glms_mat4x2_(make)(float * __restrict src) {
+glms_mat4x2_(make)(const float * __restrict src) {
   mat4x2s r;
   glm_mat4x2_make(src, r.raw);
   return r;

include/cglm/struct/mat4x3.h

@@ -12,7 +12,7 @@
 
  Functions:
    CGLM_INLINE mat4x3s glms_mat4x3_zero(void);
-   CGLM_INLINE mat4x3s glms_mat4x3_make(float * __restrict src);
+   CGLM_INLINE mat4x3s glms_mat4x3_make(const float * __restrict src);
    CGLM_INLINE mat4s   glms_mat4x3_mul(mat4x3s m1, mat3x4s m2);
    CGLM_INLINE vec4s   glms_mat4x3_mulv(mat4x3s m, vec3s v);
    CGLM_INLINE mat3x4s glms_mat4x3_transpose(mat4x3s m);
@@ -55,7 +55,7 @@ glms_mat4x3_(zero)(void) {
  */
 CGLM_INLINE
 mat4x3s
-glms_mat4x3_(make)(float * __restrict src) {
+glms_mat4x3_(make)(const float * __restrict src) {
   mat4x3s r;
   glm_mat4x3_make(src, r.raw);
   return r;

include/cglm/struct/quat.h

@@ -574,7 +574,7 @@ glms_quat_(rotate_atm)(versors q, vec3s pivot) {
  */
 CGLM_INLINE
 versors
-glms_quat_(make)(float * __restrict src) {
+glms_quat_(make)(const float * __restrict src) {
   versors dest;
   glm_quat_make(src, dest.raw);
   return dest;

include/cglm/struct/ray.h

@@ -0,0 +1,82 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglms_ray_h
+#define cglms_ray_h
+
+#include "../common.h"
+#include "../types-struct.h"
+#include "../ray.h"
+
+/* api definition */
+#define glms_ray_(NAME) CGLM_STRUCTAPI(ray, NAME)
+
+/*!
+ * @brief Möller–Trumbore ray-triangle intersection algorithm
+ *
+ * @param[in] origin         origin of ray
+ * @param[in] direction      direction of ray
+ * @param[in] v0             first vertex of triangle
+ * @param[in] v1             second vertex of triangle
+ * @param[in] v2             third vertex of triangle
+ * @param[in, out] d         distance to intersection
+ * @return whether there is intersection
+ */
+CGLM_INLINE
+bool
+glms_ray_(triangle)(vec3s  origin,
+                    vec3s  direction,
+                    vec3s  v0,
+                    vec3s  v1,
+                    vec3s  v2,
+                    float *d) {
+  return glm_ray_triangle(origin.raw, direction.raw, v0.raw, v1.raw, v2.raw, d);
+}
+
+/*!
+ * @brief ray sphere intersection
+ *
+ * - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+ * - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+ * - t1 < 0, t2 < 0: no intersection ahead of the ray
+ * - the caller can check if the intersection points (t1 and t2) fall within a
+ *   specific range (for example, tmin < t1, t2 < tmax) to determine if the
+ *   intersections are within a desired segment of the ray
+ *
+ * @param[in]  origin ray origin
+ * @param[out] dir    normalized ray direction
+ * @param[in]  s      sphere  [center.x, center.y, center.z, radii]
+ * @param[in]  t1     near point1 (closer to origin)
+ * @param[in]  t2     far point2 (farther from origin)
+ */
+CGLM_INLINE
+bool
+glms_ray_(sphere)(vec3s origin,
+                  vec3s dir,
+                  vec4s s,
+                  float * __restrict t1,
+                  float * __restrict t2) {
+  return glm_ray_sphere(origin.raw, dir.raw, s.raw, t1, t2);
+}
+
+/*!
+ * @brief point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+ *
+ * @param[in]  orig  origin of ray
+ * @param[in]  dir   direction of ray
+ * @param[in]  t     parameter
+ * @returns point point at t
+ */
+CGLM_INLINE
+vec3s
+glms_ray_(at)(vec3s orig, vec3s dir, float t) {
+  vec3s r;
+  glm_ray_at(orig.raw, orig.raw, t, r.raw);
+  return r;
+}
+
+#endif /* cglms_ray_h */

include/cglm/struct/vec2.h

@@ -54,6 +54,8 @@
    CGLM_INLINE vec2s glms_vec2_clamp(vec2s v, float minVal, float maxVal)
    CGLM_INLINE vec2s glms_vec2_lerp(vec2s from, vec2s to, float t)
    CGLM_INLINE vec2s glms_vec2_make(float * restrict src)
+   CGLM_INLINE vec2s glms_vec2_reflect(vec2s I, vec2s N)
+   CGLM_INLINE vec2s glms_vec2_refract(vec2s I, vec2s N, float eta)
  */
 
 #ifndef cglms_vec2s_h
@@ -685,10 +687,44 @@ glms_vec2_(lerp)(vec2s from, vec2s to, float t) {
  */
 CGLM_INLINE
 vec2s
-glms_vec2_(make)(float * __restrict src) {
+glms_vec2_(make)(const float * __restrict src) {
   vec2s dest;
   glm_vec2_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec2s
+glms_vec2_(reflect)(vec2s I, vec2s N) {
+  vec2s dest;
+  glm_vec2_reflect(I.raw, N.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @param[out] dest refraction result
+ */
+CGLM_INLINE
+vec2s
+glms_vec2_(refract)(vec2s I, vec2s N, float eta) {
+  vec2s dest;
+  glm_vec2_refract(I.raw, N.raw, eta, dest.raw);
+  return dest;
+}
+
 #endif /* cglms_vec2s_h */

include/cglm/struct/vec3.h

@@ -76,6 +76,9 @@
    CGLM_INLINE vec3s glms_vec3_smoothinterpc(vec3s from, vec3s to, float t);
    CGLM_INLINE vec3s glms_vec3_swizzle(vec3s v, int mask);
    CGLM_INLINE vec3s glms_vec3_make(float * restrict src);
+   CGLM_INLINE vec3s glms_vec3_faceforward(vec3s N, vec3s I, vec3s Nref);
+   CGLM_INLINE vec3s glms_vec3_reflect(vec3s I, vec3s N);
+   CGLM_INLINE vec3s glms_vec3_refract(vec3s I, vec3s N, float eta);
 
  Convenient:
    CGLM_INLINE vec3s glms_cross(vec3s a, vec3s b);
@@ -1077,10 +1080,62 @@ glms_vec3_(swizzle)(vec3s v, int mask) {
  */
 CGLM_INLINE
 vec3s
-glms_vec3_(make)(float * __restrict src) {
+glms_vec3_(make)(const float * __restrict src) {
   vec3s dest;
   glm_vec3_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief a vector pointing in the same direction as another
+ *
+ * orients a vector to point away from a surface as defined by its normal
+ *
+ * @param[in] N      vector to orient.
+ * @param[in] I      incident vector
+ * @param[in] Nref   reference vector
+ * @returns oriented vector, pointing away from the surface.
+ */
+CGLM_INLINE
+vec3s
+glms_vec3_(faceforward)(vec3s N, vec3s I, vec3s Nref) {
+  vec3s dest;
+  glm_vec3_faceforward(N.raw, I.raw, Nref.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec3s
+glms_vec3_(reflect)(vec3s I, vec3s N) {
+  vec3s dest;
+  glm_vec3_reflect(I.raw, N.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @returns refraction result
+ */
+CGLM_INLINE
+vec3s
+glms_vec3_(refract)(vec3s I, vec3s N, float eta) {
+  vec3s dest;
+  glm_vec3_refract(I.raw, N.raw, eta, dest.raw);
+  return dest;
+}
+
 #endif /* cglms_vec3s_h */

include/cglm/struct/vec4.h

@@ -67,6 +67,8 @@
    CGLM_INLINE vec4s glms_vec4_cubic(float s);
    CGLM_INLINE vec4s glms_vec4_swizzle(vec4s v, int mask);
    CGLM_INLINE vec4s glms_vec4_make(float * restrict src);
+   CGLM_INLINE vec4s glms_vec4_reflect(vec4s I, vec4s N);
+   CGLM_INLINE vec4s glms_vec4_refract(vec4s I, vec4s N, float eta);
  */
 
 #ifndef cglms_vec4s_h
@@ -921,10 +923,48 @@ glms_vec4_(swizzle)(vec4s v, int mask) {
  */
 CGLM_INLINE
 vec4s
-glms_vec4_(make)(float * __restrict src) {
+glms_vec4_(make)(const float * __restrict src) {
   vec4s dest;
   glm_vec4_make(src, dest.raw);
   return dest;
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @returns reflection result
+ */
+CGLM_INLINE
+vec4s
+glms_vec4_(reflect)(vec4s I, vec4s N) {
+  vec4s dest;
+  glm_vec4_reflect(I.raw, N.raw, dest.raw);
+  return dest;
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * this implementation does not explicitly preserve the 'w' component of the
+ * incident vector 'I' in the output 'dest', users requiring the preservation of
+ * the 'w' component should manually adjust 'dest' after calling this function.
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @returns refraction result
+ */
+CGLM_INLINE
+vec4s
+glms_vec4_(refract)(vec4s I, vec4s N, float eta) {
+  vec4s dest;
+  glm_vec4_refract(I.raw, N.raw, eta, dest.raw);
+  return dest;
+}
+
 #endif /* cglms_vec4s_h */

include/cglm/vec2.h

@@ -55,7 +55,8 @@
    CGLM_INLINE void  glm_vec2_clamp(vec2 v, float minVal, float maxVal)
    CGLM_INLINE void  glm_vec2_lerp(vec2 from, vec2 to, float t, vec2 dest)
    CGLM_INLINE void  glm_vec2_make(float * restrict src, vec2 dest)
-
+   CGLM_INLINE void  glm_vec2_reflect(vec2 I, vec2 N, vec2 dest)
+   CGLM_INLINE void  glm_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest)
  */
 
 #ifndef cglm_vec2_h
@@ -708,8 +709,52 @@ glm_vec2_lerp(vec2 from, vec2 to, float t, vec2 dest) {
  */
 CGLM_INLINE
 void
-glm_vec2_make(float * __restrict src, vec2 dest) {
+glm_vec2_make(const float * __restrict src, vec2 dest) {
   dest[0] = src[0]; dest[1] = src[1];
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[out] dest destination vector for the reflection result
+ */
+CGLM_INLINE
+void
+glm_vec2_reflect(vec2 I, vec2 N, vec2 dest) {
+  vec2 temp;
+  glm_vec2_scale(N, 2.0f * glm_vec2_dot(I, N), temp);
+  glm_vec2_sub(I, temp, dest);
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @param[out] dest refraction result
+ */
+CGLM_INLINE
+void 
+glm_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec2_dot(N, I);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec2_zero(dest);
+    return;
+  }
+
+  glm_vec2_scale(I, eta, dest);
+  glm_vec2_mulsubs(N, eni + sqrtf(k), dest);
+}
+
 #endif /* cglm_vec2_h */

include/cglm/vec3.h

@@ -80,6 +80,9 @@
    CGLM_INLINE void  glm_vec3_smoothinterpc(vec3 from, vec3 to, float t, vec3 dest);
    CGLM_INLINE void  glm_vec3_swizzle(vec3 v, int mask, vec3 dest);
    CGLM_INLINE void  glm_vec3_make(float * restrict src, vec3 dest);
+   CGLM_INLINE void  glm_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest);
+   CGLM_INLINE void  glm_vec3_reflect(vec3 I, vec3 N, vec3 dest);
+   CGLM_INLINE void  glm_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest);
 
  Convenient:
    CGLM_INLINE void  glm_cross(vec3 a, vec3 b, vec3 d);
@@ -1196,10 +1199,76 @@ glm_normalize_to(vec3 v, vec3 dest) {
  */
 CGLM_INLINE
 void
-glm_vec3_make(float * __restrict src, vec3 dest) {
+glm_vec3_make(const float * __restrict src, vec3 dest) {
   dest[0] = src[0];
   dest[1] = src[1];
   dest[2] = src[2];
 }
 
+/*!
+ * @brief a vector pointing in the same direction as another
+ *
+ * orients a vector to point away from a surface as defined by its normal
+ *
+ * @param[in] N      vector to orient
+ * @param[in] I      incident vector
+ * @param[in] Nref   reference vector
+ * @param[out] dest  oriented vector, pointing away from the surface
+ */
+CGLM_INLINE
+void
+glm_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest) {
+  if (glm_vec3_dot(I, Nref) < 0.0f) {
+    /* N is facing away from I */
+    glm_vec3_copy(N, dest);
+  } else {
+    /* N is facing towards I, negate it */
+    glm_vec3_negate_to(N, dest);
+  }
+}
+
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[out] dest reflection result
+ */
+CGLM_INLINE
+void
+glm_vec3_reflect(vec3 I, vec3 N, vec3 dest) {
+  vec3 temp;
+  glm_vec3_scale(N, 2.0f * glm_vec3_dot(I, N), temp);
+  glm_vec3_sub(I, temp, dest);
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @param[out] dest refraction result
+ */
+CGLM_INLINE
+void
+glm_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec3_dot(N, I);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec3_zero(dest);
+    return;
+  }
+
+  glm_vec3_scale(I, eta, dest);
+  glm_vec3_mulsubs(N, eni + sqrtf(k), dest);
+}
+
 #endif /* cglm_vec3_h */

include/cglm/vec4.h

@@ -65,6 +65,8 @@
    CGLM_INLINE void  glm_vec4_smoothinterpc(vec4 from, vec4 to, float t, vec4 dest);
    CGLM_INLINE void  glm_vec4_swizzle(vec4 v, int mask, vec4 dest);
    CGLM_INLINE void  glm_vec4_make(float * restrict src, vec4 dest);
+   CGLM_INLINE void  glm_vec4_reflect(vec4 I, vec4 N, vec4 dest);
+   CGLM_INLINE void  glm_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest);
 
  DEPRECATED:
    glm_vec4_dup
@@ -1299,9 +1301,61 @@ glm_vec4_swizzle(vec4 v, int mask, vec4 dest) {
  */
 CGLM_INLINE
 void
-glm_vec4_make(float * __restrict src, vec4 dest) {
+glm_vec4_make(const float * __restrict src, vec4 dest) {
   dest[0] = src[0]; dest[1] = src[1];
   dest[2] = src[2]; dest[3] = src[3];
 }
 
+/*!
+ * @brief reflection vector using an incident ray and a surface normal
+ *
+ * @param[in]  I    incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[out] dest destination vector for the reflection result
+ */
+CGLM_INLINE
+void
+glm_vec4_reflect(vec4 I, vec4 N, vec4 dest) {
+  vec4 temp;
+
+  /* TODO: direct simd touch */
+  glm_vec4_scale(N, 2.0f * glm_vec4_dot(I, N), temp);
+  glm_vec4_sub(I, temp, dest);
+
+  dest[3] = I[3];
+}
+
+/*!
+ * @brief refraction vector using entering ray, surface normal and refraction index
+ *
+ * if the angle between the entering ray I and the surface normal N is too great
+ * for a given refraction index, the return value is zero
+ *
+ * this implementation does not explicitly preserve the 'w' component of the
+ * incident vector 'I' in the output 'dest', users requiring the preservation of
+ * the 'w' component should manually adjust 'dest' after calling this function.
+ *
+ * @param[in]  I    normalized incident vector
+ * @param[in]  N    normalized normal vector
+ * @param[in]  eta  ratio of indices of refraction
+ * @param[out] dest refraction result
+ */
+CGLM_INLINE
+void 
+glm_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest) {
+  float ndi, eni, k;
+
+  ndi = glm_vec4_dot(N, I);
+  eni = eta * ndi;
+  k   = 1.0f + eta * eta - eni * eni;
+
+  if (k < 0.0f) {
+    glm_vec4_zero(dest);
+    return;
+  }
+
+  glm_vec4_scale(I, eta, dest);
+  glm_vec4_mulsubs(N, eni + sqrtf(k), dest);
+}
+
 #endif /* cglm_vec4_h */

src/mat2.c

@@ -100,6 +100,6 @@ glmc_mat2_rmc(vec2 r, mat2 m, vec2 c) {
 
 CGLM_EXPORT
 void
-glmc_mat2_make(float * __restrict src, mat2 dest) {
+glmc_mat2_make(const float * __restrict src, mat2 dest) {
   glm_mat2_make(src, dest);
 }

src/mat2x3.c

@@ -22,7 +22,7 @@ glmc_mat2x3_zero(mat2x3 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat2x3_make(float * __restrict src, mat2x3 dest) {
+glmc_mat2x3_make(const float * __restrict src, mat2x3 dest) {
   glm_mat2x3_make(src, dest);
 }
 

src/mat2x4.c

@@ -22,7 +22,7 @@ glmc_mat2x4_zero(mat2x4 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat2x4_make(float * __restrict src, mat2x4 dest) {
+glmc_mat2x4_make(const float * __restrict src, mat2x4 dest) {
   glm_mat2x4_make(src, dest);
 }
 

src/mat3.c

@@ -106,6 +106,6 @@ glmc_mat3_rmc(vec3 r, mat3 m, vec3 c) {
 
 CGLM_EXPORT
 void
-glmc_mat3_make(float * __restrict src, mat3 dest) {
+glmc_mat3_make(const float * __restrict src, mat3 dest) {
   glm_mat3_make(src, dest);
 }

src/mat3x2.c

@@ -22,7 +22,7 @@ glmc_mat3x2_zero(mat3x2 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat3x2_make(float * __restrict src, mat3x2 dest) {
+glmc_mat3x2_make(const float * __restrict src, mat3x2 dest) {
   glm_mat3x2_make(src, dest);
 }
 

src/mat3x4.c

@@ -22,7 +22,7 @@ glmc_mat3x4_zero(mat3x4 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat3x4_make(float * __restrict src, mat3x4 dest) {
+glmc_mat3x4_make(const float * __restrict src, mat3x4 dest) {
   glm_mat3x4_make(src, dest);
 }
 

src/mat4.c

@@ -166,6 +166,6 @@ glmc_mat4_rmc(vec4 r, mat4 m, vec4 c) {
 
 CGLM_EXPORT
 void
-glmc_mat4_make(float * __restrict src, mat4 dest) {
+glmc_mat4_make(const float * __restrict src, mat4 dest) {
   glm_mat4_make(src, dest);
 }

src/mat4x2.c

@@ -22,7 +22,7 @@ glmc_mat4x2_zero(mat4x2 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat4x2_make(float * __restrict src, mat4x2 dest) {
+glmc_mat4x2_make(const float * __restrict src, mat4x2 dest) {
   glm_mat4x2_make(src, dest);
 }
 

src/mat4x3.c

@@ -22,7 +22,7 @@ glmc_mat4x3_zero(mat4x3 mat) {
 
 CGLM_EXPORT
 void
-glmc_mat4x3_make(float * __restrict src, mat4x3 dest) {
+glmc_mat4x3_make(const float * __restrict src, mat4x3 dest) {
   glm_mat4x3_make(src, dest);
 }
 

src/quat.c

@@ -232,6 +232,6 @@ glmc_quat_rotate_atm(mat4 m, versor q, vec3 pivot) {
 
 CGLM_EXPORT
 void
-glmc_quat_make(float * __restrict src, versor dest) {
+glmc_quat_make(const float * __restrict src, versor dest) {
   glm_quat_make(src, dest);
 }

src/ray.c

@@ -9,5 +9,21 @@ glmc_ray_triangle(vec3   origin,
                   vec3   v1,
                   vec3   v2,
                   float *d) {
-    return glm_ray_triangle(origin, direction, v0, v1, v2, d);
+  return glm_ray_triangle(origin, direction, v0, v1, v2, d);
+}
+
+CGLM_EXPORT
+bool
+glmc_ray_sphere(vec3 origin,
+                vec3 dir,
+                vec4 s,
+                float * __restrict t1,
+                float * __restrict t2) {
+  return glm_ray_sphere(origin, dir, s, t1, t2);
+}
+
+CGLM_EXPORT
+void
+glmc_ray_at(vec3 orig, vec3 dir, float t, vec3 point) {
+  glm_ray_at(orig, dir, t, point);
 }

src/vec2.c

@@ -299,6 +299,18 @@ glmc_vec2_complex_conjugate(vec2 a, vec2 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec2_make(float * __restrict src, vec2 dest) {
+glmc_vec2_make(const float * __restrict src, vec2 dest) {
   glm_vec2_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec2_reflect(vec2 I, vec2 N, vec2 dest) {
+  glm_vec2_reflect(I, N, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest) {
+  glm_vec2_refract(I, N, eta, dest);
+}

src/vec3.c

@@ -456,6 +456,24 @@ glmc_vec3_sqrt(vec3 v, vec3 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec3_make(float * __restrict src, vec3 dest) {
+glmc_vec3_make(const float * __restrict src, vec3 dest) {
   glm_vec3_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest) {
+  glm_vec3_faceforward(N, I, Nref, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec3_reflect(vec3 I, vec3 N, vec3 dest) {
+  glm_vec3_reflect(I, N, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest) {
+  glm_vec3_refract(I, N, eta, dest);
+}

src/vec4.c

@@ -420,6 +420,18 @@ glmc_vec4_sqrt(vec4 v, vec4 dest) {
 
 CGLM_EXPORT
 void
-glmc_vec4_make(float * __restrict src, vec4 dest) {
+glmc_vec4_make(const float * __restrict src, vec4 dest) {
   glm_vec4_make(src, dest);
 }
+
+CGLM_EXPORT
+void
+glmc_vec4_reflect(vec4 I, vec4 N, vec4 dest) {
+  glm_vec4_reflect(I, N, dest);
+}
+
+CGLM_EXPORT
+void
+glmc_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest) {
+  glm_vec4_refract(I, N, eta, dest);
+}

test/src/test_ray.h

@@ -32,3 +32,43 @@ TEST_IMPL(GLM_PREFIX, ray_triangle) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, ray_sphere) {
+  vec4  sphere = {5.0f, 0.0f, 0.0f, 1.0f}; /* Sphere: center at (5, 0, 0) with radius 1 */
+  float t1, t2;
+  bool  hit;
+
+  /* Case 1: Ray misses the sphere */
+  hit = GLM(ray_sphere)((vec3){10.0f, 3.0f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(!hit); /* Expect no intersection */
+
+  /* Case 2: Ray starts inside the sphere */
+  hit = GLM(ray_sphere)((vec3){5.0f, 0.5f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(hit); /* Expect an intersection */
+  ASSERT(t1 < 0 && t2 > 0); /* Ray exits at t2 */
+
+  /* Case 3: Ray intersects the sphere from outside */
+  hit = GLM(ray_sphere)((vec3){0.0f, 0.0f, 0.0f}, (vec3){1.0f, 0.0f, 0.0f}, sphere, &t1, &t2);
+  ASSERT(hit); /* Expect an intersection */
+  ASSERT(t1 > 0 && t2 > 0); /* Intersections at t1 and t2 */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, ray_at) {
+  vec3 origin    = {0.0f, 0.0f, 0.0f};
+  vec3 direction = {1.0f, 1.0f, 1.0f};  /* Diagonal direction */
+  float distance = sqrtf(3.0f);         /* Distance along the ray; sqrt(3) for unit length due to direction normalization */
+  vec3 result;
+
+  /* Normalize the direction to ensure accurate distance measurement */
+  glm_vec3_normalize(direction);
+
+  GLM(ray_at)(origin, direction, distance, result);
+  ASSERT(fabsf(result[0] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the x-axis */
+  ASSERT(fabsf(result[1] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the y-axis */
+  ASSERT(fabsf(result[2] - 1.0f) <= 0.0000009);  /* Expecting to be 1 unit along the z-axis */
+
+  TEST_SUCCESS
+}
+

test/src/test_vec2.h

@@ -752,3 +752,70 @@ TEST_IMPL(GLM_PREFIX, vec2_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec2_reflect) {
+  vec2 dest;
+
+  /* Reflecting off a "horizontal" surface in 2D */
+  vec2 I1 = {1.0f, -1.0f}; /* Incoming vector */
+  vec2 N1 = {0.0f, 1.0f};  /* Normal vector */
+  GLM(vec2_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f &&
+         fabsf(dest[1] - 1.0f) < 0.00001f); /* Expect reflection upwards */
+
+  /* Reflecting at an angle in 2D */
+  vec2 I2 = {sqrtf(2)/2, -sqrtf(2)/2}; /* Incoming vector at 45 degrees */
+  vec2 N2 = {0.0f, 1.0f}; /* Upwards normal vector */
+  GLM(vec2_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[1] - sqrtf(2)/2) < 0.00001f); /* Expect reflection upwards */
+
+  /* Reflecting off a line in 2D representing a "vertical" surface analogy */
+  vec2 I3 = {1.0f, 0.0f};  /* Incoming vector */
+  vec2 N3 = {-1.0f, 0.0f}; /* Normal vector representing a "vertical" line */
+  GLM(vec2_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f &&
+         fabsf(dest[1]) < 0.00001f); /* Expect reflection to the left */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec2_refract) {
+  vec2 I = {sqrtf(0.5f), -sqrtf(0.5f)}; /* Incoming vector at 45 degrees to normal */
+  vec2 N = {0.0f, 1.0f};                /* Surface normal */
+  vec2 dest;
+  float eta;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  GLM(vec2_refract)(I, N, eta, dest);
+  // In 2D, we expect a similar bending behavior as in 3D, so we check dest[1]
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f)); // Refracted ray bends away from the normal
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f); // Total internal reflection
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  GLM(vec2_refract)(I, N, eta, dest);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  GLM(vec2_refract)(I, N, eta, dest);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal, less bending than air to glass
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  GLM(vec2_refract)(I, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f); // Total internal reflection
+  }
+
+  TEST_SUCCESS
+}
+

test/src/test_vec3.h

@@ -1840,3 +1840,89 @@ TEST_IMPL(GLM_PREFIX, vec3_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec3_faceforward) {
+  vec3 N = {0.0f, 1.0f, 0.0f};
+  vec3 I = {1.0f, -1.0f, 0.0f};
+  vec3 Nref = {0.0f, -1.0f, 0.0f};
+  vec3 dest;
+
+  GLM(vec3_faceforward)(N, I, Nref, dest);
+  ASSERT(dest[0] == 0.0f 
+         && dest[1] == -1.0f
+         && dest[2] == 0.0f); /* Expect N flipped */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec3_reflect) {
+  vec3 dest;
+
+  /* Original test: Reflecting off a horizontal surface */
+  vec3 I1 = {1.0f, -1.0f, 0.0f}; /* Incoming vector */
+  vec3 N1 = {0.0f, 1.0f, 0.0f};  /* Normal vector   */
+  GLM(vec3_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f
+         && fabsf(dest[1] - 1.0f) < 0.00001f
+         && fabsf(dest[2] - 0.0f) < 0.00001f); /* Expect reflection */
+
+  /* Scenario 2: Reflecting off a vertical surface  */
+  vec3 I2 = {1.0f, 0.0f, 0.0f};  /* Incoming vector */
+  vec3 N2 = {-1.0f, 0.0f, 0.0f}; /* Normal vector   */
+  GLM(vec3_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f
+         && fabsf(dest[1]) < 0.00001f
+         && fabsf(dest[2]) < 0.00001f); /* Expect reflection to the left */
+
+  /* Scenario 3: Reflecting at an angle */
+  vec3 I3 = {sqrtf(2)/2, -sqrtf(2)/2, 0.0f}; /* Incoming vector at 45 degrees */
+  vec3 N3 = {0.0f, 1.0f, 0.0f}; /* Upwards normal vector */
+  GLM(vec3_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f
+         && fabsf(dest[1] - sqrtf(2)/2) < 0.00001f
+         && fabsf(dest[2]) < 0.00001f); /* Expect reflection upwards */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec3_refract) {
+  vec3 I = {sqrtf(0.5f), -sqrtf(0.5f), 0.0f}; /* Incoming vector at 45 degrees to normal */
+  vec3 N = {0.0f, 1.0f, 0.0f};                /* Surface normal */
+  vec3 dest;
+  float eta;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  GLM(vec3_refract)(I, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f));
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f);
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  GLM(vec3_refract)(I, N, eta, dest);
+
+  /* Expect bending towards the normal */
+  ASSERT(dest[1] < -sqrtf(0.5f));
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  GLM(vec3_refract)(I, N, eta, dest);
+
+  /* Expect bending towards the normal, less bending than air to glass */
+  ASSERT(dest[1] < -sqrtf(0.5f));
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  GLM(vec3_refract)(I, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f);
+  }
+
+  TEST_SUCCESS
+}

test/src/test_vec4.h

@@ -1536,3 +1536,75 @@ TEST_IMPL(GLM_PREFIX, vec4_make) {
 
   TEST_SUCCESS
 }
+
+TEST_IMPL(GLM_PREFIX, vec4_reflect) {
+  vec4 dest;
+
+  /* Original test: Reflecting off a horizontal surface */
+  vec4 I1 = {1.0f, -1.0f, 0.0f, 0.0f}; /* Incoming vector */
+  vec4 N1 = {0.0f, 1.0f, 0.0f, 0.0f};  /* Normal vector   */
+  GLM(vec4_reflect)(I1, N1, dest);
+  ASSERT(fabsf(dest[0] - 1.0f) < 0.00001f &&
+         fabsf(dest[1] - 1.0f) < 0.00001f &&
+         fabsf(dest[2] - 0.0f) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection */
+
+  /* Scenario 2: Reflecting off a vertical surface */
+  vec4 I2 = {1.0f, 0.0f, 0.0f, 0.0f};  /* Incoming vector */
+  vec4 N2 = {-1.0f, 0.0f, 0.0f, 0.0f}; /* Normal vector   */
+  GLM(vec4_reflect)(I2, N2, dest);
+  ASSERT(fabsf(dest[0] + 1.0f) < 0.00001f &&
+         fabsf(dest[1]) < 0.00001f &&
+         fabsf(dest[2]) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection to the left */
+
+  /* Scenario 3: Reflecting at an angle */
+  vec4 I3 = {sqrtf(2)/2, -sqrtf(2)/2, 0.0f, 0.0f}; /* Incoming vector at 45 degrees */
+  vec4 N3 = {0.0f, 1.0f, 0.0f, 0.0f}; /* Upwards normal vector */
+  GLM(vec4_reflect)(I3, N3, dest);
+  ASSERT(fabsf(dest[0] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[1] - sqrtf(2)/2) < 0.00001f &&
+         fabsf(dest[2]) < 0.00001f &&
+         fabsf(dest[3] - 0.0f) < 0.00001f); /* Expect reflection upwards */
+
+  TEST_SUCCESS
+}
+
+TEST_IMPL(GLM_PREFIX, vec4_refract) {
+  vec4 I = {sqrtf(0.5f), -sqrtf(0.5f), 0.0f, 0.0f}; /* Incoming vector */
+  vec4 N = {0.0f, 1.0f, 0.0f, 0.0f};                /* Surface normal */
+  vec4 dest;
+  float eta;
+
+  /* Water to Air (eta = 1.33/1.0) */
+  eta = 1.33f / 1.0f;
+  GLM(vec4_refract)(I, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f)) {
+    ASSERT(dest[1] < -sqrtf(0.5f));
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f);
+  }
+
+  /* Air to Glass (eta = 1.0 / 1.5) */
+  eta = 1.0f / 1.5f;
+  GLM(vec4_refract)(I, N, eta, dest);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal
+
+  /* Glass to Water (eta = 1.5 / 1.33) */
+  eta = 1.5f / 1.33f;
+  GLM(vec4_refract)(I, N, eta, dest);
+  ASSERT(dest[1] < -sqrtf(0.5f)); // Expect bending towards the normal, less bending than air to glass
+
+  /* Diamond to Air (eta = 2.42 / 1.0) */
+  eta = 2.42f / 1.0f;
+  GLM(vec4_refract)(I, N, eta, dest);
+  if (!(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f)) {
+    /* High potential for total internal reflection, but if it occurs, expect significant bending */
+    ASSERT(dest[1] < -sqrtf(0.5f));
+  } else {
+    ASSERT(dest[0] == 0.0f && dest[1] == 0.0f && dest[2] == 0.0f && dest[3] == 0.0f);
+  }
+
+  TEST_SUCCESS
+}
+

test/tests.h

@@ -392,7 +392,12 @@ TEST_DECLARE(euler)
 
 /* ray */
 TEST_DECLARE(glm_ray_triangle)
+TEST_DECLARE(glm_ray_sphere)
+TEST_DECLARE(glm_ray_at)
+
 TEST_DECLARE(glmc_ray_triangle)
+TEST_DECLARE(glmc_ray_sphere)
+TEST_DECLARE(glmc_ray_at)
 
 /* quat */
 TEST_DECLARE(MACRO_GLM_QUAT_IDENTITY_INIT)
@@ -530,6 +535,8 @@ TEST_DECLARE(glm_vec2_lerp)
 TEST_DECLARE(glm_vec2_complex_mul)
 TEST_DECLARE(glm_vec2_complex_div)
 TEST_DECLARE(glm_vec2_make)
+TEST_DECLARE(glm_vec2_reflect)
+TEST_DECLARE(glm_vec2_refract)
 
 TEST_DECLARE(glmc_vec2)
 TEST_DECLARE(glmc_vec2_copy)
@@ -576,6 +583,8 @@ TEST_DECLARE(glmc_vec2_lerp)
 TEST_DECLARE(glmc_vec2_complex_mul)
 TEST_DECLARE(glmc_vec2_complex_div)
 TEST_DECLARE(glmc_vec2_make)
+TEST_DECLARE(glmc_vec2_reflect)
+TEST_DECLARE(glmc_vec2_refract)
 
 /* vec3 */
 TEST_DECLARE(MACRO_GLM_VEC3_ONE_INIT)
@@ -678,6 +687,9 @@ TEST_DECLARE(glm_vec3_fract)
 TEST_DECLARE(glm_vec3_hadd)
 TEST_DECLARE(glm_vec3_sqrt)
 TEST_DECLARE(glm_vec3_make)
+TEST_DECLARE(glm_vec3_faceforward)
+TEST_DECLARE(glm_vec3_reflect)
+TEST_DECLARE(glm_vec3_refract)
 
 TEST_DECLARE(glmc_vec3)
 TEST_DECLARE(glmc_vec3_copy)
@@ -754,6 +766,9 @@ TEST_DECLARE(glmc_vec3_fract)
 TEST_DECLARE(glmc_vec3_hadd)
 TEST_DECLARE(glmc_vec3_sqrt)
 TEST_DECLARE(glmc_vec3_make)
+TEST_DECLARE(glmc_vec3_faceforward)
+TEST_DECLARE(glmc_vec3_reflect)
+TEST_DECLARE(glmc_vec3_refract)
 
 /* vec4 */
 TEST_DECLARE(MACRO_GLM_VEC4_ONE_INIT)
@@ -842,6 +857,8 @@ TEST_DECLARE(glm_vec4_fract)
 TEST_DECLARE(glm_vec4_hadd)
 TEST_DECLARE(glm_vec4_sqrt)
 TEST_DECLARE(glm_vec4_make)
+TEST_DECLARE(glm_vec4_reflect)
+TEST_DECLARE(glm_vec4_refract)
 
 TEST_DECLARE(glmc_vec4)
 TEST_DECLARE(glmc_vec4_copy3)
@@ -914,6 +931,8 @@ TEST_DECLARE(glmc_vec4_fract)
 TEST_DECLARE(glmc_vec4_hadd)
 TEST_DECLARE(glmc_vec4_sqrt)
 TEST_DECLARE(glmc_vec4_make)
+TEST_DECLARE(glmc_vec4_reflect)
+TEST_DECLARE(glmc_vec4_refract)
 
 /* ivec2 */
 TEST_DECLARE(glm_ivec2)
@@ -1548,10 +1567,15 @@ TEST_LIST {
   TEST_ENTRY(glmc_euler_zyx_quat_lh)
 
   TEST_ENTRY(euler)
-  
+
   /* ray */
   TEST_ENTRY(glm_ray_triangle)
+  TEST_ENTRY(glm_ray_sphere)
+  TEST_ENTRY(glm_ray_at)
+
   TEST_ENTRY(glmc_ray_triangle)
+  TEST_ENTRY(glmc_ray_sphere)
+  TEST_ENTRY(glmc_ray_at)
 
   /* quat */
   TEST_ENTRY(MACRO_GLM_QUAT_IDENTITY_INIT)
@@ -1688,6 +1712,8 @@ TEST_LIST {
   TEST_ENTRY(glm_vec2_complex_mul)
   TEST_ENTRY(glm_vec2_complex_div)
   TEST_ENTRY(glm_vec2_make)
+  TEST_ENTRY(glm_vec2_reflect)
+  TEST_ENTRY(glm_vec2_refract)
 
   TEST_ENTRY(glmc_vec2)
   TEST_ENTRY(glmc_vec2_copy)
@@ -1734,6 +1760,8 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec2_complex_mul)
   TEST_ENTRY(glmc_vec2_complex_div)
   TEST_ENTRY(glmc_vec2_make)
+  TEST_ENTRY(glmc_vec2_reflect)
+  TEST_ENTRY(glmc_vec2_refract)
 
   /* vec3 */
   TEST_ENTRY(MACRO_GLM_VEC3_ONE_INIT)
@@ -1835,6 +1863,9 @@ TEST_LIST {
   TEST_ENTRY(glm_vec3_hadd)
   TEST_ENTRY(glm_vec3_sqrt)
   TEST_ENTRY(glm_vec3_make)
+  TEST_ENTRY(glm_vec3_faceforward)
+  TEST_ENTRY(glm_vec3_reflect)
+  TEST_ENTRY(glm_vec3_refract)
 
   TEST_ENTRY(glmc_vec3)
   TEST_ENTRY(glmc_vec3_copy)
@@ -1911,6 +1942,9 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec3_hadd)
   TEST_ENTRY(glmc_vec3_sqrt)
   TEST_ENTRY(glmc_vec3_make)
+  TEST_ENTRY(glmc_vec3_faceforward)
+  TEST_ENTRY(glmc_vec3_reflect)
+  TEST_ENTRY(glmc_vec3_refract)
 
   /* vec4 */
   TEST_ENTRY(MACRO_GLM_VEC4_ONE_INIT)
@@ -1999,6 +2033,8 @@ TEST_LIST {
   TEST_ENTRY(glm_vec4_hadd)
   TEST_ENTRY(glm_vec4_sqrt)
   TEST_ENTRY(glm_vec4_make)
+  TEST_ENTRY(glm_vec4_reflect)
+  TEST_ENTRY(glm_vec4_refract)
 
   TEST_ENTRY(glmc_vec4)
   TEST_ENTRY(glmc_vec4_copy3)
@@ -2071,6 +2107,8 @@ TEST_LIST {
   TEST_ENTRY(glmc_vec4_hadd)
   TEST_ENTRY(glmc_vec4_sqrt)
   TEST_ENTRY(glmc_vec4_make)
+  TEST_ENTRY(glmc_vec4_reflect)
+  TEST_ENTRY(glmc_vec4_refract)
 
   /* ivec2 */
   TEST_ENTRY(glm_ivec2)

win/cglm.vcxproj

@@ -238,6 +238,9 @@
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_lh.h" />
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_rh.h" />
     <ClInclude Include="..\include\cglm\struct\io.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec2.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec3.h" />
+    <ClInclude Include="..\include\cglm\struct\ivec4.h" />
     <ClInclude Include="..\include\cglm\struct\mat2.h" />
     <ClInclude Include="..\include\cglm\struct\mat2x3.h" />
     <ClInclude Include="..\include\cglm\struct\mat2x4.h" />
@@ -250,6 +253,7 @@
     <ClInclude Include="..\include\cglm\struct\plane.h" />
     <ClInclude Include="..\include\cglm\struct\project.h" />
     <ClInclude Include="..\include\cglm\struct\quat.h" />
+    <ClInclude Include="..\include\cglm\struct\ray.h" />
     <ClInclude Include="..\include\cglm\struct\sphere.h" />
     <ClInclude Include="..\include\cglm\struct\vec2-ext.h" />
     <ClInclude Include="..\include\cglm\struct\vec2.h" />

win/cglm.vcxproj.filters

@@ -720,5 +720,17 @@
     <ClInclude Include="..\include\cglm\struct\handed\euler_to_quat_rh.h">
       <Filter>include\cglm\struct\handed</Filter>
     </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec2.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec3.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ivec4.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\struct\ray.h">
+      <Filter>include\cglm\struct</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>