diff --git a/CgeMath.h b/CgeMath.h
index 66a4b68..a4fcd7c 100644
--- a/CgeMath.h
+++ b/CgeMath.h
@@ -126,6 +126,12 @@ int CgeRay2fIntersectBox2f(const float aStart[2], const float aDirection[2],
 int CgeSegment2fIntersectBox2f(const float aStart[2], const float aEnd[2],
                                const float bMin[2], const float bMax[2],
                                float *t, float out[2]);
+int CgeRay2fIntersectCircle(const float start[2], const float direction[2],
+                            const float center[2], float radius, float *t,
+                            float out[2]);
+int CgeSegment2fIntersectCircle(const float start[2], const float end[2],
+                                const float center[2], float radius, float *t,
+                                float out[2]);
 
 int CgeRay3fIntersectPlane(const float start[3], const float direction[3],
                            const float plane[4], float *t, float out[3]);
@@ -143,6 +149,12 @@ int CgeRay3fIntersectBox3f(const float aStart[3], const float aDirection[3],
 int CgeSegment3fIntersectBox3f(const float aStart[3], const float aEnd[3],
                                const float bMin[3], const float bMax[3],
                                float *t, float out[3]);
+int CgeRay3fIntersectSphere(const float start[3], const float direction[3],
+                            const float center[3], float radius, float *t,
+                            float out[3]);
+int CgeSegment3fIntersectSphere(const float start[3], const float end[3],
+                                const float center[3], float radius, float *t,
+                                float out[3]);
 
 void CgeVec2fAdd(const float a[2], const float b[2], float out[2]);
 void CgeVec2fSub(const float a[2], const float b[2], float out[2]);
@@ -162,6 +174,12 @@ void CgeVec2fLerp(const float a[2], const float b[2], float t, float out[2]);
 void CgeVec2fProject(const float a[2], const float b[2], float out[2]);
 void CgeVec2fBarycentric(const float a[2], const float b[2], const float c[2],
                          float v, float w, float out[2]);
+float CgeVec2fDistance(const float a[2], const float b[2]);
+void CgeVec2fReflect(const float in[2], const float norm[2], float out[2]);
+void CgeVec2fRefract(const float in[2], const float norm[2], float eta,
+                     float out[2]);
+void CgeVec2fClamp(const float in[2], float minVal, float maxVal, float out[2]);
+void CgeVec2fAbs(const float in[2], float out[2]);
 
 void CgeVec2iAdd(const int a[2], const int b[2], int out[2]);
 void CgeVec2iSub(const int a[2], const int b[2], int out[2]);
@@ -171,6 +189,8 @@ void CgeVec2iMulAdd(const int a[2], const int b[2], const int c[2], int out[2]);
 void CgeVec2iNegate(const int in[2], int out[2]);
 void CgeVec2iMin(const int a[2], const int b[2], int out[2]);
 void CgeVec2iMax(const int a[2], const int b[2], int out[2]);
+void CgeVec2iClamp(const int in[2], int minVal, int maxVal, int out[2]);
+void CgeVec2iAbs(const int in[2], int out[2]);
 
 void CgeVec3fAdd(const float a[3], const float b[3], float out[3]);
 void CgeVec3fSub(const float a[3], const float b[3], float out[3]);
@@ -190,6 +210,12 @@ void CgeVec3fLerp(const float a[3], const float b[3], float t, float out[3]);
 void CgeVec3fProject(const float a[3], const float b[3], float out[3]);
 void CgeVec3fBarycentric(const float a[3], const float b[3], const float c[3],
                          float v, float w, float out[3]);
+float CgeVec3fDistance(const float a[3], const float b[3]);
+void CgeVec3fReflect(const float in[3], const float norm[3], float out[3]);
+void CgeVec3fRefract(const float in[3], const float norm[3], float eta,
+                     float out[3]);
+void CgeVec3fClamp(const float in[3], float minVal, float maxVal, float out[3]);
+void CgeVec3fAbs(const float in[3], float out[3]);
 
 void CgeVec3iAdd(const int a[3], const int b[3], int out[3]);
 void CgeVec3iSub(const int a[3], const int b[3], int out[3]);
@@ -199,6 +225,8 @@ void CgeVec3iMulAdd(const int a[3], const int b[3], const int c[3], int out[3]);
 void CgeVec3iNegate(const int in[3], int out[3]);
 void CgeVec3iMin(const int a[3], const int b[3], int out[3]);
 void CgeVec3iMax(const int a[3], const int b[3], int out[3]);
+void CgeVec3iClamp(const int in[3], int minVal, int maxVal, int out[3]);
+void CgeVec3iAbs(const int in[3], int out[3]);
 
 void CgeVec4fAdd(const float a[4], const float b[4], float out[4]);
 void CgeVec4fSub(const float a[4], const float b[4], float out[4]);
@@ -217,6 +245,9 @@ void CgeVec4fLerp(const float a[4], const float b[4], float t, float out[4]);
 void CgeVec4fProject(const float a[4], const float b[4], float out[4]);
 void CgeVec4fBarycentric(const float a[4], const float b[4], const float c[4],
                          float v, float w, float out[4]);
+float CgeVec4fDistance(const float a[4], const float b[4]);
+void CgeVec4fClamp(const float in[4], float minVal, float maxVal, float out[4]);
+void CgeVec4fAbs(const float in[4], float out[4]);
 
 void CgeVec4iAdd(const int a[4], const int b[4], int out[4]);
 void CgeVec4iSub(const int a[4], const int b[4], int out[4]);
@@ -226,5 +257,7 @@ void CgeVec4iMulAdd(const int a[4], const int b[4], const int c[4], int out[4]);
 void CgeVec4iNegate(const int in[4], int out[4]);
 void CgeVec4iMin(const int a[4], const int b[4], int out[4]);
 void CgeVec4iMax(const int a[4], const int b[4], int out[4]);
+void CgeVec4iClamp(const int in[4], int minVal, int maxVal, int out[4]);
+void CgeVec4iAbs(const int in[4], int out[4]);
 
 #endif /* CGE_MATH_H */
diff --git a/README.md b/README.md
index 5ef02c9..fc7e5b9 100644
--- a/README.md
+++ b/README.md
@@ -51,6 +51,13 @@ int main() {
 }
 ```
 
+Checkout the `examples` folder to look at an example of the raytracer using
+this library.
+
+Here is an example of its output:
+
+![Example output](examples/demo.png)
+
 ## Portability
 
 - Written in C89 + stdint.h
diff --git a/Ray2f.c b/Ray2f.c
index 1edc295..4c80e26 100644
--- a/Ray2f.c
+++ b/Ray2f.c
@@ -134,13 +134,6 @@ int CgeRay2fIntersectBox2f(const float aStart[2], const float aDirection[2],
     timeNear = -1.0f / 0.0f;
     timeFar  =  1.0f / 0.0f;
 
-    /* Check if origin inside box */
-    if (CgeBox2fContains(bMin, bMax, aStart)) {
-        memcpy(out, aStart, sizeof(float) * 2);
-        *t = 0.0f;
-        return 1;
-    }
-
     /* Check each axis for the minimal and maximum intersection time */
     for (i = 0; i < 2; i++) {
         if (fabsf(aDirection[i]) < EPSILON) {
@@ -191,3 +184,58 @@ int CgeSegment2fIntersectBox2f(const float aStart[2], const float aEnd[2],
     *t = time;
     return 1;
 }
+
+int CgeRay2fIntersectCircle(const float start[2], const float direction[2],
+                            const float center[2], float radius, float *t,
+                            float out[2])
+{
+    float v[2], tmp[2];
+    float tCa, d2, th, t1, t2;
+
+    CgeVec2fSub(center, start, v);
+    tCa = CgeVec2fDot(v, direction);
+    d2 = CgeVec2fDot(v, v) - tCa * tCa;
+
+    if (d2 > radius * radius)
+        return 0;
+
+    th = (float)sqrt(radius * radius - d2);
+    t1 = tCa - th;
+    t2 = tCa + th;
+
+    if (t2 < EPSILON) return 0;
+
+    *t = (t1 > EPSILON) ? t1 : t2;
+    CgeVec2fScale(direction, *t, tmp);
+    CgeVec2fAdd(start, tmp, out);
+
+    return 1;
+}
+
+int CgeSegment2fIntersectCircle(const float start[2], const float end[2],
+                                const float center[2], float radius, float *t,
+                                float out[2]) {
+    float dir[2], lenSq, len;
+    float norm[2];
+
+    CgeVec2fSub(end, start, dir);
+    lenSq = CgeVec2fDot(dir, dir);
+
+    if (lenSq == 0.0f)
+        return 0;
+
+    len = (float)sqrt(lenSq);
+    CgeVec2fScale(dir, 1.0f / len, norm);
+
+    if (!CgeRay2fIntersectCircle(start, norm, center, radius, t, out)) {
+        return 0;
+    }
+
+    *t = *t / len;
+    if (*t > len)
+        return 0;
+
+    CgeVec2fScale(dir, *t, out);
+    CgeVec2fAdd(start, out, out);
+    return 1;
+}
diff --git a/Ray3f.c b/Ray3f.c
index 2279a43..a0dcfd4 100644
--- a/Ray3f.c
+++ b/Ray3f.c
@@ -135,13 +135,6 @@ int CgeRay3fIntersectBox3f(const float aStart[3], const float aDirection[3],
     timeNear = -1.0f / 0.0f;
     timeFar  =  1.0f / 0.0f;
 
-    /* Check if origin inside box */
-    if (CgeBox3fContains(bMin, bMax, aStart)) {
-        memcpy(out, aStart, sizeof(float) * 3);
-        *t = 0.0f;
-        return 1;
-    }
-
     /* Check each axis for the minimal and maximum intersection time */
     for (i = 0; i < 3; i++) {
         if (fabsf(aDirection[i]) < EPSILON) {
@@ -193,3 +186,58 @@ int CgeSegment3fIntersectBox3f(const float aStart[3], const float aEnd[3],
     *t = time;
     return 1;
 }
+
+int CgeRay3fIntersectSphere(const float start[3], const float direction[3],
+                            const float center[3], float radius, float *t,
+                            float out[3])
+{
+    float v[3], tmp[3];
+    float tCa, d2, th, t1, t2;
+
+    CgeVec3fSub(center, start, v);
+    tCa = CgeVec3fDot(v, direction);
+    d2 = CgeVec3fDot(v, v) - tCa * tCa;
+
+    if (d2 > radius * radius)
+        return 0;
+
+    th = (float)sqrt(radius * radius - d2);
+    t1 = tCa - th;
+    t2 = tCa + th;
+
+    if (t2 < EPSILON) return 0;
+
+    *t = (t1 > EPSILON) ? t1 : t2;
+    CgeVec3fScale(direction, *t, tmp);
+    CgeVec3fAdd(start, tmp, out);
+
+    return 1;
+}
+
+int CgeSegment3fIntersectSphere(const float start[3], const float end[3],
+                                const float center[3], float radius, float *t,
+                                float out[3]) {
+    float dir[3], lenSq, len;
+    float norm[3];
+
+    CgeVec3fSub(end, start, dir);
+    lenSq = CgeVec3fDot(dir, dir);
+
+    if (lenSq == 0.0f)
+        return 0;
+
+    len = (float)sqrt(lenSq);
+    CgeVec3fScale(dir, 1.0f / len, norm);
+
+    if (!CgeRay3fIntersectSphere(start, norm, center, radius, t, out)) {
+        return 0;
+    }
+
+    *t = *t / len;
+    if (*t > len)
+        return 0;
+
+    CgeVec3fScale(dir, *t, out);
+    CgeVec3fAdd(start, out, out);
+    return 1;
+}
diff --git a/Vec2f.c b/Vec2f.c
index 4f67b55..64459cc 100644
--- a/Vec2f.c
+++ b/Vec2f.c
@@ -93,3 +93,46 @@ void CgeVec2fBarycentric(const float a[2], const float b[2], const float c[2],
     SET_ROW(tmp1, v); CgeVec2fMulAdd(b, tmp1, tmp2, tmp2);
     SET_ROW(tmp1, w); CgeVec2fMulAdd(c, tmp1, tmp2, out);
 }
+
+float CgeVec2fDistance(const float a[2], const float b[2]) {
+    float tmp[2];
+
+    CgeVec2fSub(a, b, tmp);
+    return CgeVec2fLength(tmp);
+}
+
+void CgeVec2fReflect(const float in[2], const float norm[2], float out[2]) {
+    float tmp[2];
+    float dot;
+
+    dot = CgeVec2fDot(in, norm);
+    CgeVec2fScale(norm, 2.0f * dot, tmp);
+    CgeVec2fSub(in, tmp, out);
+}
+
+void CgeVec2fRefract(const float in[2], const float norm[2], float eta,
+                     float out[2]) {
+    float tmp1[2], tmp2[2];
+    float nDotI, k;
+
+    nDotI = CgeVec2fDot(norm, in);
+    k = 1.0f - eta * eta * (1.0f - nDotI * nDotI);
+    if (k < 0.0f) {
+        out[0] = out[1] = 0.0f;
+        return;
+    }
+
+    CgeVec2fScale(norm, eta * nDotI + (float)sqrt(k), tmp1);
+    CgeVec2fScale(in, eta, tmp2);
+    CgeVec2fSub(tmp2, tmp1, out);
+}
+
+void CgeVec2fClamp(const float in[2], float minVal, float maxVal, float out[2]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+}
+
+void CgeVec2fAbs(const float in[2], float out[2]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+}
diff --git a/Vec2i.c b/Vec2i.c
index 6cad67f..92598ff 100644
--- a/Vec2i.c
+++ b/Vec2i.c
@@ -40,3 +40,13 @@ void CgeVec2iMax(const int a[2], const int b[2], int out[2]) {
     if (a[0] > b[0]) out[0] = a[0]; else out[0] = b[0];
     if (a[1] > b[1]) out[1] = a[1]; else out[1] = b[1];
 }
+
+void CgeVec2iClamp(const int in[2], int minVal, int maxVal, int out[2]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+}
+
+void CgeVec2iAbs(const int in[2], int out[2]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+}
diff --git a/Vec3f.c b/Vec3f.c
index 34d8767..05faa81 100644
--- a/Vec3f.c
+++ b/Vec3f.c
@@ -107,3 +107,48 @@ void CgeVec3fBarycentric(const float a[3], const float b[3], const float c[3],
     SET_ROW(tmp1, v); CgeVec3fMulAdd(b, tmp1, tmp2, tmp2);
     SET_ROW(tmp1, w); CgeVec3fMulAdd(c, tmp1, tmp2, out);
 }
+
+float CgeVec3fDistance(const float a[3], const float b[3]) {
+    float tmp[3];
+
+    CgeVec3fSub(a, b, tmp);
+    return CgeVec3fLength(tmp);
+}
+
+void CgeVec3fReflect(const float in[3], const float norm[3], float out[3]) {
+    float tmp[3];
+    float dot;
+
+    dot = CgeVec3fDot(in, norm);
+    CgeVec3fScale(norm, 2.0f * dot, tmp);
+    CgeVec3fSub(in, tmp, out);
+}
+
+void CgeVec3fRefract(const float in[3], const float norm[3], float eta,
+                     float out[3]) {
+    float tmp1[3], tmp2[3];
+    float nDotI, k;
+
+    nDotI = CgeVec3fDot(norm, in);
+    k = 1.0f - eta * eta * (1.0f - nDotI * nDotI);
+    if (k < 0.0f) {
+        out[0] = out[1] = 0.0f;
+        return;
+    }
+
+    CgeVec3fScale(norm, eta * nDotI + (float)sqrt(k), tmp1);
+    CgeVec3fScale(in, eta, tmp2);
+    CgeVec3fSub(tmp2, tmp1, out);
+}
+
+void CgeVec3fClamp(const float in[3], float minVal, float maxVal, float out[3]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+    out[2] = in[2] < minVal ? minVal : (in[2] > maxVal ? maxVal : in[2]);
+}
+
+void CgeVec3fAbs(const float in[3], float out[3]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+    out[2] = in[2] < 0 ? -in[2] : in[2];
+}
diff --git a/Vec3i.c b/Vec3i.c
index 52532de..5fd8585 100644
--- a/Vec3i.c
+++ b/Vec3i.c
@@ -48,3 +48,15 @@ void CgeVec3iMax(const int a[3], const int b[3], int out[3]) {
     if (a[1] > b[1]) out[1] = a[1]; else out[1] = b[1];
     if (a[2] > b[2]) out[2] = a[2]; else out[2] = b[2];
 }
+
+void CgeVec3iClamp(const int in[3], int minVal, int maxVal, int out[3]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+    out[2] = in[2] < minVal ? minVal : (in[2] > maxVal ? maxVal : in[2]);
+}
+
+void CgeVec3iAbs(const int in[3], int out[3]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+    out[2] = in[2] < 0 ? -in[2] : in[2];
+}
diff --git a/Vec4f.c b/Vec4f.c
index 2d542e8..6da2492 100644
--- a/Vec4f.c
+++ b/Vec4f.c
@@ -106,3 +106,24 @@ void CgeVec4fBarycentric(const float a[4], const float b[4], const float c[4],
     SET_ROW(tmp1, v); CgeVec4fMulAdd(b, tmp1, tmp2, tmp2);
     SET_ROW(tmp1, w); CgeVec4fMulAdd(c, tmp1, tmp2, out);
 }
+
+float CgeVec4fDistance(const float a[4], const float b[4]) {
+    float tmp[4];
+
+    CgeVec4fSub(a, b, tmp);
+    return CgeVec4fLength(tmp);
+}
+
+void CgeVec4fClamp(const float in[4], float minVal, float maxVal, float out[4]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+    out[2] = in[2] < minVal ? minVal : (in[2] > maxVal ? maxVal : in[2]);
+    out[3] = in[3] < minVal ? minVal : (in[3] > maxVal ? maxVal : in[3]);
+}
+
+void CgeVec4fAbs(const float in[4], float out[4]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+    out[2] = in[2] < 0 ? -in[2] : in[2];
+    out[3] = in[3] < 0 ? -in[3] : in[3];
+}
diff --git a/Vec4i.c b/Vec4i.c
index fa9d51a..f85853c 100644
--- a/Vec4i.c
+++ b/Vec4i.c
@@ -56,3 +56,17 @@ void CgeVec4iMax(const int a[4], const int b[4], int out[4]) {
     if (a[2] > b[2]) out[2] = a[2]; else out[2] = b[2];
     if (a[3] > b[3]) out[3] = a[3]; else out[3] = b[3];
 }
+
+void CgeVec4iClamp(const int in[4], int minVal, int maxVal, int out[4]) {
+    out[0] = in[0] < minVal ? minVal : (in[0] > maxVal ? maxVal : in[0]);
+    out[1] = in[1] < minVal ? minVal : (in[1] > maxVal ? maxVal : in[1]);
+    out[2] = in[2] < minVal ? minVal : (in[2] > maxVal ? maxVal : in[2]);
+    out[3] = in[3] < minVal ? minVal : (in[3] > maxVal ? maxVal : in[3]);
+}
+
+void CgeVec4iAbs(const int in[4], int out[4]) {
+    out[0] = in[0] < 0 ? -in[0] : in[0];
+    out[1] = in[1] < 0 ? -in[1] : in[1];
+    out[2] = in[2] < 0 ? -in[2] : in[2];
+    out[3] = in[3] < 0 ? -in[3] : in[3];
+}
diff --git a/examples/demo.c b/examples/demo.c
new file mode 100644
index 0000000..e9a4af3
--- /dev/null
+++ b/examples/demo.c
@@ -0,0 +1,205 @@
+#include <stdio.h>
+#include <math.h>
+#include "CgeMath.h"
+
+#define WIDTH        640
+#define HEIGHT       480
+#define BG_R         0.05f
+#define BG_G         0.05f
+#define BG_B         0.1f
+#define MAX_BOUNCES  7
+
+static int toInt(float x) {
+    return (int)(x * 255.0f);
+}
+
+typedef struct {
+    float center[3];
+    float radius;
+    const float *color;
+    float reflectivity;
+    float shininess;
+    float specular;
+} Sphere;
+
+static const float gRed[3] = {1.0f, 0.0f, 0.0f};
+static const float gGreen[3] = {0.0f, 1.0f, 0.0f};
+static const float gBlue[3] = {0.0f, 0.0f, 1.0f};
+
+static const Sphere gSpheres[] = {
+    {{-2.0f, 0.0f, -5.0f}, 1.0f, gRed, 0.1f, 20.0f, 0.5f},
+    {{ 0.0f, 2.0f, -4.0f}, 1.0f, gGreen, 0.5f, 50.0f, 0.8f},
+    {{ 2.0f, 0.0f, -5.0f}, 1.0f, gBlue, 0.95f, 100.0f, 1.0f}
+};
+
+static const int gNumSpheres = 3;
+static const float gLightPos[3] = {0.0f, 3.0f, -2.0f};
+static const float gAmbient = 0.2f;
+static const float gOrigin[3] = {0.0f, 0.5f, 0.0f};
+
+int traceRay(const float origin[3], const float direction[3], float *t,
+             int *hitIndex, float hitPoint[3], float normal[3]) {
+    float t0, closest;
+    float p[3];
+    int i, hit;
+
+    closest = 1e10f;
+    hit = -1;
+
+    for (i = 0; i < gNumSpheres; i++) {
+        if (CgeRay3fIntersectSphere(origin, direction, gSpheres[i].center, gSpheres[i].radius, &t0, p) == 1) {
+            if (t0 > 0.001f && t0 < closest) {
+                closest = t0;
+                hit = i;
+                CgeVec3fSub(p, gSpheres[i].center, normal);
+                CgeVec3fNormal(normal, normal);
+            }
+        }
+    }
+
+    if (direction[1] != 0.0f) {
+        float t_plane;
+        t_plane = (-1.0f - origin[1]) / direction[1];
+        if (t_plane > 0.001f && t_plane < closest) {
+            closest = t_plane;
+            hit = -2;
+        }
+    }
+
+    if (closest == 1e10f) {
+        return -1;
+    }
+
+    *t = closest;
+    *hitIndex = hit;
+
+    hitPoint[0] = origin[0] + closest * direction[0];
+    hitPoint[1] = origin[1] + closest * direction[1];
+    hitPoint[2] = origin[2] + closest * direction[2];
+
+    if (hit >= 0) {
+        CgeVec3fSub(hitPoint, gSpheres[hit].center, normal);
+        CgeVec3fNormal(normal, normal);
+    } else if (hit == -2) {
+        normal[0] = 0.0f; normal[1] = 1.0f; normal[2] = 0.0f;
+    }
+
+    return 1;
+}
+
+void traceColor(const float origin[3], const float direction[3], float color[3],
+                int depth) {
+    float viewDir[3], halfVec[3], reflectDir[3], offset[3], reflectOrigin[3];
+    float hitPoint[3], normal[3], localColor[3], reflectColor[3];
+    float toLight[3], lightDir[3], dummyPoint[3], dummyNormal[3];
+    float t, NdotL, diffuse, specular, brightness, tile, NdotH, reflectivity;
+    float shadow, tShadow;
+    const float *surf;
+    int hitIndex, dummyIdx, ix, iz;
+
+    if (depth >= MAX_BOUNCES) {
+        color[0] = color[1] = color[2] = 0.0f;
+        return;
+    }
+
+    if (traceRay(origin, direction, &t, &hitIndex, hitPoint, normal) >= 0) {
+        CgeVec3fSub(gLightPos, hitPoint, toLight);
+        CgeVec3fNormal(toLight, lightDir);
+
+        shadow = 1.0f;
+        if (traceRay(hitPoint, lightDir, &tShadow, &dummyIdx, dummyPoint, dummyNormal) >= 0) {
+            float lightDist = CgeVec3fDistance(gLightPos, hitPoint);
+            if (tShadow < lightDist) {
+                shadow = 0.0f;
+            }
+        }
+
+        NdotL = CgeVec3fDot(normal, lightDir);
+        diffuse = (NdotL > 0.0f) ? NdotL : 0.0f;
+
+        specular = 0.0f;
+        if (hitIndex >= 0 && NdotL > 0.0f) {
+            CgeVec3fSub(gOrigin, hitPoint, viewDir);
+            CgeVec3fNormal(viewDir, viewDir);
+
+            CgeVec3fAdd(viewDir, lightDir, halfVec);
+            CgeVec3fNormal(halfVec, halfVec);
+
+            NdotH = CgeVec3fDot(normal, halfVec);
+            if (NdotH > 0.0f) {
+                specular = powf(NdotH, gSpheres[hitIndex].shininess) * gSpheres[hitIndex].specular;
+            }
+        }
+
+        brightness = gAmbient + diffuse * shadow + specular * shadow;
+
+        if (hitIndex == -2) {
+            ix = (int)floorf(hitPoint[0]);
+            iz = (int)floorf(hitPoint[2]);
+            tile = ((ix + iz) % 2 == 0) ? 1.0f : 0.3f;
+            localColor[0] = localColor[1] = localColor[2] = tile * brightness;
+        } else {
+            surf = gSpheres[hitIndex].color;
+            localColor[0] = surf[0] * brightness;
+            localColor[1] = surf[1] * brightness;
+            localColor[2] = surf[2] * brightness;
+        }
+
+        reflectivity = (hitIndex >= 0) ? gSpheres[hitIndex].reflectivity : 0.0f;
+        if (reflectivity > 0.0f) {
+            CgeVec3fReflect(direction, normal, reflectDir);
+            CgeVec3fScale(reflectDir, 0.001f, offset);
+            CgeVec3fAdd(hitPoint, offset, reflectOrigin);
+
+            traceColor(reflectOrigin, reflectDir, reflectColor, depth + 1);
+
+            localColor[0] = localColor[0] * (1.0f - reflectivity) + reflectColor[0] * reflectivity;
+            localColor[1] = localColor[1] * (1.0f - reflectivity) + reflectColor[1] * reflectivity;
+            localColor[2] = localColor[2] * (1.0f - reflectivity) + reflectColor[2] * reflectivity;
+        }
+
+        color[0] = localColor[0];
+        color[1] = localColor[1];
+        color[2] = localColor[2];
+    } else {
+        color[0] = BG_R;
+        color[1] = BG_G;
+        color[2] = BG_B;
+    }
+}
+
+int main() {
+    int x, y;
+    float u, v, aspect, fov;
+    float direction[3], color[3];
+
+    printf("P3\n%d %d\n255\n", WIDTH, HEIGHT);
+
+    aspect = (float)WIDTH / (float)HEIGHT;
+    fov = 1.0f;
+
+    for (y = 0; y < HEIGHT; y++) {
+        for (x = 0; x < WIDTH; x++) {
+            u = (2.0f * (x + 0.5f) / (float)WIDTH - 1.0f) * aspect;
+            v = (2.0f * (HEIGHT - y - 0.5f) / (float)HEIGHT - 1.0f);
+
+            direction[0] = u;
+            direction[1] = v;
+            direction[2] = -fov;
+            CgeVec3fNormal(direction, direction);
+            traceColor(gOrigin, direction, color, 0);
+
+            if (color[0] < 0.0f) color[0] = 0.0f; else if (color[0] > 1.0f) color[0] = 1.0f;
+            if (color[1] < 0.0f) color[1] = 0.0f; else if (color[1] > 1.0f) color[1] = 1.0f;
+            if (color[2] < 0.0f) color[2] = 0.0f; else if (color[2] > 1.0f) color[2] = 1.0f;
+
+            color[0] = powf(color[0], 1.0f / 2.2f);
+            color[1] = powf(color[1], 1.0f / 2.2f);
+            color[2] = powf(color[2], 1.0f / 2.2f);
+
+            printf("%d %d %d\n", toInt(color[0]), toInt(color[1]), toInt(color[2]));
+        }
+    }
+
+    return 0;
+}
diff --git a/examples/demo.png b/examples/demo.png
new file mode 100644
index 0000000..8c8f63b
Binary files /dev/null and b/examples/demo.png differ