diff --git a/CREDITS b/CREDITS
index 00cc83235..f875d1a5b 100644
--- a/CREDITS
+++ b/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
diff --git a/Makefile.am b/Makefile.am
index 474569943..2a8ce0eca 100644
--- a/Makefile.am
+++ b/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 \
diff --git a/docs/source/ray.rst b/docs/source/ray.rst
index c5faf336b..ab89b7af6 100644
--- a/docs/source/ray.rst
+++ b/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
diff --git a/docs/source/vec2.rst b/docs/source/vec2.rst
index 102f4fa69..8b1bdf02e 100644
--- a/docs/source/vec2.rst
+++ b/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
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -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
diff --git a/docs/source/vec3.rst b/docs/source/vec3.rst
index c7d92c815..07412c603 100644
--- a/docs/source/vec3.rst
+++ b/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
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -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
diff --git a/docs/source/vec4.rst b/docs/source/vec4.rst
index 4281816f3..df270e75e 100644
--- a/docs/source/vec4.rst
+++ b/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
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -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
diff --git a/include/cglm/call/ray.h b/include/cglm/call/ray.h
index 1fff0554c..e529fdf91 100644
--- a/include/cglm/call/ray.h
+++ b/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
diff --git a/include/cglm/call/vec2.h b/include/cglm/call/vec2.h
index f95166394..bb6e629dc 100644
--- a/include/cglm/call/vec2.h
+++ b/include/cglm/call/vec2.h
@@ -197,6 +197,14 @@ CGLM_EXPORT
 void
 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
 }
 #endif
diff --git a/include/cglm/call/vec3.h b/include/cglm/call/vec3.h
index f3ab653fd..b24529c73 100644
--- a/include/cglm/call/vec3.h
+++ b/include/cglm/call/vec3.h
@@ -334,6 +334,18 @@ CGLM_EXPORT
 void
 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
 }
 #endif
diff --git a/include/cglm/call/vec4.h b/include/cglm/call/vec4.h
index 644facfaf..207bc794e 100644
--- a/include/cglm/call/vec4.h
+++ b/include/cglm/call/vec4.h
@@ -311,6 +311,14 @@ CGLM_EXPORT
 void
 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
 }
 #endif
diff --git a/include/cglm/ray.h b/include/cglm/ray.h
index ced1ad6af..38e171a41 100644
--- a/include/cglm/ray.h
+++ b/include/cglm/ray.h
@@ -7,12 +7,18 @@
 
 /*
  Functions:
-   CGLM_INLINE bool glm_line_triangle_intersect(vec3   origin,
-                                                vec3   direction,
-                                                vec3   v0,
-                                                vec3   v1,
-                                                vec3   v2,
-                                                float *d);
+   CGLM_INLINE bool glm_ray_triangle(vec3   origin,
+                                     vec3   direction,
+                                     vec3   v0,
+                                     vec3   v1,
+                                     vec3   v2,
+                                     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
diff --git a/include/cglm/struct.h b/include/cglm/struct.h
index 030f45a5b..1426589da 100644
--- a/include/cglm/struct.h
+++ b/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
 }
diff --git a/include/cglm/struct/ray.h b/include/cglm/struct/ray.h
new file mode 100644
index 000000000..1341fe97a
--- /dev/null
+++ b/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 */
diff --git a/include/cglm/struct/vec2.h b/include/cglm/struct/vec2.h
index 99db56a64..7a48bce76 100644
--- a/include/cglm/struct/vec2.h
+++ b/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
@@ -691,4 +693,38 @@ glms_vec2_(make)(const float * __restrict src) {
   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 */
diff --git a/include/cglm/struct/vec3.h b/include/cglm/struct/vec3.h
index 2cd208f9a..3be409919 100644
--- a/include/cglm/struct/vec3.h
+++ b/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);
@@ -1083,4 +1086,56 @@ glms_vec3_(make)(const float * __restrict src) {
   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 */
diff --git a/include/cglm/struct/vec4.h b/include/cglm/struct/vec4.h
index c4f200dfd..71e441fb0 100644
--- a/include/cglm/struct/vec4.h
+++ b/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
@@ -927,4 +929,42 @@ glms_vec4_(make)(const float * __restrict src) {
   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 */
diff --git a/include/cglm/vec2.h b/include/cglm/vec2.h
index 5abf61672..58878f551 100644
--- a/include/cglm/vec2.h
+++ b/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
@@ -712,4 +713,48 @@ 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 */
diff --git a/include/cglm/vec3.h b/include/cglm/vec3.h
index a53cd2850..b3dcfeab4 100644
--- a/include/cglm/vec3.h
+++ b/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);
@@ -1202,4 +1205,70 @@ glm_vec3_make(const float * __restrict src, vec3 dest) {
   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 */
diff --git a/include/cglm/vec4.h b/include/cglm/vec4.h
index 8dd5b8578..d8c2ea07e 100644
--- a/include/cglm/vec4.h
+++ b/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
@@ -1304,4 +1306,56 @@ glm_vec4_make(const float * __restrict src, vec4 dest) {
   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 */
diff --git a/src/ray.c b/src/ray.c
index 973c05992..9592b83a5 100644
--- a/src/ray.c
+++ b/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);
 }
diff --git a/src/vec2.c b/src/vec2.c
index 5a23ff57b..27716c4af 100644
--- a/src/vec2.c
+++ b/src/vec2.c
@@ -302,3 +302,15 @@ void
 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);
+}
diff --git a/src/vec3.c b/src/vec3.c
index 86660c83f..33f3bc477 100644
--- a/src/vec3.c
+++ b/src/vec3.c
@@ -459,3 +459,21 @@ void
 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);
+}
diff --git a/src/vec4.c b/src/vec4.c
index 2fcd4ca37..be0cbbce0 100644
--- a/src/vec4.c
+++ b/src/vec4.c
@@ -423,3 +423,15 @@ void
 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);
+}
diff --git a/test/src/test_ray.h b/test/src/test_ray.h
index c1b0281bd..6f2bec508 100644
--- a/test/src/test_ray.h
+++ b/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
+}
+
diff --git a/test/src/test_vec2.h b/test/src/test_vec2.h
index aa75b1c0d..322316532 100644
--- a/test/src/test_vec2.h
+++ b/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
+}
+
diff --git a/test/src/test_vec3.h b/test/src/test_vec3.h
index 65570c02f..6fe4df974 100644
--- a/test/src/test_vec3.h
+++ b/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
+}
diff --git a/test/src/test_vec4.h b/test/src/test_vec4.h
index cf73b67d3..4df471c6a 100644
--- a/test/src/test_vec4.h
+++ b/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
+}
+
diff --git a/test/tests.h b/test/tests.h
index e9509bdef..d85a15596 100644
--- a/test/tests.h
+++ b/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)
diff --git a/win/cglm.vcxproj b/win/cglm.vcxproj
index 949c994f6..dcca926af 100644
--- a/win/cglm.vcxproj
+++ b/win/cglm.vcxproj
@@ -253,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" />
diff --git a/win/cglm.vcxproj.filters b/win/cglm.vcxproj.filters
index 71143f5f7..f200fb29b 100644
--- a/win/cglm.vcxproj.filters
+++ b/win/cglm.vcxproj.filters
@@ -729,5 +729,8 @@
     <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>
\ No newline at end of file