diff --git a/.gitattributes b/.gitattributes
index 43118ac..8010e9f 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -1 +1,2 @@
-*.csr eol=lf
\ No newline at end of file
+*.csr eol=lf
+*.sh text eol=lf
\ No newline at end of file
diff --git a/README.md b/README.md
index edf17b5..a27664d 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,10 @@
 <h1 align="center">The Caitlyn Renderer :camera:</h1>
-<p align="center"><img width="600" alt="Render1" src="https://github.com/cypraeno/caitlyn/assets/25397938/9f93e7a7-37d0-43e4-bea1-e81859f75f00"></p>
+<p align="center"><img width="600" alt="Render1" src="https://github.com/user-attachments/assets/a61be07b-fab8-4af2-a34e-0a4faac66713"></p>
 
 
-Caitlyn is an in-development Monte Carlo ray tracer built in C++ by a team of students from the University of Waterloo and Wilfrid Laurier University. Caitlyn is actively working towards being a beautiful tool for bridging the accuracy and lighting of raytracing and the breathtaking visuals of pixel art (see our [portfolio](#our-portfolio) and Odd Tales' "The Last Night"!)
+Caitlyn is an reverse path tracer built for implementing raytracing lighting to 'pixellax' 3D visuals. It is built in C++ by a team of students from the University of Waterloo and Wilfrid Laurier University. 
+
+'Pixellax' is an art style with a focus on using 2D pixel art in 3D environments with an emphasis on beautiful, colourful lighting. (see our [portfolio](#our-portfolio) and Odd Tales' "The Last Night"!)
 
 _Interested in getting involved? Contact [Connor Loi](ctloi@uwaterloo.ca) or [Samuel Bai](sbai@uwaterloo.ca)._
 
@@ -20,24 +22,24 @@ Caitlyn MCRT is built on Debian 12. It may work on other distros, but we recomme
 ### Setup
 Before continuing:
 - Install `Docker Desktop`.
-- Pull the latest `docker pull connortbot/caitlyn-mcrt:base-vX.X.X`
+- Pull the latest `docker pull connortbot/caitlyn-mcrt:base-v0.2.1`
 
 ### Build
 You may pull the repository from within the container or mount a volume. Either works!
-Run `cmake -B build/ -S .` to create files in the `build` folder. `cd build`, and `make`. Don't forget to initialize the submodules.
+Run `cmake -B build/ -S .` to create files in the `build` folder. `cd build`, and `make`. _Don't forget to initialize the submodules!_
 
 ### Basic Rendering
-Caitlyn renders scenes from our custom filetype `.csr`. By default, the `caitlyn` executable will read the scene from a `scene.csr` file, so you need to have one before running. In this guide, we'll just run the `example.csr`, which you can copy from [here](https://github.com/cypraeno/csr-schema/blob/main/examples/example.csr).
+Caitlyn renders scenes from our custom filetype `.csr`. By default, the `caitlyn` executable will read the scene from a `scene.csr` file, so you need to have one before running. In this guide, we'll just run the `0.1.5-showcase.csr`, which you can copy from the `tests` folder.
 
 To learn how to write CSR files, check out the [Basic Guide](https://github.com/cypraeno/csr-schema/blob/main/docs/basic-guide.md).
 
-Caitlyn has a user-friendly command line interface, allowing you to customize samples, depth, type of multithreading, and more. Once you have the executable, you can run `./caitlyn --help` to see all the options at your disposal.
+Caitlyn has a user-friendly command line interface, allowing you to customize samples, depth, use multithreading, and more. Once you have the executable, you can run `./caitlyn --help` to see all the options at your disposal.
 
 Let's render a PNG file of the example scene! Ensure that you have your CSR file in the same directory.
 ```
-./caitlyn -i example.csr -t png -r 600 600
+./caitlyn -i 0.1.5-showcase.csr -o 0.1.5-showcase.png -t png -s 50
 ```
-This will read the scene from `example.csr` and output as a `png`.
+This will read the scene from `0.1.5-showcase.csr` and output as a `png`.
 And now you have your first caitlyn-rendered scene!
 
 ## Our Portfolio
@@ -57,13 +59,10 @@ Flags like `--samples` and `--depth` control the amount of time spent on the ren
 
 Sometimes, CSR files will have features not supported in your version of `caitlyn`. You can check this with the version indicator at the top of the CSR file and with `./caitlyn --version`.
 
-For users who have a better understanding of their computer's resources, the `--threads` and `--vectorization` flags control the use of more efficient architecture. While `threads` dictate the amount of CPU threads to split the workloads on, the `vectorization` flag will dictate the type of SIMD batching. `[NONE|SSE|AVX|AVX512]`.
-
-
 ## Contribute
 For contribution or general inquiries, please email one of us at [Connor Loi](ctloi@uwaterloo.ca) or [Samuel Bai](sbai@uwaterloo.ca).
 
-The people the made it happen:
+Acknowledgements:
 
 <div align="center">
 <a href="https://github.com/connortbot">Connor Loi</a>,
@@ -73,5 +72,16 @@ The people the made it happen:
 <a href="https://github.com/ASharpMarble">Jonathan Wang</a>,
 <a href="https://github.com/18gen">Gen Ichihashi</a>,
 <a href="https://github.com/maxtan84">Max Tan</a>,
-<a href="https://github.com/rickyhuangjh">Ricky Huang</a>,
+<a href="https://github.com/rickyhuangjh">Ricky Huang</a>
 </div>
+
+Citation:
+```bibtex
+@software{Caitlyn,
+    title = {Caitlyn MCRT},
+    author = {Connor Loi and Samuel Bai},
+    note = {https://github.com/cypraeno/caitlyn},
+    version = {0.1.5},
+    year = 2024
+}
+```
diff --git a/csr-schema b/csr-schema
index d60bad3..983fced 160000
--- a/csr-schema
+++ b/csr-schema
@@ -1 +1 @@
-Subproject commit d60bad35643631184f81c7ac037d09952f79aed7
+Subproject commit 983fced58e645cd300abbc69d5a778e1756ea624
diff --git a/include/intersection/hit_info.hh b/include/intersection/hit_info.hh
index e956e2b..ae4a62b 100644
--- a/include/intersection/hit_info.hh
+++ b/include/intersection/hit_info.hh
@@ -6,10 +6,16 @@
 struct HitInfo {
     point3 pos;
     vec3 normal;
+    vec3 microfacet_normal;
     bool front_face;
     float t;
     double u;
     double v;
+    bool transparent = false;
+
+    bool medium = false;
+
+    float rand = 0.5; // used in MixtureBSDF
 
     /** @brief Given a face's outward normal and the initial ray, sets front_face to represent
     if collision hits it from the front or not. */
diff --git a/include/intersection/intersects.h b/include/intersection/intersects.h
index e1a35f1..b825577 100644
--- a/include/intersection/intersects.h
+++ b/include/intersection/intersects.h
@@ -19,4 +19,9 @@ void setupRayHit8(struct RTCRayHit8& rayhit, std::vector<ray>& rays);
 /** @brief modifies given RTCRayHit object to be ready for rtcIntersect16 usage*/
 void setupRayHit16(struct RTCRayHit16& rayhit, std::vector<ray>& rays);
 
+/**
+ * @brief In a given scene, fires a continuous ray and fills information for multiple hits.
+*/
+void MultiIntersect(int max, ray r_in, RTCScene& rtc_scene, std::vector<int>& ids, std::vector<float>& tfars);
+
 #endif
\ No newline at end of file
diff --git a/include/materials/material.h b/include/materials/material.h
index 8197f4b..0a2de23 100644
--- a/include/materials/material.h
+++ b/include/materials/material.h
@@ -4,77 +4,125 @@
 #include "general.h"
 #include "hit_info.hh"
 #include "texture.h"
+#include "onb.h"
+
+#include <complex>
+#include "microfacet.h"
 
 class hit_record;
+class Medium;
+
+// CONSTANTS
+const float SPECULAR_ROUGHNESS_SAMPLING_CUTOFF = 0.1;
+
+enum BSDF_TYPE {
+    DIFFUSE,
+    GLOSSY,
+    SPECULAR,
+    TRANSMISSION,
+    TRANSPARENT
+};
+struct BSDFSample {
+
+    bool scatter;
+    vec3 scatter_direction;
+    color bsdf_value;
+    float pdf_value;
+    BSDF_TYPE type = BSDF_TYPE::DIFFUSE;
+};
 
 class material {
 
     public:
-        virtual color emitted(double u, double v, const point3& p) const {
-            return color(0,0,0);
-        }
+        virtual color emitted(double u, double v, const point3& p) const;
+        virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const;
+        virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+	    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const = 0;
+        virtual BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const;
 };
 
+/**
+ * @class lambertian
+ * @brief Implements basic lambertian material with cosine direction sampling.
+ * @deprecated Use Oren-Nayar for diffuse if possible. At some point, CSR schema should use Diffuse and defualt to Oren-Nayar anyways.
+*/
 class lambertian : public material {
 
     public:
         lambertian(const color& a) : albedo(make_shared<solid_color>(a)) {}
         lambertian(shared_ptr<texture> a) : albedo(a) {}
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override {
-
-            auto scatter_direction = rec.normal + random_unit_vector();
-
-            if (scatter_direction.near_zero()) {
-                scatter_direction = rec.normal;
-            }
+        virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override {
+            onb uvw;
+            uvw.build_from_w(rec.normal);
+            auto scatter_direction = uvw.local(random_cosine_direction());
             scattered = ray(rec.pos, scatter_direction, r_in.time());
-            attenuation = albedo->value(rec.u, rec.v, rec.pos);
             
             return true;
         }
 
+        // A lambertians BRDF value is its albedo / pi
+        virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            return albedo->value(rec.u, rec.v, rec.pos) / pi;
+        }
+        
+        virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            auto cos_theta = dot(rec.normal, scattered.direction().unit_vector());
+            return fmax(0.0, cos_theta / pi);
+        }
+
     private:
     shared_ptr<texture> albedo;
 };
 
-class hemispheric : public material {
 
-    public:
+/**
+ * @class metal
+ * @brief Implements simple coloured perfect specular with fuzz, which is done by randomly warping output direction. Is not physically accurate
+ * since fuzz is not taken into account in f or pdf. Fresnel is not used.
+ * 
+ * @deprecated Use CookTorrance instead. CSR should at some point use Metal, which is CookTorrance and NOT this class.
+ * 
+*/
+class metal : public material {
 
-        hemispheric(const color& a) : albedo(a) {}
+    public:
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override {
-            auto scatter_direction = random_in_hemisphere(rec.normal);
+        metal(const color& a, double f) : albedo(a), fuzz(f < 1 ? f : 1) {}
 
-            if (scatter_direction.near_zero()) {
-                scatter_direction = rec.normal;
-            }
-            scattered = ray(rec.pos,scatter_direction, r_in.time());
-            attenuation = albedo;
+        virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override {
+            vec3 reflected = reflect(r_in.direction().unit_vector(), rec.normal);
+            scattered = ray(rec.pos, reflected + fuzz*random_in_unit_sphere(), r_in.time());
 
-            return true;
+            return (dot(scattered.direction(), rec.normal) > 0);
         }
 
-    public:
-
-        color albedo;
-};
+        virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            return albedo;
+        }
 
-class metal : public material {
+        virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            return 1.0;
+        }
 
-    public:
+        virtual BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override {
+            BSDFSample sample_data;
+            // Sample the microfacet distribution to get the scatter direction.
+            color attenuation; // placeholder until it gets removed from the scatter function header
+            sample_data.scatter = scatter(r_in, rec, attenuation, scattered);
+            sample_data.scatter_direction = scattered.direction().unit_vector();
 
-        metal(const color& a, double f) : albedo(a), fuzz(f < 1 ? f : 1) {}
+            // Sample the BRDF for the value
+            sample_data.bsdf_value = generate(r_in, scattered, rec);
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override {
-            vec3 reflected = reflect(r_in.direction().unit_vector(), rec.normal);
-            scattered = ray(rec.pos, reflected + fuzz*random_in_unit_sphere(), r_in.time());
-            attenuation = albedo;
+            // Find the PDF for the MDF
+            sample_data.pdf_value = pdf(r_in, scattered, rec);
 
-            return (dot(scattered.direction(), rec.normal) > 0);
+            // Provide type for sample
+            if (fuzz < SPECULAR_ROUGHNESS_SAMPLING_CUTOFF) { sample_data.type = BSDF_TYPE::SPECULAR; } // 0.05 was picked arbitrarily, should experiment
+            else { sample_data.type = BSDF_TYPE::GLOSSY; }
+            return sample_data;
         }
 
     public:
@@ -83,14 +131,21 @@ class metal : public material {
         double fuzz;
 };
 
+/**
+ * @class dielectric
+ * @brief Implements simple coloured perfect dielectric (without implementing roughness). Is physically incorrect, does not use fresnel or
+ * any worthwhile techniques.
+ * 
+ * @deprecated Use CookTorranceDielectric instead. CSR should at some point use Transmission, which is CookTorranceDielectric and NOT this class.
+ * 
+*/
 class dielectric : public material {
 
     public:
 
         dielectric(double index_of_refraction) : ir(index_of_refraction) {}
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override {
-            attenuation = color(1.0, 1.0, 1.0);
+        virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override {
             // If the hit is on the front face, ir is the refracted index.
             // If the hit comes from the outside, then 1.0 is the refracted index (air)
             double refraction_ratio = rec.front_face ? (1.0/ir) : ir;
@@ -109,7 +164,34 @@ class dielectric : public material {
             scattered = ray(rec.pos, direction, r_in.time());
             return true;
         }
-        
+
+        virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            return color(1.0, 1.0, 1.0);
+        }
+
+        virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override {
+            return 1.0;
+        }
+
+        virtual BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override {
+            BSDFSample sample_data;
+            // Sample the microfacet distribution to get the scatter direction.
+            color attenuation; // placeholder until it gets removed from the scatter function header
+            sample_data.scatter = scatter(r_in, rec, attenuation, scattered);
+            sample_data.scatter_direction = scattered.direction().unit_vector();
+
+            // Sample the BRDF for the value
+            sample_data.bsdf_value = generate(r_in, scattered, rec);
+
+            // Find the PDF for the MDF
+            sample_data.pdf_value = pdf(r_in, scattered, rec);
+
+            // Provide type for sample
+            sample_data.type = BSDF_TYPE::SPECULAR; 
+            // is incorrect. only works because render function uses direct light sampling + bsdf sampling the same way for specular and transmission.
+
+            return sample_data;
+        }
 
     public:
 
@@ -127,33 +209,226 @@ class dielectric : public material {
         }
 };
 
+/**
+ * @class OrenNayar
+ * @brief Implements the Oren-Nayar reflectance model for simulating the appearance of rough diffuse surfaces.
+ * The Oren-Nayar reflectance model is an extension of the Lambertian model that accounts for the roughness of the
+ * surface, providing a more accurate representation of diffuse reflection from surfaces that are not perfectly smooth.
+ * 
+ * @note The correctness of the pdf and scatter functions, which use cosine-weighted sampling similar to the Lambertian
+ * class, may need further verification.
+ * 
+*/
+class OrenNayar : public material {
+
+    public:
+    OrenNayar(color albedo, float roughness) : albedo{albedo}, roughness{roughness} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override;
+
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+
+    private:
+    color albedo;
+    float roughness;
+};
+
+/**
+ * @class CookTorrance
+ * @brief Implements the Cook-Torrance BRDF model for simulating the specular reflection of a conductor.
+ * This is a more complex model of the original `metal` material. 
+ * This implementation uses the GGX (Trowbridge-Reitz) microfacet distribution to simulate the roughness.
+ * 
+ * @param albedo [OPTIONAL] Base colour.
+ * @param roughness In range [0-1] defines how rough the surface of the material becomes (less shiny).
+ * @param absorption [OPTIONAL] RGB value of how much to not absorb. The higher the color channel, the less that one is absorbed.
+ * @param refraction [OPTIONAL] RGB value of how much to refract (i.e NOT reflect). The higher the color channel, the less it shows.
+ * 
+ * @note by default, if no MDF is specified in the constructor, GGX is used.
+*/
+class CookTorrance : public material {
+
+    public:
+    CookTorrance(color albedo, float roughness)
+        : complex{false}, albedo{albedo}, MDF{std::make_shared<GGX>(roughness)} {}
+
+    CookTorrance(color albedo, std::shared_ptr<Microfacet> mdf)
+        : complex(false), albedo{albedo}, MDF{mdf} {}
+
+    CookTorrance(color absorption, color refraction, float roughness)
+        : complex(true), absorption_coefficient{absorption}, eta{refraction}, MDF{std::make_shared<GGX>(roughness)} {}
+
+    CookTorrance(color absorption, color refraction, std::shared_ptr<Microfacet> mdf)
+        : complex(true), absorption_coefficient{absorption}, eta{refraction}, MDF{mdf} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override;
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+    virtual BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override;
+
+    private:
+    bool complex;
+    color albedo;
+    color absorption_coefficient;
+    color eta;
+    std::shared_ptr<Microfacet> MDF;
+
+    // F, G, D functions
+    vec3 fresnelSchlick(float cosTheta, vec3 F0) const;
+    float FrComplex(float cosTheta_i, std::complex<float> eta) const;
+    vec3 FrComplex(float cosTheta_v, vec3 k, vec3 eta) const;
+
+};
+
+/**
+ * @class CookTorranceDielectric
+ * @brief Implements the Cook-Torrance Dielectric BxDF model.
+ * This implementation uses the GGX (Trowbridge-Reitz) microfacet distribution to simulate the roughness.
+ * 
+ * @param albedo
+ * @param eta Index of refraction (e.g 1.5 for glass)
+ * @param roughness In range [0-1] defines how rough the surface of the material becomes (less shiny).
+ * @param complexFresnel Indicates type of F term to calculate. 0 uses FrComplex, and any positive integer is used as the exponent to the
+ * Schlick approximation.
+ * @note by default, if no MDF is specified in the constructor, GGX is used.
+*/
+class CookTorranceDielectric : public material {
+
+    public:
+    CookTorranceDielectric(color albedo, float eta, float roughness, int complexFresnel = 0) 
+        : albedo{albedo}, eta{(eta == 0.0f) ? 0.0f : (float)fmax(eta, 1.0001f)}, 
+        MDF{std::make_shared<GGX>(roughness)}, complexFresnel{(int)fmax(complexFresnel, 0)} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const;
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+    BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override;
+
+    //private:
+    color albedo;
+    float eta;
+    std::shared_ptr<Microfacet> MDF;
+    int complexFresnel;
+
+    float FrDielectric(float cosTheta_i) const;
+
+    float fresnelSchlick(float cosTheta, int exponent) const;
+
+    private:
+
+    color f_r(const ray& r_in, const HitInfo& rec, const ray& scattered, float R) const;
+    float pdf_r(const ray& r_in, const HitInfo& rec, const ray& scattered, float R) const;
+    float pdf_t(const ray& r_in, const HitInfo& rec, const ray& scattered, float T) const;
+    color f_t(const ray& r_in, const HitInfo& rec, const ray& scattered, float T) const;
+};
+
+class isotropic : public material {
+    public:
+    color albedo;
+
+    isotropic(const color& albedo) : albedo{albedo} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const;
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+};
+
+
 class pixel_lambertian : public material {
 
     public:
         pixel_lambertian(shared_ptr<PixelImageTexture> a) : albedo(a) {}
 
-        virtual bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override {
-            float t = random_double();
-            color4 val = albedo->value(rec.u, rec.v);
-            if (t > val.A) {
-                scattered = ray(rec.pos, r_in.direction(), 0.0);
-                attenuation = color(1.0, 1.0, 1.0);
-            } else {
-                auto scatter_direction = rec.normal + random_unit_vector();
-
-                if (scatter_direction.near_zero()) {
-                    scatter_direction = rec.normal;
-                }
-                scattered = ray(rec.pos, scatter_direction, r_in.time());
-                attenuation = val.RGB;
-            }
-            
-            
-            return true;
-        }
+        virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override;
+        virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+        virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const override;
+        virtual BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const;
 
     private:
     shared_ptr<PixelImageTexture> albedo;
 };
 
+/**
+ * @class MixtureBSDF
+ * @brief A linearly interpolated mixture of N materials, using a vector of weights and materials.
+ * The mixture is done by simple randomization, where the weights decide the frequency that a certain mixed
+ * material is used. This means that there is no confusing interpolation between outgoing directions or sampling.
+ * 
+ * @param weights Vector of floats representing weights of each material. Should sum to 1 to retain energy conservation.
+ * @param mats Vector of materials.
+ * 
+ * @note It is assumed that the given vectors are of the same length and that the weights sum to 1.
+*/
+class MixtureBSDF : public material {
+    public:
+    MixtureBSDF(std::vector<float> weights, std::vector<std::shared_ptr<material>> mats) : weights{weights}, mats{mats} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const;
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const; // assumes that scatter has already been called or sample has already been called, and thus rand is already generated.
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const;
+    BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override; // NOTE: ASSUMES WEIGHTS IS AT LEAST OF SIZE 1 OTHERWISE BEHAVIOUR IS UNDEFINED
+
+    private:
+    std::vector<float> weights;
+    std::vector<std::shared_ptr<material>> mats;
+
+    // Returns negative if weights vector has nothing.
+    int chooseSampleMaterial(float rand) const;
+};
+
+/**
+ * @class LayeredBSDF
+ * @brief A representation of a top BSDF layered over a bottom BSDF with no medium inside. All sampling is done by tracing ray stochastically
+ * through the layers and sampling each time.
+ * 
+ * @param top Top Layer BSDF.
+ * @param bottom Bottom Layer BSDF.
+ * @param termination Russian Roulette termination condition for number of bounces. If exceeded, pretends light is absorbed.
+ * 
+ * @note SCATTER, GENERATE, AND PDF ARE NOT READY.
+ * USE SAMPLE INSTEAD.
+ * 
+ * @bug Bright noise still present
+ * => Replicate with a diffuse overlaid by dielectric (e.g diamond) and notice that
+ * bright noise does not reduce with samples. inf checks do not remove it, its possible that there are
+ * very "large" numbers occurring!
+*/
+class LayeredBSDF : public material {
+    public:
+    LayeredBSDF(std::shared_ptr<material> top, std::shared_ptr<material> bottom, std::shared_ptr<Medium> medium, int termination = 10)
+        : top{top}, bottom{bottom}, medium{medium}, termination{termination} {}
+
+    virtual bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+        return true;
+    }
+
+    virtual color generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+        return color(1,1,1);
+    }
+
+    virtual double pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+        return 1.0;
+    }
+
+    // NOTES:
+    // - The D term in GGX is known to scale at ridiculous amounts to overflow to inf when multiple products, as seen in layering.
+    //   For now, since we know that D exists in the f and pdf, it is safe to arbitrarily set it to 1 or omit it completely, but a better solution is needed.
+    // - Russian Roulette termination does not return black. This is to avoid black specks, but is PHYSICALLY IMPLAUSIBLE.
+    //   There may be a better solution!
+    //   We can check if the ray is bouncing back towards the INITIAL LAYER by if rec.front_face = on_top
+    //   This means that it can never exit via the non-initial layer as a result of Russian Roulette. This is a sacrifice becasue
+    //   there are currently no flags to check if a layer is transmissible or not to exit. However, we know that the initial layer must be.
+    BSDFSample sample(const ray& r_in, HitInfo& rec, ray& scattered) const override;
+
+    private:
+    float thickness = 1.0;
+
+    int termination;
+    std::shared_ptr<material> top;
+    std::shared_ptr<Medium> medium;
+    std::shared_ptr<material> bottom;
+};
+
 #endif
diff --git a/include/materials/medium.h b/include/materials/medium.h
new file mode 100644
index 0000000..3f43fa9
--- /dev/null
+++ b/include/materials/medium.h
@@ -0,0 +1,73 @@
+#ifndef MEDIUM_H
+#define MEDIUM_H
+
+#include <algorithm>
+#include <map>
+#include "general.h"
+#include "material.h"
+
+class Medium {
+    public:
+    Medium(float density, shared_ptr<material> phase_function);
+
+    virtual float particleDistance() const;
+
+    /**
+     * @brief calculates transmittance coefficient between point x and y, assuming a constant medium.
+     * @note Assumes density as equivalent to extinction coefficient.
+     * May not be accurate!
+    */
+    virtual float transmittance(const vec3& x, const vec3& y) const;
+    virtual float transmittance(float dist) const;
+
+    float density;
+    shared_ptr<material> phase;
+};
+
+
+/**
+ * @class MediumRecord
+ * @brief Used for two purposes: tracking of mediums when ray tracing/marching and for calculating transmittance in DLS.
+*/
+class MediumRecord {
+    public:
+    std::vector<shared_ptr<Medium>> mediums;
+    
+    // Transmittance calculation variables
+    float transmittance = 1.0;
+    point3 recent_position;
+    shared_ptr<Medium> highest_density_volume = nullptr;
+
+
+
+    MediumRecord(point3 initial_position);
+    
+    bool contains(std::shared_ptr<Medium> vol) const;
+    
+    /**
+     * @brief Used to update the mediums list and update transmittance.
+     * Called whenever a new volume is hit in which we travel through the participating mediums.
+     * 
+     * 1. Let D be the distance between new point and previous point. This is distance traveled.
+     * 2. Multiply transmittance by the transmittance applied by the highest_density_volume in D distance.
+     * 3A. If volume is in mediums, then we are exiting - remove from mediums.
+     * 3B. If volume is NOT in mediums, then we are entering. Check if vol->density > highest_density_volume->density.
+    */
+    void hitVolume(shared_ptr<Medium> vol, point3 hitPoint);
+
+    /**
+     * @brief Only used for direct light sampling. Given
+    */
+
+    /**
+     * @brief Iterates through mediums and calls particleDistance to simulate the nearest particle of all
+     * participating mediums.
+     * @returns ptr to the particle's corresponding medium and updates passed in dist float.
+    */
+    shared_ptr<Medium> particleDistance(float& dist) const;
+
+    
+};
+
+
+#endif
\ No newline at end of file
diff --git a/include/microfacet.h b/include/microfacet.h
new file mode 100644
index 0000000..b67e24e
--- /dev/null
+++ b/include/microfacet.h
@@ -0,0 +1,72 @@
+#ifndef MICROFACET_H
+#define MICROFACET_H
+
+#include "general.h"
+#include "onb.h"
+
+/**
+ * @class MicrofacetDF
+ * @brief Classes containing definitions to sample and weight for
+ * microfacet distributions in BRDFs.
+*/
+class Microfacet {
+    public:
+    float roughness;
+
+    Microfacet(float roughness) : roughness{roughness} {}
+
+    virtual float D(float cosTheta) const = 0;
+    virtual float G(float cosTheta_V, float cosTheta_L) const = 0;
+    virtual vec3 sample(vec3 normal) const = 0;
+};
+
+class GGX : public Microfacet {
+    public:
+    GGX(float roughness) : Microfacet(roughness) {}
+
+    float D(float NoH) const override {
+        float r = fmax(0.0001, roughness);
+        // float alpha = r * r;
+        // float alpha2 = alpha * alpha;
+        float alpha2 = r * r;
+        float NoH2 = NoH * NoH;
+        float b = (NoH2 * (alpha2 - 1.0) + 1.0);
+        // return (alpha2 / pi) / (b * b);
+        return 1;
+        // - The D term in GGX is known to scale at ridiculous amounts to overflow to inf when multiple products, as seen in layering.
+        //   For now, since we know that D exists in the f and pdf, it is safe to arbitrarily set it to 1 or omit it completely, but a better solution is needed.
+    }
+
+    float G(float NoV, float NoL) const override {
+        return G1_GGX_Schlick(NoL) * G1_GGX_Schlick(NoV);
+    }
+
+    vec3 sample(vec3 N) const override {
+        onb ortho;
+        ortho.build_from_w(N);
+
+        auto e1 = random_double();
+        auto e2 = random_double();
+
+        float theta = atan(roughness * sqrt(e1 / (1 - e1)));
+        float phi = 2 * pi * e2;
+
+        // Calculate normal with y as the up vector
+        float x = sin(theta) * cos(phi);
+        float y = sin(theta) * sin(phi);
+        float z = cos(theta);
+        vec3 microfacet_normal = vec3(x, y, z).unit_vector();
+
+        vec3 adjusted = ortho.local(microfacet_normal);
+        return adjusted;
+    }
+
+    private:
+    float G1_GGX_Schlick(float NoV) const {
+        float alpha = roughness * roughness;
+        float k = alpha / 2.0;
+        return fmax(NoV, 0.001) / (NoV * (1.0 - k) + k);
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/include/objects/geometry.h b/include/objects/geometry.h
index a89b48b..683fe07 100644
--- a/include/objects/geometry.h
+++ b/include/objects/geometry.h
@@ -15,6 +15,9 @@ class Geometry : public Visual {
     RTCGeometry geom;
     Geometry(vec3 position, RTCGeometry geom);
 
+    // Copy constructor, only copies over RTCGeometry field
+    Geometry(shared_ptr<Geometry> geom);
+
     /** @brief given a geomID (referring to ID given by scene attachment), find the material pointer. Usually called by renderer. */
     virtual shared_ptr<material> materialById(unsigned int geomID) const = 0;
 
@@ -23,6 +26,13 @@ class Geometry : public Visual {
      * @return HitInfo A structure containing details about the intersection (e.g., hit point, normal at the hit, material properties).
      */
     virtual HitInfo getHitInfo(const ray& r, const vec3& p, const float t, unsigned int geomID) const = 0;
+
+    virtual point3 sample(const HitInfo& rec) const {
+        return point3(0,0,0);
+    };
+    virtual double pdf(const HitInfo& rec, ray sample_ray) const {
+        return 0.0;
+    };
 };
 
 #endif
\ No newline at end of file
diff --git a/include/objects/light.h b/include/objects/light.h
index fe0cbe3..f3227db 100644
--- a/include/objects/light.h
+++ b/include/objects/light.h
@@ -10,7 +10,7 @@ class emissive : public material {
     public:
     color emission_color;
     emissive(color emission_color);
-    bool scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const override;
+    bool scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const override;
     color emitted(double u, double v, const point3& p) const override;
 };
 
diff --git a/include/objects/primitives/quad_primitive.h b/include/objects/primitives/quad_primitive.h
index abaa171..13998a1 100644
--- a/include/objects/primitives/quad_primitive.h
+++ b/include/objects/primitives/quad_primitive.h
@@ -21,6 +21,8 @@ class QuadPrimitive : public Primitive {
 
         HitInfo getHitInfo(const ray& r, const vec3& p, const float t, unsigned int geomID) const override;
 
+        point3 sample(const HitInfo& rec) const override;
+        double pdf(const HitInfo& rec, ray sample_ray) const override;
         vec3 getV();
         vec3 getU();
 };
diff --git a/include/objects/primitives/sphere_primitive.h b/include/objects/primitives/sphere_primitive.h
index 4ad8448..2861e01 100644
--- a/include/objects/primitives/sphere_primitive.h
+++ b/include/objects/primitives/sphere_primitive.h
@@ -16,8 +16,12 @@ class SpherePrimitive : public Primitive {
 
     HitInfo getHitInfo(const ray& r, const vec3& p, const float t, unsigned int geomID) const override;
 
+    point3 sample(const HitInfo& rec) const override;
+    double pdf(const HitInfo& rec, ray sample_ray) const override;
+
     private:
     static void get_sphere_uv(const point3& p, double& u, double& v);
+    static vec3 random_to_sphere(double radius, double distance_squared);
 };
 
 #endif
\ No newline at end of file
diff --git a/include/objects/volume.h b/include/objects/volume.h
new file mode 100644
index 0000000..a541cb1
--- /dev/null
+++ b/include/objects/volume.h
@@ -0,0 +1,24 @@
+#ifndef VOLUME_H
+#define VOLUME_H
+
+#include "geometry.h"
+#include "medium.h"
+
+class Volume : public Geometry {
+    public:
+    // Releases geometry pointer of the geom ptr, and keeps for itself. Thus, the used geometry cannot be added itself.
+    Volume(shared_ptr<Medium> medium, shared_ptr<Geometry> geom, RTCDevice device);
+
+    virtual shared_ptr<material> materialById(unsigned int geomID) const override;
+
+    virtual HitInfo getHitInfo(const ray& r, const vec3& p, const float t, unsigned int geomID) const override;
+
+    // sample and pdf are not implemented as volumes should not be sampled for direct light sampling
+    // however, reconsider with emissive volumes
+
+    RTCScene volume_scene; // object scene only containing the volume
+    shared_ptr<Medium> medium;
+    shared_ptr<Geometry> geom_ptr;
+};
+
+#endif
\ No newline at end of file
diff --git a/include/onb.h b/include/onb.h
new file mode 100644
index 0000000..13dae3d
--- /dev/null
+++ b/include/onb.h
@@ -0,0 +1,40 @@
+#ifndef ONB_H
+#define ONB_H
+
+#include "general.h"
+
+class onb {
+  public:
+    onb() {}
+
+    vec3 operator[](int i) const { return axis[i]; }
+    vec3& operator[](int i) { return axis[i]; }
+
+    vec3 u() const { return axis[0]; }
+    vec3 v() const { return axis[1]; }
+    vec3 w() const { return axis[2]; }
+
+    vec3 local(double a, double b, double c) const {
+        return a*u() + b*v() + c*w();
+    }
+
+    vec3 local(const vec3& a) const {
+        return a.x()*u() + a.y()*v() + a.z()*w();
+    }
+
+    void build_from_w(const vec3& w) {
+        vec3 unit_w = w.unit_vector();
+        vec3 a = (fabs(unit_w.x()) > 0.9) ? vec3(0,1,0) : vec3(1,0,0);
+        vec3 v = cross(unit_w, a).unit_vector();
+        vec3 u = cross(unit_w, v);
+        axis[0] = u;
+        axis[1] = v;
+        axis[2] = unit_w;
+    }
+
+  public:
+    vec3 axis[3];
+};
+
+
+#endif
\ No newline at end of file
diff --git a/include/render.h b/include/render.h
index 6bc79d5..c9e73f2 100644
--- a/include/render.h
+++ b/include/render.h
@@ -6,6 +6,25 @@
 #include "scene.h"
 #include "vec3.h"
 
+#include "sampling.h"
+
+/**
+ * @brief Most updated integrator for path tracing through scenes
+ * 
+ * 
+ * @bug Known Issues:
+ * - MultiIntersect returns 0 length vector during direct light sampling
+ * - This occurs when sampling FROM a quad, whereby the sampled point of the light yields
+ * a direction PARALLEL to a u or v vector of the quad itself.
+ * e.g if the quad is perpendicular to the y axis and the direction sampled is y=0, then:
+ * -> MultiIntersect returns 0 length vector which normally leads to SEGFAULT.
+ * -> (Unstable?) fix is added, which is to apply a small offset by the normal of the hit progressively until
+ * MultiIntersect succeeds.
+ * 
+ * For now, direct is disabled.
+*/
+color trace_ray(const ray& r, std::shared_ptr<Scene> scene, int depth);
+
 struct RenderData {
     int image_width;
     int image_height;
diff --git a/include/scene.h b/include/scene.h
index 7c496a2..920af0f 100644
--- a/include/scene.h
+++ b/include/scene.h
@@ -5,9 +5,12 @@
 #include <map>
 #include "camera.h"
 #include "material.h"
+#include "light.h"
 #include "sphere_primitive.h"
 #include "instances.h"
+#include "volume.h"
 #include "hit_info.hh"
+#include "color.h"
 
 // SCENE INTERFACE
 // The scene class object covers all relevant objects in a scene:
@@ -30,6 +33,17 @@ class Scene {
     std::map<unsigned int, std::shared_ptr<Geometry>> geom_map;
     RTCScene rtc_scene;
 
+    // Sky Colours
+    color sky_bottom;
+    color sky_top;
+
+    // Lights
+    std::vector<std::shared_ptr<Geometry>> physical_lights;
+    std::vector<std::shared_ptr<Light>> lights;
+
+    // Relevant storage
+    std::vector<std::shared_ptr<Volume>> volumes; // used to check initial mediums
+
     // Default Constructor
     // requires a device to initialize RTCScene
     Scene(RTCDevice device, Camera cam);
@@ -37,7 +51,12 @@ class Scene {
     void commitScene();
     void releaseScene();
     unsigned int add_primitive(std::shared_ptr<Primitive> prim);
+    unsigned int add_volume(std::shared_ptr<Volume> vol);
+
+    void add_physical_light(std::shared_ptr<Geometry> geom_ptr);
     unsigned int add_primitive_instance(std::shared_ptr<PrimitiveInstance> pi_ptr, RTCDevice device);
+
+    void set_sky_colour(color bottom, color top);
 };
 
 void add_sphere(RTCDevice device, RTCScene scene);
diff --git a/include/util/general.h b/include/util/general.h
index 5bd1ea3..96573e6 100644
--- a/include/util/general.h
+++ b/include/util/general.h
@@ -13,6 +13,7 @@ using std::make_shared;
 using std::sqrt;
 
 const double infinity = std::numeric_limits<double>::infinity();
+const double euler = std::exp(1.0);
 const double pi = 3.1415926535897932385;
 
 double degrees_to_radians(double degrees);
diff --git a/include/util/sampling.h b/include/util/sampling.h
new file mode 100644
index 0000000..83de147
--- /dev/null
+++ b/include/util/sampling.h
@@ -0,0 +1,22 @@
+#ifndef SAMPLING_H
+#define SAMPLING_H
+
+// ATTENTION:
+// credit: OSL testrender ©
+// NOT OUR CODE, IN RELEASE, THIS MUST BE REPLACED WITH A PROPER CREDIT!!!
+
+
+struct MIS {
+    // for the function below, enumerate the cases for:
+    // the sampled function being a weight or eval,
+    // the "other" function being a weight or eval
+    enum MISMode { WEIGHT_WEIGHT, WEIGHT_EVAL, EVAL_WEIGHT };
+
+    template<MISMode mode> static inline float power_heuristic(float pdf1, float pdf2) {
+        float num = pow(pdf1, mode);
+        float denom = pow(pdf1, mode) + pow(pdf2, mode);
+        return num / denom;
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/include/util/vec3.h b/include/util/vec3.h
index e8cd975..ca32144 100644
--- a/include/util/vec3.h
+++ b/include/util/vec3.h
@@ -71,6 +71,7 @@ vec3 operator/(vec3 v, float t);
 // vector multiplication
 float dot(const vec3 &u, const vec3 &v);
 vec3 cross(const vec3 &u, const vec3 &v);
+vec3 mix(vec3 x, vec3 y, float a);
 
 /** @brief overloads std::ostream& operator<< to support vec3s */
 std::ostream& operator<<(std::ostream &out, const vec3 &v);
@@ -80,9 +81,14 @@ vec3 random_unit_vector();
 vec3 random_in_unit_sphere();
 vec3 random_in_hemisphere(const vec3& normal);
 vec3 random_in_unit_disk();
+vec3 random_cosine_direction();
 
 // reflection and refraction
+
+// reflects the incoming vector v across the normal n. v is not outward.
 vec3 reflect(const vec3& v, const vec3& n);
+
+// refracts the incoming vector v across the normal n. v is not outward.
 vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat);
 
 #endif
diff --git a/main.cc b/main.cc
index 2e9f84b..e1cdc3b 100644
--- a/main.cc
+++ b/main.cc
@@ -20,4 +20,3 @@ int main(int argc, char* argv[]) {
 
     output(render_data, scene_ptr->cam, scene_ptr, config);
 }
-
diff --git a/src/intersection/intersects.cc b/src/intersection/intersects.cc
index 364f9eb..0dc3b6f 100644
--- a/src/intersection/intersects.cc
+++ b/src/intersection/intersects.cc
@@ -71,3 +71,26 @@ void setupRayHit16(struct RTCRayHit16& rayhit, std::vector<ray>& rays) {
         ix += 1;
     }
 }
+
+void MultiIntersect(int max, ray r_in, RTCScene& rtc_scene, std::vector<int>& ids, std::vector<float>& tfars) {
+    struct RTCRayHit rayhit;
+    ray r = r_in;
+    setupRayHit1(rayhit, r);
+    for (int i=0; i<max; i++) {
+        rtcIntersect1(rtc_scene, &rayhit);
+        int targetID;
+        if (rayhit.hit.instID[0] != RTC_INVALID_GEOMETRY_ID) {
+            targetID = rayhit.hit.instID[0];
+        } else if (rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) {
+            targetID = rayhit.hit.geomID;
+        } else {
+            break;
+        }
+        ids.push_back(targetID);
+        // Calculate tfar for the original ray
+        float time = (r.at(rayhit.ray.tfar) - r_in.origin()).length() / (r_in.direction().length());
+        tfars.push_back(time);
+        r = ray(r.at(rayhit.ray.tfar), r.direction(), 0.0);
+        setupRayHit1(rayhit, r);
+    }
+}
\ No newline at end of file
diff --git a/src/materials/material.cc b/src/materials/material.cc
new file mode 100644
index 0000000..2c5d5ba
--- /dev/null
+++ b/src/materials/material.cc
@@ -0,0 +1,570 @@
+#include "material.h"
+
+#include "medium.h"
+
+color material::emitted(double u, double v, const point3& p) const { return color(0,0,0); }
+bool material::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const { return true; }
+color material::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const { return color(0,0,0); }
+double material::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const { return 1.0; };
+
+// Default behaviour for sample unless overriden by material
+BSDFSample material::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    BSDFSample sample_data;
+    // Sample the microfacet distribution to get the scatter direction.
+    color attenuation; // placeholder until it gets removed from the scatter function header
+    sample_data.scatter = scatter(r_in, rec, attenuation, scattered);
+    sample_data.scatter_direction = scattered.direction().unit_vector();
+
+    // Sample the BRDF for the value
+    sample_data.bsdf_value = generate(r_in, scattered, rec);
+
+    // Find the PDF for the MDF
+    sample_data.pdf_value = pdf(r_in, scattered, rec);
+    return sample_data;
+}
+
+
+// ===================== OREN NAYAR ===============================
+bool OrenNayar::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+    onb uvw;
+    uvw.build_from_w(rec.normal);
+    auto scatter_direction = uvw.local(random_cosine_direction());
+    scattered = ray(rec.pos, scatter_direction, r_in.time());
+    
+    return true;
+}
+
+color OrenNayar::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    vec3 w_i = scattered.direction().unit_vector();
+    vec3 w_o = -(r_in.direction().unit_vector());
+
+    // Calculate azimuthal angles.
+    vec3 projected_i = (w_i - (dot(w_i, rec.normal) * rec.normal)).unit_vector();
+    vec3 projected_o = (w_o - (dot(w_o, rec.normal) * rec.normal)).unit_vector();
+    float cos_azimuth = dot(projected_i, projected_o);
+
+
+    float theta_i = acos(dot(w_i, rec.normal));
+    float theta_o = acos(dot(w_o, rec.normal));
+    
+    float sigma2 = roughness * roughness;
+    float A = 1 - (sigma2 / (2 * (sigma2 + 0.33)));
+
+    float B = (0.45 * sigma2) / (sigma2 + 0.09);
+
+    float alpha = fmax(theta_i, theta_o);
+    float beta = fmin(theta_i, theta_o);
+
+    color diffuse_term = albedo / pi;
+
+
+    return diffuse_term * (A + B * (fmax(0, cos_azimuth) * sin(alpha) * tan(beta)));
+}
+
+double OrenNayar::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    auto cos_theta = dot(rec.normal, scattered.direction().unit_vector());
+    return fmax(0.0, cos_theta / pi);
+}
+
+// ===================== COOK TORRANCE ===============================
+
+bool CookTorrance::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+    vec3 microfacet_normal = MDF->sample(rec.normal);
+    rec.microfacet_normal = microfacet_normal;
+    vec3 scatter_direction = reflect(r_in.direction().unit_vector(), microfacet_normal);
+    scattered = ray(rec.pos, scatter_direction, r_in.time());
+
+    return (dot(scattered.direction(), rec.normal) > 0);
+}
+
+color CookTorrance::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    vec3 L = scattered.direction().unit_vector();
+    vec3 N = rec.normal;
+    vec3 V = -(r_in.direction().unit_vector());
+    vec3 H = (V + L).unit_vector();
+
+    float NoV = clamp(dot(N, V), 0.0, 1.0);
+    float NoL = clamp(dot(N, L), 0.0, 1.0);
+    float NoH = clamp(dot(N, H), 0.0, 1.0);
+    float VoH = clamp(dot(V, H), 0.0, 1.0);
+
+    vec3 f0 = albedo; vec3 F;
+    if (complex) { F = FrComplex(fabs(dot(V,rec.microfacet_normal)), absorption_coefficient, eta); }
+    else { F = fresnelSchlick(VoH, f0); }
+
+    float D = MDF->D(NoH);
+    float G = MDF->G(NoV, NoL);
+
+    vec3 spec = (F * D * G) / (4.0 * fmax(NoV, 0.001) * fmax(NoL, 0.001));
+
+    return spec;
+}
+
+double CookTorrance::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    vec3 V = -r_in.direction().unit_vector();
+    vec3 L = scattered.direction().unit_vector();
+    vec3 H = (V + L).unit_vector();
+    vec3 N = rec.normal;
+
+    float NoH = clamp(dot(N, H), 0.0, 1.0);
+    float VoH = clamp(dot(V, H), 0.0, 1.0);
+
+    float D = MDF->D(NoH);
+    // Convert D(N·H) to pdf based on the microfacet normal distribution.
+    // The Jacobian of the half-vector reflection transformation is |4 * (V·H)|.
+    // This accounts for the change in area density when mapping from H to L.
+    float jacobian = 4.0 * abs(dot(V, H));
+    if (jacobian < 0.0001) return 0;
+
+    return D / jacobian;
+}
+
+BSDFSample CookTorrance::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    BSDFSample sample_data;
+    // Sample the microfacet distribution to get the scatter direction.
+    color attenuation; // placeholder until it gets removed from the scatter function header
+    sample_data.scatter = scatter(r_in, rec, attenuation, scattered);
+    sample_data.scatter_direction = scattered.direction().unit_vector();
+
+    // Sample the BRDF for the value
+    sample_data.bsdf_value = generate(r_in, scattered, rec);
+
+    // Find the PDF for the MDF
+    sample_data.pdf_value = pdf(r_in, scattered, rec);
+    if (MDF->roughness < SPECULAR_ROUGHNESS_SAMPLING_CUTOFF) { sample_data.type = BSDF_TYPE::SPECULAR; } // 0.05 was picked arbitrarily, should experiment
+    else { sample_data.type = BSDF_TYPE::GLOSSY; }
+    return sample_data;
+}
+
+vec3 CookTorrance::fresnelSchlick(float cosTheta, vec3 F0) const {
+    return F0 + (color(1.0, 1.0, 1.0) - F0) * pow(1.0 - cosTheta, 5.0);
+}
+
+float CookTorrance::FrComplex(float cosTheta_i, std::complex<float> eta) const {
+    using Complex = std::complex<float>;
+    cosTheta_i = clamp(cosTheta_i, 0, 1);
+    float sin2Theta_i = 1 - (cosTheta_i * cosTheta_i);
+    Complex sin2Theta_t = sin2Theta_i / (eta * eta);
+    Complex val(1, -2);
+    Complex cosTheta_t = std::sqrt(val - sin2Theta_t);
+    
+    Complex r_parl = (eta * cosTheta_i - cosTheta_t) /
+                    (eta * cosTheta_i + cosTheta_t);
+    Complex r_perp = (cosTheta_i - eta * cosTheta_t) /
+                    (cosTheta_i + eta * cosTheta_t);
+    return (std::norm(r_parl) + std::norm(r_perp)) / 2;
+}
+
+vec3 CookTorrance::FrComplex(float cosTheta_v, vec3 k, vec3 eta) const {
+    float x = FrComplex(cosTheta_v, std::complex<float>(eta.x(), k.x()));
+    float y = FrComplex(cosTheta_v, std::complex<float>(eta.y(), k.y()));
+    float z = FrComplex(cosTheta_v, std::complex<float>(eta.z(), k.z()));
+    return vec3(x,y,z);
+}
+
+// ===================== COOK DIELECTRIC ===============================
+
+bool CookTorranceDielectric::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+    vec3 wo = -r_in.direction().unit_vector();
+    vec3 N = rec.normal;
+    vec3 wm = MDF->sample(N); // outward microfacet normal.
+    rec.microfacet_normal = wm;
+
+    float cosTheta_i = dot(wo, wm);
+    float R;
+    if (complexFresnel == 0) { R = FrDielectric(cosTheta_i); }
+    else { R = fresnelSchlick(cosTheta_i, complexFresnel); }
+    float T = 1 - R;
+
+    float u = random_double();
+    double refraction_ratio = rec.front_face ? (1.0/eta) : eta;
+    double sinTheta_i = sqrt(1.0 - cosTheta_i*cosTheta_i);
+
+    if (u < (R / (R + T)) || refraction_ratio * sinTheta_i > 1.0) { // reflectance
+        vec3 wi = reflect(-wo, wm);
+        scattered = ray(rec.pos, wi, r_in.time());
+    } else {
+        vec3 wi = refract(-wo, wm, refraction_ratio);
+        scattered = ray(rec.pos, wi, r_in.time());
+    }
+    return true;
+}
+color CookTorranceDielectric::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const { // assumes wm has been defined in rec
+    vec3 wo = -r_in.direction().unit_vector();
+    vec3 N = rec.normal;
+    vec3 wi = scattered.direction();
+    vec3 wm = rec.microfacet_normal;
+    float cosTheta_i = dot(wo, wm);
+    float R;
+    if (complexFresnel == 0) { R = FrDielectric(cosTheta_i); }
+    else { R = fresnelSchlick(cosTheta_i, complexFresnel); }
+    float T = 1 - R;
+    if (cosTheta_i > 0) { // reflectance
+        return f_r(r_in, rec, scattered, R);
+    } else { // refractance
+        return f_t(r_in, rec, scattered, T);
+    }
+}
+double CookTorranceDielectric::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const { // assumes wm has been defined in rec
+    vec3 wo = -r_in.direction().unit_vector();
+    vec3 N = rec.normal;
+    vec3 wi = scattered.direction();
+    vec3 wm = rec.microfacet_normal;
+    float cosTheta_i = dot(wo, wm);
+    float R;
+    if (complexFresnel == 0) { R = FrDielectric(cosTheta_i); }
+    else { R = fresnelSchlick(cosTheta_i, complexFresnel); }
+    float T = 1 - R;
+    if (cosTheta_i > 0) { // reflectance
+        return pdf_r(r_in, rec, scattered, R);
+    } else { // refractance
+        return pdf_t(r_in, rec, scattered, T);
+    }
+};
+
+BSDFSample CookTorranceDielectric::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    // Vectors wo and wi are the outgoing and incident directions respectively.
+    vec3 wo = -r_in.direction().unit_vector();
+
+    BSDFSample sample_data;
+    vec3 N = rec.normal;
+    vec3 wm = MDF->sample(N); // outward microfacet normal.
+    rec.microfacet_normal = wm;
+
+    float cosTheta_i = dot(wo, wm);
+    float R;
+    if (complexFresnel == 0) { R = FrDielectric(cosTheta_i); }
+    else { R = fresnelSchlick(cosTheta_i, complexFresnel); }
+    float T = 1 - R;
+
+    float u = random_double();
+    double refraction_ratio = rec.front_face ? (1.0/eta) : eta;
+    double sinTheta_i = sqrt(1.0 - cosTheta_i*cosTheta_i);
+
+    if (u < (R / (R + T)) || refraction_ratio * sinTheta_i > 1.0) { // reflectance
+        vec3 wi = reflect(-wo, wm);
+        scattered = ray(rec.pos, wi, r_in.time());
+        sample_data.scatter_direction = wi;
+        sample_data.scatter = (dot(wi, N) > 0);
+
+        sample_data.bsdf_value = f_r(r_in, rec, scattered, R);
+        sample_data.pdf_value = pdf_r(r_in, rec, scattered, R);
+        if (MDF->roughness < SPECULAR_ROUGHNESS_SAMPLING_CUTOFF) { sample_data.type = BSDF_TYPE::SPECULAR; } // 0.05 was picked arbitrarily, should experiment
+        else { sample_data.type = BSDF_TYPE::GLOSSY; }
+    } else { // transmission
+        vec3 wi = refract(-wo, wm, refraction_ratio);
+        scattered = ray(rec.pos, wi, r_in.time());
+        sample_data.scatter_direction = wi;
+        sample_data.scatter = (dot(wi, N) < 0);
+
+        sample_data.bsdf_value = f_t(r_in, rec, scattered, T);
+        sample_data.pdf_value = pdf_t(r_in, rec, scattered, T);
+        sample_data.type = BSDF_TYPE::TRANSMISSION;
+    }
+
+    return sample_data;
+}
+
+float CookTorranceDielectric::FrDielectric(float cosTheta_i) const {
+    float temp_eta = eta;
+    cosTheta_i = clamp(cosTheta_i, -1.0, 1.0);
+    if (cosTheta_i < 0) {
+        temp_eta = 1 / eta;
+        cosTheta_i = -cosTheta_i;
+    }
+
+    float sin2Theta_i = 1 - (cosTheta_i * cosTheta_i);
+    float sin2Theta_t = sin2Theta_i / (temp_eta * temp_eta);
+    if (sin2Theta_t >= 1.0) {
+        return 1.0;
+    }
+    float cosTheta_t = sqrt(1 - sin2Theta_t);
+    float r_parallel = (temp_eta * cosTheta_i - cosTheta_t) / (temp_eta * cosTheta_i + cosTheta_t);
+    float r_perp = (cosTheta_i - (temp_eta * cosTheta_t)) / (cosTheta_i + (eta * cosTheta_t));
+    return ((r_parallel * r_parallel) + (r_perp * r_perp)) / 2;
+}
+
+float CookTorranceDielectric::fresnelSchlick(float cosTheta, int exponent) const {
+    float F0 = pow(((1 - eta) / (1 + eta)), 2);
+    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, exponent);
+}
+
+color CookTorranceDielectric::f_r(const ray& r_in, const HitInfo& rec, const ray& scattered, float R) const {
+    vec3 V = -r_in.direction().unit_vector();
+    vec3 L = scattered.direction().unit_vector();
+    vec3 H = (V + L).unit_vector();
+    vec3 N = rec.normal;
+
+    float NoH = clamp(dot(N, H), 0.0, 1.0);
+    float NoV = clamp(dot(N, V), 0.0, 1.0);
+    float NoL = clamp(dot(N, L), 0.0, 1.0);
+
+    float D = MDF->D(NoH);
+    float G = MDF->G(NoV, NoL);
+    
+    color R_col = color(R, R, R) * albedo;
+
+    color num = D * G * R_col;
+    float denom = (4.0 * fmax(NoV, 0.001) * fmax(NoL, 0.001));
+
+    return num / denom;
+}
+
+float CookTorranceDielectric::pdf_r(const ray& r_in, const HitInfo& rec, const ray& scattered, float R) const {
+    vec3 V = -r_in.direction().unit_vector();
+    vec3 L = scattered.direction().unit_vector();
+    vec3 H = (V + L).unit_vector();
+    vec3 N = rec.normal;
+
+    float NoH = clamp(dot(N, H), 0.0, 1.0);
+    float VoH = clamp(dot(V, H), 0.0, 1.0);
+
+    float D = MDF->D(NoH);
+    // Convert D(N·H) to pdf based on the microfacet normal distribution.
+    // The Jacobian of the half-vector reflection transformation is |4 * (V·H)|.
+    // This accounts for the change in area density when mapping from H to L.
+    float jacobian = 4.0 * abs(dot(V, H));
+    if (jacobian < 0.0001) return 0;
+
+    return (D * R) / jacobian;
+}
+
+float CookTorranceDielectric::pdf_t(const ray& r_in, const HitInfo& rec, const ray& scattered, float T) const {
+    double etap = rec.front_face ? (1.0/eta) : eta;
+
+    vec3 wo = -r_in.direction().unit_vector();
+    vec3 wi = scattered.direction().unit_vector();
+    vec3 wn = rec.normal;
+    vec3 wm = rec.microfacet_normal;
+    vec3 h = (wo + wi).unit_vector();
+
+    float denom = (dot(wi, wm) + dot(wo, wm) / etap) * (dot(wi, wm) + dot(wo, wm) / etap);
+    float dwm_dwi = fabs(dot(wi, wm)) / denom;
+    float NoM = dot(wm, wn);
+    float D = MDF->D(NoM);
+    return D * dwm_dwi * T;
+}
+
+color CookTorranceDielectric::f_t(const ray& r_in, const HitInfo& rec, const ray& scattered, float T) const {
+    double etap = rec.front_face ? (1.0/eta) : eta;
+
+    vec3 wo = -r_in.direction().unit_vector();
+    vec3 wi = scattered.direction().unit_vector();
+    vec3 wn = rec.normal;
+    vec3 wm = rec.microfacet_normal;
+    vec3 h = (wo + wi).unit_vector();
+
+    float NoM = dot(wm, wn);
+    float NoO = dot(wn, wo);
+    float NoI = dot(wn, wi);
+    float D = MDF->D(NoM);
+    float G = MDF->G(fabs(NoO), fabs(NoI));
+    color T_col = color(T, T, T) * albedo;
+    color num = D * G * T_col;
+
+    float IoM = dot(wi, wm);
+    float OoM = dot(wo, wm);
+    float denom = (IoM + OoM / etap) * (IoM + OoM / etap);
+    float dotabs = fabs(IoM * OoM / (dot(wi, wn) * dot(wo, wn) * denom)); // 1: e+14, 2: inf
+    return num * dotabs;
+} 
+// ===================== ISOTROPIC ===============================
+bool isotropic::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+    scattered = ray(rec.pos, random_unit_vector(), r_in.time());
+    return true;
+}
+color isotropic::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    return albedo / (4 * pi);
+}
+double isotropic::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    return 1 / (4 * pi);
+};
+// ===================== PIXEL LAMBERTIAN ===============================
+bool pixel_lambertian::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const{
+    float t = random_double();
+    color4 val = albedo->value(rec.u, rec.v);
+    if (t > val.A) { // transparent
+        rec.transparent = true;
+        scattered = ray(rec.pos, r_in.direction(), 0.0);
+    } else {
+        onb uvw;
+        uvw.build_from_w(rec.normal);
+        auto scatter_direction = uvw.local(random_cosine_direction());
+        scattered = ray(rec.pos, scatter_direction, r_in.time());
+    }
+    
+    
+    return true;
+}
+
+color pixel_lambertian::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    if (!rec.transparent) {
+        return albedo->value(rec.u, rec.v).RGB / pi;
+    } else {
+        return color(1.0, 1.0, 1.0) / pi;
+    }
+}
+
+double pixel_lambertian::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    auto cos_theta = dot(rec.normal, scattered.direction().unit_vector());
+    if (!rec.transparent) { return fmax(0.0, cos_theta / pi); }
+    else { return fabs(cos_theta) / pi; }
+}
+
+BSDFSample pixel_lambertian::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    BSDFSample sample_data;
+    // Sample the microfacet distribution to get the scatter direction.
+    color attenuation; // placeholder until it gets removed from the scatter function header
+    sample_data.scatter = scatter(r_in, rec, attenuation, scattered);
+    sample_data.scatter_direction = scattered.direction().unit_vector();
+
+    // Sample the BRDF for the value
+    sample_data.bsdf_value = generate(r_in, scattered, rec);
+
+    // Find the PDF for the MDF
+    sample_data.pdf_value = pdf(r_in, scattered, rec);
+    if (rec.transparent) { sample_data.type = BSDF_TYPE::TRANSPARENT; }
+    return sample_data;
+}
+
+// ===================== MIXTURE ===============================
+bool MixtureBSDF::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
+    rec.rand = random_double();
+    int mat_ix = chooseSampleMaterial(rec.rand);
+    if (mat_ix >= 0) { // weights is valid
+        return mats[mat_ix]->scatter(r_in, rec, attenuation, scattered);
+    } else {
+        return false;
+    }
+}
+
+// assumes that scatter has already been called or sample has already been called, and thus rand is already generated.
+color MixtureBSDF::generate(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    int mat_ix = chooseSampleMaterial(rec.rand);
+    if (mat_ix >= 0) { return mats[mat_ix]->generate(r_in, scattered, rec);
+    } else { return color(1,1,1); }
+}
+
+double MixtureBSDF::pdf(const ray& r_in, const ray& scattered, const HitInfo& rec) const {
+    int mat_ix = chooseSampleMaterial(rec.rand);
+    if (mat_ix >= 0) { return mats[mat_ix]->pdf(r_in, scattered, rec);
+    } else { return 1.0; }
+}
+
+// NOTE: ASSUMES WEIGHTS IS AT LEAST OF SIZE 1 OTHERWISE BEHAVIOUR IS UNDEFINED
+BSDFSample MixtureBSDF::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    rec.rand = random_double();
+    int mat_ix = chooseSampleMaterial(rec.rand);
+    if (mat_ix >= 0) { return mats[mat_ix]->sample(r_in, rec, scattered);
+    } else {
+        BSDFSample sample_data;
+        return sample_data;
+    }
+}
+int MixtureBSDF::chooseSampleMaterial(float rand) const {
+    float cumulative_weight = 0.0f;
+    for (size_t i = 0; i < weights.size(); i++) {
+        cumulative_weight += weights[i];
+        if (rand < cumulative_weight) {
+            return i;
+        }
+    }
+    return (int)weights.size() - 1;
+}
+// ===================== LAYERED ===============================
+BSDFSample LayeredBSDF::sample(const ray& r_in, HitInfo& rec, ray& scattered) const {
+    // Return in case calculating a full simulation becomes impossible or irrelevant
+    BSDFSample absorbed; absorbed.scatter = false;
+    
+    ray r = r_in;
+    HitInfo rec_manip = rec;
+    vec3 outward_normal = rec.front_face ? rec.normal : -rec.normal;
+
+    bool on_top = rec_manip.front_face;
+
+    std::shared_ptr<material> current = on_top ? top : bottom;
+    BSDFSample bs = current->sample(r, rec_manip, scattered);
+
+    color f = bs.bsdf_value;
+    float pdf = bs.pdf_value;
+
+    float distance_in_from_top = rec.front_face ? 0.0 : 1.0; // amount of distance ray has traveled into between the layers
+    bool prev_medium = false; // previous hit was a medium
+
+    if (!bs.scatter) { return absorbed; }
+    if (bs.type != BSDF_TYPE::TRANSMISSION && bs.type != BSDF_TYPE::TRANSPARENT) {
+        bs.type = (dot(bs.scatter_direction, rec.normal) < 0) ? BSDF_TYPE::TRANSMISSION : BSDF_TYPE::SPECULAR;
+        return bs;
+    }
+    for (int depth = 0; depth < termination; depth++) {
+        
+        // Follow random walk through layers to sample layered BSDF
+        on_top = !on_top;
+        current = on_top ? top : bottom;
+
+        // Possibly terminate layered BSDF sampling with Russian Roulette
+        float rrBeta = fmax(fmax(f.x(), f.y()), f.z()) / bs.pdf_value;
+        if (depth > 3 && rrBeta < 0.25) { // rrBeta < 0.50 reduces by more, but probably reduces accuracy
+            float q = fmax(0, 1-rrBeta);
+            if (random_double() < q) {
+                if (on_top == rec.front_face) {
+                    bs.type = (dot(bs.scatter_direction, rec.normal) < 0) ? BSDF_TYPE::TRANSMISSION : BSDF_TYPE::SPECULAR;
+                    return bs;
+                }
+            }
+            pdf *= 1 - q;
+        }
+
+        if (!prev_medium) { f = f * fabs(dot(rec_manip.normal, bs.scatter_direction.unit_vector())); }
+        prev_medium = false;
+        r = ray(rec_manip.pos - bs.scatter_direction, bs.scatter_direction, r.time());
+
+        if (medium) {
+            float distance = medium->particleDistance();
+            float distance_to_barrier;
+            if (on_top) { distance_to_barrier = distance_in_from_top; }
+            else { distance_to_barrier = thickness - distance_in_from_top; }
+            if (distance < distance_to_barrier) { // collide with particle
+                distance_in_from_top += (on_top * -distance) + (!on_top * distance);
+                bs = medium->phase->sample(r, rec_manip, scattered);
+                prev_medium = true;
+
+                f = f * bs.bsdf_value;
+                pdf = pdf * bs.pdf_value;
+                bs.bsdf_value = f;
+                bs.pdf_value = pdf;
+
+                if (dot(bs.scatter_direction, outward_normal) > 0) {
+                    on_top = false; // because it will flip at the start of the next loop
+                    rec_manip.front_face = false;
+                    rec_manip.normal = -outward_normal;
+                } else {
+                    on_top = true; // because it will flip at the start of the next loop
+                    rec_manip.front_face = true;
+                    rec_manip.normal = outward_normal;
+                }
+
+                continue;
+            } else {
+                if (on_top) { distance_in_from_top = 0.0;} else { distance_in_from_top = 1.0; }
+            }
+        }
+
+        bs = current->sample(r, rec_manip, scattered);
+        f = f * bs.bsdf_value;
+        pdf = pdf * bs.pdf_value;
+        bs.bsdf_value = f;
+        bs.pdf_value = pdf;
+
+        if (bs.type == BSDF_TYPE::TRANSMISSION) {
+            bs.type = (dot(bs.scatter_direction, rec.normal) < 0) ? BSDF_TYPE::TRANSMISSION : BSDF_TYPE::SPECULAR;
+            return bs;
+        }
+
+        // Flip since coming from the bottom!
+        rec_manip.front_face = !rec_manip.front_face;
+        rec_manip.normal = -rec_manip.normal;
+    }
+    bs.type = (dot(bs.scatter_direction, rec.normal) < 0) ? BSDF_TYPE::TRANSMISSION : BSDF_TYPE::SPECULAR;
+    return bs;
+}
\ No newline at end of file
diff --git a/src/materials/medium.cc b/src/materials/medium.cc
new file mode 100644
index 0000000..9970d1b
--- /dev/null
+++ b/src/materials/medium.cc
@@ -0,0 +1,77 @@
+#include "medium.h"
+
+Medium::Medium(float density, shared_ptr<material> phase_function) 
+    : density{density}, phase{phase_function} {}
+
+float Medium::particleDistance() const {
+    float neg_inv_density = -1 / density;
+    auto hit_distance = neg_inv_density * log(random_double());
+    return hit_distance;
+}
+
+float Medium::transmittance(const vec3& x, const vec3& y) const {
+    float exponent = -(density * fabs((x - y).length()));
+    return pow(euler, exponent);
+}
+
+float Medium::transmittance(float dist) const {
+    float exponent = -(density * dist);
+    return pow(euler, exponent);
+}
+
+
+// MediumStack
+
+MediumRecord::MediumRecord(point3 initial_position) : recent_position{initial_position} {}
+
+bool MediumRecord::contains(std::shared_ptr<Medium> vol) const {
+    return std::find(mediums.begin(), mediums.end(), vol) != mediums.end();
+}
+
+void MediumRecord::hitVolume(shared_ptr<Medium> vol, point3 hitPoint) {
+    float distance_traveled = (recent_position - hitPoint).length();
+    if (!mediums.empty()) {
+        transmittance *= highest_density_volume->transmittance(distance_traveled);
+        if (!vol) { return; }
+        if (!contains(vol)) { // entering
+            mediums.push_back(vol);
+            if (vol->density > highest_density_volume->density) { highest_density_volume = vol; }
+        } else { // exiting
+            mediums.erase(std::remove(mediums.begin(), mediums.end(), vol), mediums.end());
+
+            float highestDensity = -1.0f;
+            if (vol == highest_density_volume) { // reset highest density volume by finding it
+                for (const auto& ptr : mediums) {
+                    if (ptr->density > highestDensity) {
+                        highestDensity = ptr->density;
+                        highest_density_volume = ptr;
+                    }
+                }
+            }
+        }
+    } else { // mediums empty
+        if (!vol) { return; }
+        highest_density_volume = vol;
+        mediums.push_back(vol);
+        recent_position = hitPoint;
+    }
+}
+
+shared_ptr<Medium> MediumRecord::particleDistance(float& dist) const {
+    float min_dist;
+    shared_ptr<Medium> p = nullptr;
+    for (size_t i = 0; i < mediums.size(); ++i) {
+        float ndist = mediums[i]->particleDistance();
+        if (i == 0) {
+            min_dist = ndist;
+            p = mediums[i];
+        } else {
+            if (ndist < min_dist) {
+                min_dist = ndist;
+                p = mediums[i];
+            }
+        }
+    }
+    dist = min_dist;
+    return p;
+}
\ No newline at end of file
diff --git a/src/objects/geometry.cc b/src/objects/geometry.cc
index 5db1d15..0429c03 100644
--- a/src/objects/geometry.cc
+++ b/src/objects/geometry.cc
@@ -1,3 +1,8 @@
 #include "geometry.h"
 
-Geometry::Geometry(vec3 position, RTCGeometry geom) : geom{geom}, Visual(position) {}
\ No newline at end of file
+Geometry::Geometry(vec3 position, RTCGeometry geom) : geom{geom}, Visual(position) {}
+
+Geometry::Geometry(shared_ptr<Geometry> geom_ptr) : Visual(geom_ptr->position) {
+    geom = geom_ptr->geom;
+    rtcRetainGeometry(geom);
+}
\ No newline at end of file
diff --git a/src/objects/light.cc b/src/objects/light.cc
index 84a0dda..35d2cd1 100644
--- a/src/objects/light.cc
+++ b/src/objects/light.cc
@@ -2,7 +2,7 @@
 
 emissive::emissive(color emission_color) : emission_color{emission_color} {}
 
-bool emissive::scatter(const ray& r_in, const HitInfo& rec, color& attenuation, ray& scattered) const {
+bool emissive::scatter(const ray& r_in, HitInfo& rec, color& attenuation, ray& scattered) const {
     return false;
 }
 
diff --git a/src/objects/primitives/quad_primitive.cc b/src/objects/primitives/quad_primitive.cc
index 1b5c456..2bc87df 100644
--- a/src/objects/primitives/quad_primitive.cc
+++ b/src/objects/primitives/quad_primitive.cc
@@ -41,10 +41,15 @@ HitInfo QuadPrimitive::getHitInfo(const ray& r, const vec3& p, const float t, un
     return record;
 }
 
-vec3 QuadPrimitive::getV() {
-    return v;
+point3 QuadPrimitive::sample(const HitInfo& rec) const {
+    vec3 random_point = this->position + random_double() * this->u + random_double() * this->v;
+    return random_point;
 }
-
-vec3 QuadPrimitive::getU() {
-    return u;
-}
\ No newline at end of file
+double QuadPrimitive::pdf(const HitInfo& light_record, ray sample_ray) const {
+    auto distance_squared = light_record.t * light_record.t * sample_ray.direction().length_squared();
+    auto cosine = fabs(dot(sample_ray.direction().unit_vector(), light_record.normal));
+    double area = cross(u,v).length();
+    return distance_squared / (cosine * area);
+}
+vec3 QuadPrimitive::getV() { return v; }
+vec3 QuadPrimitive::getU() { return u; }
diff --git a/src/objects/primitives/sphere_primitive.cc b/src/objects/primitives/sphere_primitive.cc
index c6090c3..4a18b7f 100644
--- a/src/objects/primitives/sphere_primitive.cc
+++ b/src/objects/primitives/sphere_primitive.cc
@@ -33,6 +33,23 @@ HitInfo SpherePrimitive::getHitInfo(const ray& r, const vec3& p, const float t,
     return record;
 }
 
+point3 SpherePrimitive::sample(const HitInfo& rec) const {
+    vec3 direction = position - rec.pos;
+    vec3 to_sphere = random_in_hemisphere(direction.unit_vector());
+    return position + (to_sphere * radius);
+    // auto distance_squared = direction.length_squared();
+    // onb uvw;
+    // uvw.build_from_w(direction);
+    // return uvw.local(random_to_sphere(radius, distance_squared));
+}
+
+double SpherePrimitive::pdf(const HitInfo& rec, ray sample_ray) const {
+    auto cos_theta_max = sqrt(1 - radius*radius/(position - sample_ray.origin()).length_squared());
+    auto solid_angle = 2*pi*(1-cos_theta_max);
+    return 1 / solid_angle;
+    //return (4 * pi);
+}
+
 void SpherePrimitive::get_sphere_uv(const point3& p, double& u, double& v) {
     // p: a given point on the sphere of radius one, centered at the origin.
     // u: returned value [0,1] of angle around the Y axis from X=-1.
@@ -46,4 +63,16 @@ void SpherePrimitive::get_sphere_uv(const point3& p, double& u, double& v) {
 
     u = phi / (2*pi);
     v = theta / pi;
+}
+
+vec3 SpherePrimitive::random_to_sphere(double radius, double distance_squared) {
+    auto r1 = random_double();
+    auto r2 = random_double();
+    auto z = 1 + r2*(sqrt(1-radius*radius/distance_squared) - 1);
+
+    auto phi = 2*pi*r1;
+    auto x = cos(phi)*sqrt(1-z*z);
+    auto y = sin(phi)*sqrt(1-z*z);
+
+    return vec3(x, y, z);
 }
\ No newline at end of file
diff --git a/src/objects/volume.cc b/src/objects/volume.cc
new file mode 100644
index 0000000..98e9e10
--- /dev/null
+++ b/src/objects/volume.cc
@@ -0,0 +1,18 @@
+#include "volume.h"
+
+Volume::Volume(shared_ptr<Medium> medium, shared_ptr<Geometry> geom_ptr, RTCDevice device) : medium{medium}, geom_ptr{geom_ptr}, Geometry(geom_ptr) {
+    volume_scene = rtcNewScene(device);
+    unsigned int geomID = rtcAttachGeometry(volume_scene, geom_ptr->geom);
+    rtcReleaseGeometry(geom_ptr->geom);
+    rtcCommitScene(volume_scene);
+}
+
+shared_ptr<material> Volume::materialById(unsigned int geomID) const {
+    return medium->phase;
+}
+
+HitInfo Volume::getHitInfo(const ray& r, const vec3& p, const float t, unsigned int geomID) const {
+    HitInfo rec = geom_ptr->getHitInfo(r, p, t, geomID);
+    rec.medium = true;
+    return rec;
+}
\ No newline at end of file
diff --git a/src/output/output.cc b/src/output/output.cc
index fc3af68..b166b55 100644
--- a/src/output/output.cc
+++ b/src/output/output.cc
@@ -91,7 +91,8 @@ void output(RenderData& render_data, Camera& cam, std::shared_ptr<Scene> scene_p
         }
     }
 
-    if (config.verbose) {
+    if (true) {
+    // if (config.verbose) {
         auto current_time = std::chrono::high_resolution_clock::now();
         auto elapsed_time = std::chrono::duration_cast<std::chrono::milliseconds>(current_time - start_time).count();
         double time_seconds = elapsed_time / 1000.0;
diff --git a/src/parsers/cli_parser.cc b/src/parsers/cli_parser.cc
index 42c39aa..f4ff089 100644
--- a/src/parsers/cli_parser.cc
+++ b/src/parsers/cli_parser.cc
@@ -1,3 +1,9 @@
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <cstdlib>
+#include "render.h"
+
 #include "cli_parser.hh"
 
 void outputHelpGuide(std::ostream& out) {
@@ -14,7 +20,7 @@ void outputHelpGuide(std::ostream& out) {
         << " -h,  --help                           Show this help message.\n"
         << " -V,  --verbose                        Enables more descriptive messages of scenes and rendering process.\n"
         << " -T,  --threads <amt>                  If multithreading is enabled, sets amount of threads used.\n"
-        << " -Vx, --vectorization <batch_size>     Set SIMD vectorization batch size [0|4|8|16]. If NONE = 0, do not enable the flag.\n";
+        << " -Vx, --vectorization <batch_size>     [UNSUPPORTED/DEPRECATED] Set SIMD vectorization batch size [0|4|8|16]. If NONE = 0, do not enable the flag.\n";
     exit(0);
 }
 
@@ -96,14 +102,16 @@ Config parseArguments(int argc, char* argv[]) {
         } 
 
         else if(arg == "-Vx" || arg == "--vectorization") {
-            int choice = checkValidIntegerInput(i, argc, argv, "-Vx/--vectorization");
-            if (choice == 0 || choice == 4 || choice == 8 || choice == 16) config.vectorization = choice;
-            else throw std::invalid_argument("Error: Invalid option for --vectorization [1|4|8|16]. Use '--help' for more information.");
+            throw std::runtime_error("Error (DEPRECATION): -Vx/--vectorization not supported in v0.1.5.");
+            // CURRENTLY DEPRECATED AND NEEDS UPDATE TO 0.1.5
+            // int choice = checkValidIntegerInput(i, argc, argv, "-Vx/--vectorization");
+            // if (choice == 0 || choice == 4 || choice == 8 || choice == 16) config.vectorization = choice;
+            // else throw std::invalid_argument("Error: Invalid option for --vectorization [1|4|8|16]. Use '--help' for more information.");
         } 
 
         else if(arg == "-v" || arg == "--version") {
             config.showVersion = true;
-            std::cout << "caitlyn version 0.1.3" << std::endl;
+            std::cout << "caitlyn version 0.1.5" << std::endl;
         } 
         
         else if(arg == "-h" || arg == "--help") {
diff --git a/src/parsers/csr_parser.cc b/src/parsers/csr_parser.cc
index b6e5514..604c277 100644
--- a/src/parsers/csr_parser.cc
+++ b/src/parsers/csr_parser.cc
@@ -7,8 +7,13 @@ std::shared_ptr<Scene> CSRParser::parseCSR(std::string& filePath, RTCDevice devi
     file = std::ifstream(filePath);
     std::string line;
     std::map<std::string, std::shared_ptr<material>> materials;
+    std::map<std::string, std::shared_ptr<emissive>> emissives;
     std::map<std::string, std::shared_ptr<texture>> textures;
     std::map<std::string, std::shared_ptr<Primitive>> primitives;
+    std::map<std::string, std::shared_ptr<Primitive>> medium_primitives;
+
+    std::map<std::string, std::shared_ptr<Medium>> mediums;
+    std::map<std::string, std::shared_ptr<Volume>> volumes;
     
     if (!file.is_open() || !file.good()) {
         rtcReleaseDevice(device);
@@ -27,19 +32,23 @@ std::shared_ptr<Scene> CSRParser::parseCSR(std::string& filePath, RTCDevice devi
 
     while (getNextLine(file, line)) {
         line = trim(line);
-        if (startsWith(line, "Material")) {
+        if (startsWith(line, "Sky")) {
+            std::string top, bottom;
+            getNextLine(file, top); getNextLine(file, bottom);
+            scene_ptr->set_sky_colour(readXYZProperty(bottom), readXYZProperty(top));
+        } else if (startsWith(line, "Material")) {
             // Extract material ID from brackets (e.g., Material[Lambertian] -> Lambertian)
             auto idStart = line.find('[') + 1;
             auto idEnd = line.find(']');
             std::string materialType = line.substr(idStart, idEnd - idStart);
-            if (materialType == "Lambertian") {
+            if (materialType == "Diffuse") {
                 std::string materialId, texture;
                 getNextLine(file, materialId); getNextLine(file, texture);
                 std::string texture_id = readStringProperty(texture);
                 if (texture_id == "no") {
-                    std::string albedo;
-                    getNextLine(file, albedo);
-                    materials[readStringProperty(materialId)] = std::make_shared<lambertian>(readXYZProperty(albedo));  
+                    std::string albedo, roughness;
+                    getNextLine(file, albedo); getNextLine(file, roughness);
+                    materials[readStringProperty(materialId)] = std::make_shared<OrenNayar>(readXYZProperty(albedo), readDoubleProperty(roughness));  
                 } else {
                     std::shared_ptr<PixelImageTexture> plamb = std::dynamic_pointer_cast<PixelImageTexture>(textures[texture_id]);
                     if (plamb) { // texture is a pixel lambert
@@ -49,17 +58,41 @@ std::shared_ptr<Scene> CSRParser::parseCSR(std::string& filePath, RTCDevice devi
                     }
                 }
             } else if (materialType == "Metal") {
-                std::string materialId, albedo, fuzz;
-                getNextLine(file, materialId); getNextLine(file, albedo); getNextLine(file, fuzz);
-                materials[readStringProperty(materialId)] = std::make_shared<metal>(readXYZProperty(albedo), readDoubleProperty(fuzz));
+                std::string materialId, albedo, roughness;
+                getNextLine(file, materialId); getNextLine(file, albedo); getNextLine(file, roughness);
+                materials[readStringProperty(materialId)] = std::make_shared<CookTorrance>(readXYZProperty(albedo), readDoubleProperty(roughness));
             } else if (materialType == "Dielectric") {
-                std::string materialId, ir;
-                getNextLine(file, materialId); getNextLine(file, ir);
-                materials[readStringProperty(materialId)] = std::make_shared<dielectric>(readDoubleProperty(ir));
+                std::string materialId, albedo, eta, roughness, sheen;
+                getNextLine(file, materialId); getNextLine(file, albedo); getNextLine(file, eta); getNextLine(file, roughness); getNextLine(file, sheen);
+                materials[readStringProperty(materialId)] = std::make_shared<CookTorranceDielectric>(
+                    readXYZProperty(albedo), readDoubleProperty(eta), readDoubleProperty(roughness), (int)readDoubleProperty(sheen)
+                );
             } else if (materialType == "Emissive") {
                 std::string materialId, rgb, strength;
                 getNextLine(file, materialId); getNextLine(file, rgb); getNextLine(file, strength);
                 materials[readStringProperty(materialId)] = std::make_shared<emissive>( (readDoubleProperty(strength) * readXYZProperty(rgb)) );
+                // emissives[readStringProperty(materialId)] = std::make_shared<emissive>( (readDoubleProperty(strength) * readXYZProperty(rgb)) );
+            } else if (materialType == "Mixture") {
+                std::vector<float> weights;
+                std::vector<std::shared_ptr<material>> mats;
+                std::string materialId, number;
+                getNextLine(file, materialId); getNextLine(file, number);
+                int num_mats = (int)readDoubleProperty(number);
+                for (int i=0; i<num_mats; i++) {
+                    std::string mixed, weight;
+                    getNextLine(file, mixed); getNextLine(file, weight);
+                    weights.push_back(readDoubleProperty(weight));
+                    mats.push_back(materials[readStringProperty(mixed)]);
+                }
+                materials[readStringProperty(materialId)] = make_shared<MixtureBSDF>(weights, mats);
+            } else if (materialType == "Layered") {
+                std::string materialId, top, bottom, medium;
+                getNextLine(file, materialId); getNextLine(file, top); getNextLine(file, bottom); getNextLine(file, medium);
+                materials[readStringProperty(materialId)] = make_shared<LayeredBSDF>(
+                    materials[readStringProperty(top)],
+                    materials[readStringProperty(bottom)],
+                    (readStringProperty(medium) == "no") ? nullptr : mediums[readStringProperty(medium)]
+                );
             } else {
                 rtcReleaseDevice(device);
                 throw std::runtime_error("Material type UNDEFINED: Material[Lambertian|Metal|Dielectric|Emissive]");
@@ -86,23 +119,47 @@ std::shared_ptr<Scene> CSRParser::parseCSR(std::string& filePath, RTCDevice devi
                 throw std::runtime_error("Texture type UNDEFINED: Texture[Checker|Image|Noise]");
             }
         } else if (startsWith(line, "Sphere")) {
-            std::string id, position, material, radius;
+            std::string id, position, material, radius, medium;
             getNextLine(file, id); getNextLine(file, position); getNextLine(file, material); getNextLine(file, radius);
-            auto sphere = make_shared<SpherePrimitive>(readXYZProperty(position), materials[readStringProperty(material)], readDoubleProperty(radius), device);
-            primitives[readStringProperty(id)] = sphere;
-            scene_ptr->add_primitive(sphere);
+            getNextLine(file, medium);
+            // bool usesEmissive = (emissives.find(readStringProperty(material)) != emissives.end());
+            auto sphere = make_shared<SpherePrimitive>(
+                readXYZProperty(position), 
+                // usesEmissive ? emissives[readStringProperty(material)] : 
+                materials[readStringProperty(material)], 
+                readDoubleProperty(radius), device
+            );
+            if (!readBooleanProperty(medium)) { 
+                primitives[readStringProperty(id)] = sphere;
+                scene_ptr->add_primitive(sphere);
+            } else {
+                medium_primitives[readStringProperty(id)] = sphere;
+            }
+            // if (usesEmissive) { scene_ptr->add_physical_light(sphere); }
         } else if (startsWith(line, "Quad")) {
             std::string id, position, u, v, material;
             getNextLine(file, id); getNextLine(file, position); getNextLine(file, u); getNextLine(file, v); getNextLine(file, material);
-            auto quad = make_shared<QuadPrimitive>(readXYZProperty(position), readXYZProperty(u), readXYZProperty(v), materials[readStringProperty(material)], device);
+            // bool usesEmissive = (emissives.find(readStringProperty(material)) != emissives.end());
+            auto quad = make_shared<QuadPrimitive>(
+                readXYZProperty(position), readXYZProperty(u), readXYZProperty(v), 
+                // usesEmissive ? emissives[readStringProperty(material)] : 
+                materials[readStringProperty(material)], 
+                device
+            );
             primitives[readStringProperty(id)] = quad;
             scene_ptr->add_primitive(quad);
+            // if (usesEmissive) { scene_ptr->add_physical_light(quad); }
         } else if (startsWith(line, "Box")) {
-            std::string id, position, a, b, c, material;
+            std::string id, position, a, b, c, material, medium;
             getNextLine(file, id); getNextLine(file, position); getNextLine(file, a); getNextLine(file, b); getNextLine(file, c); getNextLine(file, material);
+            getNextLine(file, medium);
             auto box = make_shared<BoxPrimitive>(readXYZProperty(position), readXYZProperty(a), readXYZProperty(b), readXYZProperty(c), materials[readStringProperty(material)], device);
-            primitives[readStringProperty(id)] = box;
-            scene_ptr->add_primitive(box);
+            if (!readBooleanProperty(medium)) { 
+                primitives[readStringProperty(id)] = box;
+                scene_ptr->add_primitive(box);
+            } else {
+                medium_primitives[readStringProperty(id)] = box;
+            }
         } else if (startsWith(line, "Instance")) {
             auto idStart = line.find('[') + 1;
             auto idEnd = line.find(']');
@@ -140,25 +197,40 @@ std::shared_ptr<Scene> CSRParser::parseCSR(std::string& filePath, RTCDevice devi
                 auto instance = make_shared<QuadPrimitiveInstance>(instance_ptr, transform, device);
                 scene_ptr->add_primitive_instance(instance, device);
             } else if (instanceType == "BoxPrimitive") {
-                    std::string prim_id, translate;
-                    getNextLine(file, prim_id); getNextLine(file, translate);
-                    vec3 translateVector = readXYZProperty(translate);
-                    float transform[12] = {
-                        1, 0, 0, translateVector.x(),
-                        0, 1, 0, translateVector.y(),
-                        0, 0, 1, translateVector.z()
-                    };
-                    std::shared_ptr<BoxPrimitive> instance_ptr = std::dynamic_pointer_cast<BoxPrimitive>(primitives[readStringProperty(prim_id)]);
-                    if (!instance_ptr) {
-                        rtcReleaseDevice(device);
-                        throw std::runtime_error("Instance key ERROR: " + readStringProperty(prim_id) + " is not a BoxPrimitive!");
-                    }
-                    auto instance = make_shared<BoxPrimitiveInstance>(instance_ptr, transform, device);
-                    scene_ptr->add_primitive_instance(instance, device);
+                std::string prim_id, translate;
+                getNextLine(file, prim_id); getNextLine(file, translate);
+                vec3 translateVector = readXYZProperty(translate);
+                float transform[12] = {
+                    1, 0, 0, translateVector.x(),
+                    0, 1, 0, translateVector.y(),
+                    0, 0, 1, translateVector.z()
+                };
+                std::shared_ptr<BoxPrimitive> instance_ptr = std::dynamic_pointer_cast<BoxPrimitive>(primitives[readStringProperty(prim_id)]);
+                if (!instance_ptr) {
+                    rtcReleaseDevice(device);
+                    throw std::runtime_error("Instance key ERROR: " + readStringProperty(prim_id) + " is not a BoxPrimitive!");
+                }
+                auto instance = make_shared<BoxPrimitiveInstance>(instance_ptr, transform, device);
+                scene_ptr->add_primitive_instance(instance, device);
             } else {
                 rtcReleaseDevice(device);
                 throw std::runtime_error("Instance type UNDEFINED: Instance[SpherePrimitive|QuadPrimitive|BoxPrimitive]");
             }
+        } else if (startsWith(line, "Medium")) {
+            std::string medium_id, density, albedo;
+            getNextLine(file, medium_id); getNextLine(file, density); getNextLine(file, albedo);
+            auto medium = make_shared<Medium>(readDoubleProperty(density), make_shared<isotropic>(readXYZProperty(albedo)));
+            mediums[readStringProperty(medium_id)] = medium;
+        } else if (startsWith(line, "Volume")) {
+            std::string volume_id, medium_id, prim_id;
+            getNextLine(file, volume_id); getNextLine(file, medium_id); getNextLine(file, prim_id);
+            auto volume = make_shared<Volume>(
+                mediums[readStringProperty(medium_id)],
+                medium_primitives[readStringProperty(prim_id)],
+                device
+            );
+            volumes[readStringProperty(volume_id)] = volume;
+            scene_ptr->add_volume(volume);
         }
     }
 
diff --git a/src/render.cc b/src/render.cc
index 6d705b3..21521ba 100644
--- a/src/render.cc
+++ b/src/render.cc
@@ -1,5 +1,226 @@
 #include "render.h"
 
+color INVALID_SAMPLE = color(INT16_MIN, INT16_MIN, INT16_MIN);
+
+color trace_ray(const ray& r, std::shared_ptr<Scene> scene, int depth) {
+    HitInfo record;
+
+    color weight = color(1.0, 1.0, 1.0);
+    color accumulated_color = color(0,0,0);
+
+    ray r_in = r;
+    BSDF_TYPE incoming_type = BSDF_TYPE::DIFFUSE;
+
+    MediumRecord med_rec(r.origin());
+    // Trace rays in all volume scenes to check which ones we reside in
+    struct RTCRayHit vol_rayhit;
+    HitInfo vol_record;
+    for (const auto& ptr : scene->volumes) {
+        setupRayHit1(vol_rayhit, r_in);
+        rtcIntersect1(ptr->volume_scene, &vol_rayhit);
+        int targetID;
+        if (vol_rayhit.hit.instID[0] != RTC_INVALID_GEOMETRY_ID) {
+            targetID = vol_rayhit.hit.instID[0];
+        } else if (vol_rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) {
+            targetID = vol_rayhit.hit.geomID;
+        } else {
+            continue;
+        }
+        vol_record = ptr->getHitInfo(r_in, r_in.at(vol_rayhit.ray.tfar), vol_rayhit.ray.tfar, targetID);
+        if (!vol_record.front_face) { // inside the volume!
+            record.medium = false;
+            record = vol_record;
+            if (record.medium) {
+                med_rec.hitVolume(ptr->medium, r_in.origin());
+                incoming_type = BSDF_TYPE::SPECULAR;
+            }
+        }
+    }
+
+    for (int i=0; i<depth; i++) {
+        // Enable of disable direct light sampling (debug only, should always be enabled)
+        bool direct = false;
+        bool raymarched = false; // set to true if we are colliding with a medium particle and not a surface
+        std::shared_ptr<material> mat_ptr = nullptr;
+        ray scattered;
+        color attenuation;
+        struct RTCRayHit rayhit;
+        setupRayHit1(rayhit, r_in);
+
+        rtcIntersect1(scene->rtc_scene, &rayhit);
+
+        // Check for volume intersections
+        float particleDist;
+        shared_ptr<Medium> m_ptr = med_rec.particleDistance(particleDist);
+        if (m_ptr) { // we are in a medium
+            float istDist = rayhit.ray.tfar * r_in.direction().length();
+            if (particleDist < istDist) { // volume intersection
+                // Update record
+                record.pos = r_in.at(particleDist / r_in.direction().length());
+                raymarched = true;
+            }
+        }
+
+        if (!raymarched) { // process information of next geometry hit
+            int targetID;
+            if (rayhit.hit.instID[0] != RTC_INVALID_GEOMETRY_ID) {
+                targetID = rayhit.hit.instID[0];
+            } else if (rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) {
+                targetID = rayhit.hit.geomID;
+            } else {
+                vec3 unit_direction = r_in.direction().unit_vector();
+                auto t = 0.5*(unit_direction.y() + 1.0);
+
+                color sky = (1.0-t)*(scene->sky_top) + t*(scene->sky_bottom); // lerp formula (1.0-t)*start + t*endval
+                accumulated_color += weight * sky;
+                return accumulated_color;
+            }
+
+            std::shared_ptr<Geometry> geomhit = scene->geom_map[targetID];
+            mat_ptr = geomhit->materialById(targetID);
+            record.medium = false;
+            record = geomhit->getHitInfo(r_in, r_in.at(rayhit.ray.tfar), rayhit.ray.tfar, targetID);
+            if (record.medium) {
+                std::shared_ptr<Volume> volhit = std::dynamic_pointer_cast<Volume>(geomhit);
+                if (volhit) {
+                    med_rec.hitVolume(volhit->medium, r_in.at(rayhit.ray.tfar));
+                    r_in = ray(r_in.at(rayhit.ray.tfar), r_in.direction(), 0.0);
+                    incoming_type = BSDF_TYPE::SPECULAR;
+                    continue; // ignore edges of volumes? move to next bounce
+                }
+            }
+        } else {
+            mat_ptr = m_ptr->phase;
+        }
+        // Get emission contribution
+        color color_from_emission = mat_ptr->emitted(record.u, record.v, record.pos);
+
+        BSDFSample sample_data = mat_ptr->sample(r_in, record, scattered);
+        if (sample_data.type != BSDF_TYPE::DIFFUSE) { direct = false; }
+        if (incoming_type != BSDF_TYPE::DIFFUSE || sample_data.type == BSDF_TYPE::TRANSMISSION) {
+            accumulated_color += weight * color_from_emission;
+        } else if (direct == false) { accumulated_color += weight * color_from_emission; } else {
+            // To prevent double contribution of emission, only directly add if and only if:
+            // => we are directly hitting the light, i.e (i==0)
+            if (i == 0) { accumulated_color += weight * color_from_emission; }
+        }
+
+        // Direct Light Sampling
+        // Disabled as of v0.1.5, needs improvements!
+        // Additionally:
+        // => may be using wrong w in cos term
+        // => check if pdfs are 0 are not in place (invalid sample if so)
+        if (direct && color_from_emission.length() == 0.0) {
+            int N = (int)scene->physical_lights.size(); // amount of lights
+            for (auto& light_ptr : scene->physical_lights) { // only accounts for physical lights currently
+                // Create ray from hit point to the light
+                point3 sampled_point = light_ptr->sample(record);
+                vec3 light_dir = (sampled_point - record.pos).unit_vector(); // direction from hit point to the light
+                ray light_ray = ray(record.pos, light_dir, 0.0);
+
+                float distWithinMedium = (sampled_point - record.pos).length(); // distance to light
+
+                // MultiIntersect to light to capture medium boundaries
+                MediumRecord light_med_rec(record.pos);
+                light_med_rec.mediums = med_rec.mediums;
+                light_med_rec.highest_density_volume = med_rec.highest_density_volume;
+                std::vector<int> ids;
+                std::vector<float> tfars;
+                // errors warning: overflow in conversion from 'float' to 'int' changes value from '+Inff' to '2147483647' [-Woverflow]
+                // MultiIntersect(std::numeric_limits<float>::infinity(), light_ray, scene->rtc_scene, ids, tfars);
+                int intersections_to_accept = 10;
+                MultiIntersect(intersections_to_accept, light_ray, scene->rtc_scene, ids, tfars); // assume output ids.length() == tfars.length()
+                
+                float epsilon = distWithinMedium * 1e-5;
+                int count = 1;
+                while ((int)ids.size() == 0) {
+                    // This case may occur if the ray sampled runs parallel to a quad.
+                    // To remedy, we apply a small offset by the hit normal.
+                    point3 new_pos = record.pos + count*epsilon * record.normal;
+                    light_ray = ray(new_pos, (sampled_point - new_pos).unit_vector(), 0.0);
+                    MultiIntersect(intersections_to_accept, light_ray, scene->rtc_scene, ids, tfars);
+                    count++;
+                }
+                
+                // In this section, we fire trace each recorded intersection from pos -> light
+                // It is hoped/assumed that we find it within 'intersections_to_accept' intersections or we find some obscurement
+                bool non_medium_encountered = false;
+                std::shared_ptr<Geometry> light_geomhit;
+                int light_id;
+                int light_tfar;
+                for (size_t j = 0; j < ids.size(); j++) {
+                    int id = ids[j];
+                    float tfar = tfars[j];
+                    light_geomhit = scene->geom_map[id];
+                    std::shared_ptr<Volume> possible_volume_hit = std::dynamic_pointer_cast<Volume>(light_geomhit);
+                    if (!possible_volume_hit) { // non volume encountered
+                        if (light_geomhit == light_ptr) { // hit the light
+                            light_id = id;
+                            light_tfar = tfar;
+                            light_med_rec.hitVolume(nullptr, light_ray.at(tfar));
+                            break;
+                        }
+                        non_medium_encountered = true;
+                        break;
+                    } else { // one of the intersections was a volume
+                        light_med_rec.hitVolume(possible_volume_hit->medium, light_ray.at(tfar));
+                    }
+                }
+               
+                if (!non_medium_encountered) { // if it is the light, we are not obscured from the light
+                    // Store hit data of tracing the ray from here to the light
+                    HitInfo light_record;
+                    light_record = light_geomhit->getHitInfo(light_ray, light_ray.at(light_tfar), light_tfar, light_id);
+                    
+                    // Get the light's material
+                    std::shared_ptr<material> light_mat_ptr = light_geomhit->materialById(light_id);
+
+                    // Sample BSDF of hit point with incoming light
+                    ray light_scattered;
+                    
+                    BSDFSample light_sample_data;
+                    color att;
+                    // We do NOT call the above line because it would sample a possibly different microfacet normal
+                    // than what is already sampled previous to the Direct Light Sampling (for complex BSDFs that use microfacets)
+                    // Both generate and pdf assume that, if a microfacet normal is needed, it is already defined. Thus, we use the previous
+                    // and pass in the same HitInfo.
+                    light_sample_data.bsdf_value = mat_ptr->generate(r_in, light_ray, record);
+                    light_sample_data.pdf_value = mat_ptr->pdf(r_in, light_ray, record);
+                    
+                    // Find pdf for the light hit point
+                    double light_pdf_value = light_ptr->pdf(light_record, light_ray);
+
+                    // Find contribution of light using MIS power heuristic of light_pdf and sample pdf
+                    double light_cos_theta = fabs(dot(record.normal, light_dir));
+                    color light_contribution = weight * light_sample_data.bsdf_value * light_cos_theta
+                            * MIS::power_heuristic<MIS::EVAL_WEIGHT>(light_pdf_value, light_sample_data.pdf_value) / light_pdf_value;
+                    float transmittance_coeff = light_med_rec.transmittance;
+                    light_contribution *= transmittance_coeff;
+                    // Get emission of light
+                    color light_Le = light_mat_ptr->emitted(light_record.u, light_record.v, light_record.pos);
+                    accumulated_color += light_contribution * light_Le / N;
+                }
+            }
+        }
+
+        // Indirect ray contribution
+        if (!sample_data.scatter) {
+            return accumulated_color;
+        }
+        double cos_theta = fabs(dot(record.normal, sample_data.scatter_direction.unit_vector()));
+        if (!raymarched) {
+            if (sample_data.pdf_value == 0) { return INVALID_SAMPLE; }
+            weight = weight * (sample_data.bsdf_value * cos_theta / sample_data.pdf_value);
+        } else {
+            if (sample_data.pdf_value == 0) { return INVALID_SAMPLE; }
+            weight = weight * (sample_data.bsdf_value / sample_data.pdf_value);
+        }
+        r_in = scattered;
+        incoming_type = sample_data.type;
+	}
+    return accumulated_color;
+}
+
 void setRenderData(RenderData& render_data, const float aspect_ratio, const int image_width, const int samples_per_pixel, const int max_depth) {
     const int image_height = static_cast<int>(image_width / aspect_ratio);
     render_data.image_width = image_width;
@@ -62,17 +283,53 @@ void render_scanlines(int lines, int start_line, std::shared_ptr<Scene> scene_pt
     int samples_per_pixel   = data.samples_per_pixel;
     int max_depth           = data.max_depth;
 
+    int sqrt_samples = int(sqrt(samples_per_pixel));
+
+
     for (int j=start_line; j>=start_line - (lines - 1); --j) {
 
         for (int i=0; i<image_width; ++i) {
 
             color pixel_color(0, 0, 0);
+            color valid_sample_accum(0, 0, 0);  // Accumulate only valid samples
+            int valid_sample_count = 0;         // Track number of valid samples
+
+            for (int py = 0; py < sqrt_samples; ++py) {
+                for (int px = 0; px < sqrt_samples; ++px) {
+                    // Stratified sampling within the pixel
+                    auto u = (i + (px + random_double()) / sqrt_samples) / (image_width - 1);
+                    auto v = (j + (py + random_double()) / sqrt_samples) / (image_height - 1);
+                    ray r = cam.get_ray(u, v);
+                    color curr_sample = trace_ray(r, scene_ptr, max_depth);
+                    
+                    // Check if the sample is invalid (PDF = 0 or another condition)
+                    // trace_ray should catch invalid sampled, but we do another check here
+                    bool nan_or_inf = (
+                        !std::isfinite(curr_sample.x()) ||
+                        !std::isfinite(curr_sample.y()) ||
+                        !std::isfinite(curr_sample.z())
+                    );
+                    //if (curr_sample == INVALID_SAMPLE) {
+                    if (
+                        (curr_sample.x() == INVALID_SAMPLE.x() &&
+                        curr_sample.y() == INVALID_SAMPLE.y() &&
+                        curr_sample.z() == INVALID_SAMPLE.z()) ||
+                        nan_or_inf
+                    ) {
+                        // Replace the invalid sample with a "fake" sample that is the average
+                        if (valid_sample_count > 0) {
+                            curr_sample = valid_sample_accum / valid_sample_count;  // Use average of valid samples so far
+                        } else {
+                            curr_sample = color(0, 0, 0);  // No valid samples yet, use default color
+                        }
+                    } else {
+                        // Accumulate the valid sample and increase the count
+                        valid_sample_accum += curr_sample;
+                        valid_sample_count++;
+                    }
 
-            for (int s=0; s < samples_per_pixel; s++) {
-                auto u = (i + random_double()) / (image_width-1);
-                auto v = (j + random_double()) / (image_height-1);
-                ray r = cam.get_ray(u, v);
-                pixel_color += colorize_ray(r, scene_ptr, max_depth);
+                    pixel_color += curr_sample;
+                }
             }
 
             int buffer_index = j * image_width + i;
diff --git a/src/scene.cc b/src/scene.cc
index 8bf61d2..d87caf3 100644
--- a/src/scene.cc
+++ b/src/scene.cc
@@ -18,6 +18,18 @@ unsigned int Scene::add_primitive(std::shared_ptr<Primitive> prim) {
     return primID;
 }
 
+unsigned int Scene::add_volume(std::shared_ptr<Volume> vol) {
+    volumes.push_back(vol);
+    unsigned int primID = rtcAttachGeometry(rtc_scene, vol->geom);
+    rtcReleaseGeometry(vol->geom);
+    geom_map[primID] = vol;
+    return primID;
+}
+
+void Scene::add_physical_light(std::shared_ptr<Geometry> geom_ptr) {
+    physical_lights.push_back(geom_ptr);
+}
+
 unsigned int Scene::add_primitive_instance(std::shared_ptr<PrimitiveInstance> pi_ptr, RTCDevice device) {
     RTCGeometry instance_geom = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_INSTANCE);
     rtcSetGeometryInstancedScene(instance_geom, pi_ptr->instance_scene);
@@ -87,3 +99,8 @@ void add_triangle(RTCDevice device, RTCScene scene) {
     unsigned int triangleID = rtcAttachGeometry(scene, geom);
     rtcReleaseGeometry(geom);
 }
+
+void Scene::set_sky_colour(color bottom, color top) {
+    sky_bottom = bottom;
+    sky_top = top;
+}
\ No newline at end of file
diff --git a/src/util/vec3.cc b/src/util/vec3.cc
index 0e58088..cf5d6c6 100644
--- a/src/util/vec3.cc
+++ b/src/util/vec3.cc
@@ -88,6 +88,10 @@ vec3 cross(const vec3 &u, const vec3 &v) {
                 u.x() * v.y() - u.y() * v.x()};
 }
 
+vec3 mix(vec3 x, vec3 y, float a) {
+    return x * (1 - a) + y * a;
+}
+
 std::ostream& operator<<(std::ostream &out, const vec3 &v) {
     return out << v.x() << ' ' << v.y() << ' ' << v.z();
 }
@@ -119,6 +123,18 @@ vec3 random_in_unit_disk() {
     return rand_vec.unit_vector();
 }
 
+vec3 random_cosine_direction() {
+    auto r1 = random_double();
+    auto r2 = random_double();
+
+    auto phi = 2*pi*r1;
+    auto x = cos(phi)*sqrt(r2);
+    auto y = sin(phi)*sqrt(r2);
+    auto z = sqrt(1-r2);
+
+    return vec3(x, y, z);
+}
+
 vec3 reflect(const vec3& v, const vec3& n) { return v - 2*n * dot(v, n); }
 
 vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat) {
diff --git a/tests/0.1.5-showcase.csr b/tests/0.1.5-showcase.csr
new file mode 100644
index 0000000..a94ee3d
--- /dev/null
+++ b/tests/0.1.5-showcase.csr
@@ -0,0 +1,197 @@
+version 0.1.5
+
+Camera
+lookfrom 30 8 0
+lookat 0 4 0
+vup 0 1 0
+vfov 50
+aspect_ratio 16/9
+aperture 0.0001
+focus_dist 10.0
+
+# NEW FEATURE: We can set the sky colours!
+Sky
+top 1.0 0.9 0.5
+bottom 1.0 0.0 0.0
+
+# Scene Setup with red ground
+Material[Diffuse]
+id dark_brown
+texture no
+albedo 0.9 0.9 0.9
+roughness 0.0001
+
+Sphere
+id ground
+position 0 -2000 0
+material dark_brown
+radius 2000
+medium false
+
+
+# NEW FEATURE: Making a Volume
+# First, we create a sphere.
+Sphere
+id volume_sphere
+position 0 4 -8
+material dark_brown
+radius 4
+medium true
+
+# Then, we make a Medium (the properties of the volume)
+Medium
+id medium_example
+density 0.75
+albedo 0.23 0.1 0.6
+
+# Finally, we build the physical volume with the 
+# shape of the Sphere and properties of the Medium
+Volume
+id volume_example
+medium_id medium_example
+prim_id volume_sphere
+
+# NEW FEATURE: Making a Box with the new Diamond!
+
+# Lets start with making our new Dielectrics (like a diamond :p)
+Material[Dielectric]
+id cyan_diamond
+albedo 0.1 0.8 0.8
+eta 2.4
+roughness 0.0001
+sheen 0
+
+# And now we make the box with our dielectric material!
+Box
+id diamond_box
+position 0 4 4
+a 8 0 0
+b 0 -4 0
+c 0 0 8
+material cyan_diamond
+medium false
+
+# We can also make volumes with this one
+Box
+id volume_box
+position 0 8 4
+a 8 0 0
+b 0 -8 0
+c 0 0 8
+material cyan_diamond
+medium true
+
+# Then, we make another Medium (the properties of the volume)
+Medium
+id medium_box_example
+density 0.31 # lets make this one a little thinner!
+albedo 0.5 0.1 0.5
+
+# Now make a Medium out of this
+Volume
+id volume_box_example
+medium_id medium_box_example
+prim_id volume_box
+
+
+# We already covered the new dielectric model, but let's also make the two new types!
+
+# First, the new diffuse! (It's actually already visible on the ground, but lets make another).
+Material[Diffuse]
+id green_diffuse
+texture no
+albedo 0.13 0.9 0.1
+roughness 0.0001
+
+Sphere
+id green_sphere
+position -8 6 0
+material green_diffuse
+radius 6
+medium false
+
+# Now, lets check out the new metal!
+# Let's make a very shiny mirror...which involves a Metal and a Quad!
+Material[Metal]
+id mirror_example
+albedo 0.9 0.89 1.0 # to make it easy to see, we'll make it a little scratched up
+roughness 0.005
+
+Quad
+id mirror_quad
+position -20 18 -8
+u 15 0 -15
+v 0 -17 0
+material mirror_example
+
+# A very cool new feature is the Pixel Lambertian
+# You can finally add pixel art into your scenes! In fact, let's make Kylo Ren in front of
+# all the objects. And we'll make it a bit more cinematic...
+
+# Make the Kylo Ren
+Texture[Image]
+id kylo_texture
+transparency true
+path kylo.png
+
+Material[Diffuse]
+id kylo_material
+texture kylo_texture
+
+Quad
+id 0
+position 5 0 2.5
+u 0 0 -7.5
+v 0 7.5 0
+material kylo_material
+
+# Now, let's make it a little more cinematic...
+# Emissive
+Material[Emissive]
+id key_light
+rgb 0.6 0.6 0.1
+strength 10
+
+Sphere
+id key_sphere
+position -4000 0 4000
+material key_light
+radius 2000
+medium false
+
+# Let's also make a small floor underneath Kylo thats a MIX of multiple different materials!
+# NEW FEATURE: Mixture
+
+Material[Mixture]
+id mixture_example
+num 3
+mixed mirror_example
+weight 0.333
+mixed green_diffuse
+weight 0.333
+mixed cyan_diamond
+weight 0.334
+
+Quad
+id floor_quad
+position 1 0 2.5
+u 0 0 -7.5
+v 8 0 0
+material mixture_example
+
+# NEW FEATURE: Layered materials
+# Lets make a cool combination that isn't a mixture, but rather two layers of materials.
+# What if we want to give the diffuse a diamond exterior?
+
+Material[Layered]
+id layered_example
+top cyan_diamond
+bottom green_diffuse
+medium medium_example # it could also be 'no' if we didn't want to use one
+
+Sphere
+id layered_sphere
+position 7 2 -8
+material layered_example
+radius 2
+medium false
\ No newline at end of file
diff --git a/tests/cornell_box.csr b/tests/cornell_box.csr
index 1a3308d..1830eca 100644
--- a/tests/cornell_box.csr
+++ b/tests/cornell_box.csr
@@ -10,20 +10,23 @@ aspect_ratio 1/1
 aperture 0.0001
 focus_dist 10.0
 
-Material[Lambertian]
+Material[Diffuse]
 id red
 texture no
 albedo 0.65 0.05 0.05
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id white
 texture no
 albedo 0.73 0.73 0.73
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id green
 texture no
 albedo 0.12 0.45 0.15
+roughness 0.0001
 
 Material[Emissive]
 id lightmaterial
diff --git a/tests/dielectric_grid.csr b/tests/dielectric_grid.csr
new file mode 100644
index 0000000..4e935dd
--- /dev/null
+++ b/tests/dielectric_grid.csr
@@ -0,0 +1,3249 @@
+version 0.1.5
+
+Camera
+lookfrom -30 -10 0
+lookat 0 -10 0
+vup 0 1 0
+vfov 50
+aspect_ratio 16/9
+aperture 0.0001
+focus_dist 10.0
+
+Sky
+top 1.0 0.1 0.1
+bottom 0.1 1.0 0.1
+
+
+Material[Dielectric]
+id 1
+albedo 1 1 1
+eta 1
+roughness 0.0001
+sheen 0
+
+Sphere
+id 1
+position 0 0 -20
+material 1
+radius 1
+medium false
+
+Material[Dielectric]
+id 2
+albedo 1 1 1
+eta 1
+roughness 0.05
+sheen 0
+
+Sphere
+id 2
+position 0 0 -18
+material 2
+radius 1
+medium false
+
+Material[Dielectric]
+id 3
+albedo 1 1 1
+eta 1
+roughness 0.1
+sheen 0
+
+Sphere
+id 3
+position 0 0 -16
+material 3
+radius 1
+medium false
+
+Material[Dielectric]
+id 4
+albedo 1 1 1
+eta 1
+roughness 0.15
+sheen 0
+
+Sphere
+id 4
+position 0 0 -14
+material 4
+radius 1
+medium false
+
+Material[Dielectric]
+id 5
+albedo 1 1 1
+eta 1
+roughness 0.2
+sheen 0
+
+Sphere
+id 5
+position 0 0 -12
+material 5
+radius 1
+medium false
+
+Material[Dielectric]
+id 6
+albedo 1 1 1
+eta 1
+roughness 0.25
+sheen 0
+
+Sphere
+id 6
+position 0 0 -10
+material 6
+radius 1
+medium false
+
+Material[Dielectric]
+id 7
+albedo 1 1 1
+eta 1
+roughness 0.3
+sheen 0
+
+Sphere
+id 7
+position 0 0 -8
+material 7
+radius 1
+medium false
+
+Material[Dielectric]
+id 8
+albedo 1 1 1
+eta 1
+roughness 0.35
+sheen 0
+
+Sphere
+id 8
+position 0 0 -6
+material 8
+radius 1
+medium false
+
+Material[Dielectric]
+id 9
+albedo 1 1 1
+eta 1
+roughness 0.4
+sheen 0
+
+Sphere
+id 9
+position 0 0 -4
+material 9
+radius 1
+medium false
+
+Material[Dielectric]
+id 10
+albedo 1 1 1
+eta 1
+roughness 0.45
+sheen 0
+
+Sphere
+id 10
+position 0 0 -2
+material 10
+radius 1
+medium false
+
+Material[Dielectric]
+id 11
+albedo 1 1 1
+eta 1
+roughness 0.5
+sheen 0
+
+Sphere
+id 11
+position 0 0 0
+material 11
+radius 1
+medium false
+
+Material[Dielectric]
+id 12
+albedo 1 1 1
+eta 1
+roughness 0.55
+sheen 0
+
+Sphere
+id 12
+position 0 0 2
+material 12
+radius 1
+medium false
+
+Material[Dielectric]
+id 13
+albedo 1 1 1
+eta 1
+roughness 0.6
+sheen 0
+
+Sphere
+id 13
+position 0 0 4
+material 13
+radius 1
+medium false
+
+Material[Dielectric]
+id 14
+albedo 1 1 1
+eta 1
+roughness 0.65
+sheen 0
+
+Sphere
+id 14
+position 0 0 6
+material 14
+radius 1
+medium false
+
+Material[Dielectric]
+id 15
+albedo 1 1 1
+eta 1
+roughness 0.7
+sheen 0
+
+Sphere
+id 15
+position 0 0 8
+material 15
+radius 1
+medium false
+
+Material[Dielectric]
+id 16
+albedo 1 1 1
+eta 1
+roughness 0.75
+sheen 0
+
+Sphere
+id 16
+position 0 0 10
+material 16
+radius 1
+medium false
+
+Material[Dielectric]
+id 17
+albedo 1 1 1
+eta 1
+roughness 0.8
+sheen 0
+
+Sphere
+id 17
+position 0 0 12
+material 17
+radius 1
+medium false
+
+Material[Dielectric]
+id 18
+albedo 1 1 1
+eta 1
+roughness 0.85
+sheen 0
+
+Sphere
+id 18
+position 0 0 14
+material 18
+radius 1
+medium false
+
+Material[Dielectric]
+id 19
+albedo 1 1 1
+eta 1
+roughness 0.9
+sheen 0
+
+Sphere
+id 19
+position 0 0 16
+material 19
+radius 1
+medium false
+
+Material[Dielectric]
+id 20
+albedo 1 1 1
+eta 1
+roughness 0.95
+sheen 0
+
+Sphere
+id 20
+position 0 0 18
+material 20
+radius 1
+medium false
+
+Material[Dielectric]
+id 21
+albedo 1 1 1
+eta 1
+roughness 1
+sheen 0
+
+Sphere
+id 21
+position 0 0 20
+material 21
+radius 1
+medium false
+
+Material[Dielectric]
+id 22
+albedo 1 1 1
+eta 1.2
+roughness 0.0001
+sheen 0
+
+Sphere
+id 22
+position 0 -2 -20
+material 22
+radius 1
+medium false
+
+Material[Dielectric]
+id 23
+albedo 1 1 1
+eta 1.2
+roughness 0.05
+sheen 0
+
+Sphere
+id 23
+position 0 -2 -18
+material 23
+radius 1
+medium false
+
+Material[Dielectric]
+id 24
+albedo 1 1 1
+eta 1.2
+roughness 0.1
+sheen 0
+
+Sphere
+id 24
+position 0 -2 -16
+material 24
+radius 1
+medium false
+
+Material[Dielectric]
+id 25
+albedo 1 1 1
+eta 1.2
+roughness 0.15
+sheen 0
+
+Sphere
+id 25
+position 0 -2 -14
+material 25
+radius 1
+medium false
+
+Material[Dielectric]
+id 26
+albedo 1 1 1
+eta 1.2
+roughness 0.2
+sheen 0
+
+Sphere
+id 26
+position 0 -2 -12
+material 26
+radius 1
+medium false
+
+Material[Dielectric]
+id 27
+albedo 1 1 1
+eta 1.2
+roughness 0.25
+sheen 0
+
+Sphere
+id 27
+position 0 -2 -10
+material 27
+radius 1
+medium false
+
+Material[Dielectric]
+id 28
+albedo 1 1 1
+eta 1.2
+roughness 0.3
+sheen 0
+
+Sphere
+id 28
+position 0 -2 -8
+material 28
+radius 1
+medium false
+
+Material[Dielectric]
+id 29
+albedo 1 1 1
+eta 1.2
+roughness 0.35
+sheen 0
+
+Sphere
+id 29
+position 0 -2 -6
+material 29
+radius 1
+medium false
+
+Material[Dielectric]
+id 30
+albedo 1 1 1
+eta 1.2
+roughness 0.4
+sheen 0
+
+Sphere
+id 30
+position 0 -2 -4
+material 30
+radius 1
+medium false
+
+Material[Dielectric]
+id 31
+albedo 1 1 1
+eta 1.2
+roughness 0.45
+sheen 0
+
+Sphere
+id 31
+position 0 -2 -2
+material 31
+radius 1
+medium false
+
+Material[Dielectric]
+id 32
+albedo 1 1 1
+eta 1.2
+roughness 0.5
+sheen 0
+
+Sphere
+id 32
+position 0 -2 0
+material 32
+radius 1
+medium false
+
+Material[Dielectric]
+id 33
+albedo 1 1 1
+eta 1.2
+roughness 0.55
+sheen 0
+
+Sphere
+id 33
+position 0 -2 2
+material 33
+radius 1
+medium false
+
+Material[Dielectric]
+id 34
+albedo 1 1 1
+eta 1.2
+roughness 0.6
+sheen 0
+
+Sphere
+id 34
+position 0 -2 4
+material 34
+radius 1
+medium false
+
+Material[Dielectric]
+id 35
+albedo 1 1 1
+eta 1.2
+roughness 0.65
+sheen 0
+
+Sphere
+id 35
+position 0 -2 6
+material 35
+radius 1
+medium false
+
+Material[Dielectric]
+id 36
+albedo 1 1 1
+eta 1.2
+roughness 0.7
+sheen 0
+
+Sphere
+id 36
+position 0 -2 8
+material 36
+radius 1
+medium false
+
+Material[Dielectric]
+id 37
+albedo 1 1 1
+eta 1.2
+roughness 0.75
+sheen 0
+
+Sphere
+id 37
+position 0 -2 10
+material 37
+radius 1
+medium false
+
+Material[Dielectric]
+id 38
+albedo 1 1 1
+eta 1.2
+roughness 0.8
+sheen 0
+
+Sphere
+id 38
+position 0 -2 12
+material 38
+radius 1
+medium false
+
+Material[Dielectric]
+id 39
+albedo 1 1 1
+eta 1.2
+roughness 0.85
+sheen 0
+
+Sphere
+id 39
+position 0 -2 14
+material 39
+radius 1
+medium false
+
+Material[Dielectric]
+id 40
+albedo 1 1 1
+eta 1.2
+roughness 0.9
+sheen 0
+
+Sphere
+id 40
+position 0 -2 16
+material 40
+radius 1
+medium false
+
+Material[Dielectric]
+id 41
+albedo 1 1 1
+eta 1.2
+roughness 0.95
+sheen 0
+
+Sphere
+id 41
+position 0 -2 18
+material 41
+radius 1
+medium false
+
+Material[Dielectric]
+id 42
+albedo 1 1 1
+eta 1.2
+roughness 1
+sheen 0
+
+Sphere
+id 42
+position 0 -2 20
+material 42
+radius 1
+medium false
+
+Material[Dielectric]
+id 43
+albedo 1 1 1
+eta 1.4
+roughness 0.0001
+sheen 0
+
+Sphere
+id 43
+position 0 -4 -20
+material 43
+radius 1
+medium false
+
+Material[Dielectric]
+id 44
+albedo 1 1 1
+eta 1.4
+roughness 0.05
+sheen 0
+
+Sphere
+id 44
+position 0 -4 -18
+material 44
+radius 1
+medium false
+
+Material[Dielectric]
+id 45
+albedo 1 1 1
+eta 1.4
+roughness 0.1
+sheen 0
+
+Sphere
+id 45
+position 0 -4 -16
+material 45
+radius 1
+medium false
+
+Material[Dielectric]
+id 46
+albedo 1 1 1
+eta 1.4
+roughness 0.15
+sheen 0
+
+Sphere
+id 46
+position 0 -4 -14
+material 46
+radius 1
+medium false
+
+Material[Dielectric]
+id 47
+albedo 1 1 1
+eta 1.4
+roughness 0.2
+sheen 0
+
+Sphere
+id 47
+position 0 -4 -12
+material 47
+radius 1
+medium false
+
+Material[Dielectric]
+id 48
+albedo 1 1 1
+eta 1.4
+roughness 0.25
+sheen 0
+
+Sphere
+id 48
+position 0 -4 -10
+material 48
+radius 1
+medium false
+
+Material[Dielectric]
+id 49
+albedo 1 1 1
+eta 1.4
+roughness 0.3
+sheen 0
+
+Sphere
+id 49
+position 0 -4 -8
+material 49
+radius 1
+medium false
+
+Material[Dielectric]
+id 50
+albedo 1 1 1
+eta 1.4
+roughness 0.35
+sheen 0
+
+Sphere
+id 50
+position 0 -4 -6
+material 50
+radius 1
+medium false
+
+Material[Dielectric]
+id 51
+albedo 1 1 1
+eta 1.4
+roughness 0.4
+sheen 0
+
+Sphere
+id 51
+position 0 -4 -4
+material 51
+radius 1
+medium false
+
+Material[Dielectric]
+id 52
+albedo 1 1 1
+eta 1.4
+roughness 0.45
+sheen 0
+
+Sphere
+id 52
+position 0 -4 -2
+material 52
+radius 1
+medium false
+
+Material[Dielectric]
+id 53
+albedo 1 1 1
+eta 1.4
+roughness 0.5
+sheen 0
+
+Sphere
+id 53
+position 0 -4 0
+material 53
+radius 1
+medium false
+
+Material[Dielectric]
+id 54
+albedo 1 1 1
+eta 1.4
+roughness 0.55
+sheen 0
+
+Sphere
+id 54
+position 0 -4 2
+material 54
+radius 1
+medium false
+
+Material[Dielectric]
+id 55
+albedo 1 1 1
+eta 1.4
+roughness 0.6
+sheen 0
+
+Sphere
+id 55
+position 0 -4 4
+material 55
+radius 1
+medium false
+
+Material[Dielectric]
+id 56
+albedo 1 1 1
+eta 1.4
+roughness 0.65
+sheen 0
+
+Sphere
+id 56
+position 0 -4 6
+material 56
+radius 1
+medium false
+
+Material[Dielectric]
+id 57
+albedo 1 1 1
+eta 1.4
+roughness 0.7
+sheen 0
+
+Sphere
+id 57
+position 0 -4 8
+material 57
+radius 1
+medium false
+
+Material[Dielectric]
+id 58
+albedo 1 1 1
+eta 1.4
+roughness 0.75
+sheen 0
+
+Sphere
+id 58
+position 0 -4 10
+material 58
+radius 1
+medium false
+
+Material[Dielectric]
+id 59
+albedo 1 1 1
+eta 1.4
+roughness 0.8
+sheen 0
+
+Sphere
+id 59
+position 0 -4 12
+material 59
+radius 1
+medium false
+
+Material[Dielectric]
+id 60
+albedo 1 1 1
+eta 1.4
+roughness 0.85
+sheen 0
+
+Sphere
+id 60
+position 0 -4 14
+material 60
+radius 1
+medium false
+
+Material[Dielectric]
+id 61
+albedo 1 1 1
+eta 1.4
+roughness 0.9
+sheen 0
+
+Sphere
+id 61
+position 0 -4 16
+material 61
+radius 1
+medium false
+
+Material[Dielectric]
+id 62
+albedo 1 1 1
+eta 1.4
+roughness 0.95
+sheen 0
+
+Sphere
+id 62
+position 0 -4 18
+material 62
+radius 1
+medium false
+
+Material[Dielectric]
+id 63
+albedo 1 1 1
+eta 1.4
+roughness 1
+sheen 0
+
+Sphere
+id 63
+position 0 -4 20
+material 63
+radius 1
+medium false
+
+Material[Dielectric]
+id 64
+albedo 1 1 1
+eta 1.6
+roughness 0.0001
+sheen 0
+
+Sphere
+id 64
+position 0 -6 -20
+material 64
+radius 1
+medium false
+
+Material[Dielectric]
+id 65
+albedo 1 1 1
+eta 1.6
+roughness 0.05
+sheen 0
+
+Sphere
+id 65
+position 0 -6 -18
+material 65
+radius 1
+medium false
+
+Material[Dielectric]
+id 66
+albedo 1 1 1
+eta 1.6
+roughness 0.1
+sheen 0
+
+Sphere
+id 66
+position 0 -6 -16
+material 66
+radius 1
+medium false
+
+Material[Dielectric]
+id 67
+albedo 1 1 1
+eta 1.6
+roughness 0.15
+sheen 0
+
+Sphere
+id 67
+position 0 -6 -14
+material 67
+radius 1
+medium false
+
+Material[Dielectric]
+id 68
+albedo 1 1 1
+eta 1.6
+roughness 0.2
+sheen 0
+
+Sphere
+id 68
+position 0 -6 -12
+material 68
+radius 1
+medium false
+
+Material[Dielectric]
+id 69
+albedo 1 1 1
+eta 1.6
+roughness 0.25
+sheen 0
+
+Sphere
+id 69
+position 0 -6 -10
+material 69
+radius 1
+medium false
+
+Material[Dielectric]
+id 70
+albedo 1 1 1
+eta 1.6
+roughness 0.3
+sheen 0
+
+Sphere
+id 70
+position 0 -6 -8
+material 70
+radius 1
+medium false
+
+Material[Dielectric]
+id 71
+albedo 1 1 1
+eta 1.6
+roughness 0.35
+sheen 0
+
+Sphere
+id 71
+position 0 -6 -6
+material 71
+radius 1
+medium false
+
+Material[Dielectric]
+id 72
+albedo 1 1 1
+eta 1.6
+roughness 0.4
+sheen 0
+
+Sphere
+id 72
+position 0 -6 -4
+material 72
+radius 1
+medium false
+
+Material[Dielectric]
+id 73
+albedo 1 1 1
+eta 1.6
+roughness 0.45
+sheen 0
+
+Sphere
+id 73
+position 0 -6 -2
+material 73
+radius 1
+medium false
+
+Material[Dielectric]
+id 74
+albedo 1 1 1
+eta 1.6
+roughness 0.5
+sheen 0
+
+Sphere
+id 74
+position 0 -6 0
+material 74
+radius 1
+medium false
+
+Material[Dielectric]
+id 75
+albedo 1 1 1
+eta 1.6
+roughness 0.55
+sheen 0
+
+Sphere
+id 75
+position 0 -6 2
+material 75
+radius 1
+medium false
+
+Material[Dielectric]
+id 76
+albedo 1 1 1
+eta 1.6
+roughness 0.6
+sheen 0
+
+Sphere
+id 76
+position 0 -6 4
+material 76
+radius 1
+medium false
+
+Material[Dielectric]
+id 77
+albedo 1 1 1
+eta 1.6
+roughness 0.65
+sheen 0
+
+Sphere
+id 77
+position 0 -6 6
+material 77
+radius 1
+medium false
+
+Material[Dielectric]
+id 78
+albedo 1 1 1
+eta 1.6
+roughness 0.7
+sheen 0
+
+Sphere
+id 78
+position 0 -6 8
+material 78
+radius 1
+medium false
+
+Material[Dielectric]
+id 79
+albedo 1 1 1
+eta 1.6
+roughness 0.75
+sheen 0
+
+Sphere
+id 79
+position 0 -6 10
+material 79
+radius 1
+medium false
+
+Material[Dielectric]
+id 80
+albedo 1 1 1
+eta 1.6
+roughness 0.8
+sheen 0
+
+Sphere
+id 80
+position 0 -6 12
+material 80
+radius 1
+medium false
+
+Material[Dielectric]
+id 81
+albedo 1 1 1
+eta 1.6
+roughness 0.85
+sheen 0
+
+Sphere
+id 81
+position 0 -6 14
+material 81
+radius 1
+medium false
+
+Material[Dielectric]
+id 82
+albedo 1 1 1
+eta 1.6
+roughness 0.9
+sheen 0
+
+Sphere
+id 82
+position 0 -6 16
+material 82
+radius 1
+medium false
+
+Material[Dielectric]
+id 83
+albedo 1 1 1
+eta 1.6
+roughness 0.95
+sheen 0
+
+Sphere
+id 83
+position 0 -6 18
+material 83
+radius 1
+medium false
+
+Material[Dielectric]
+id 84
+albedo 1 1 1
+eta 1.6
+roughness 1
+sheen 0
+
+Sphere
+id 84
+position 0 -6 20
+material 84
+radius 1
+medium false
+
+Material[Dielectric]
+id 85
+albedo 1 1 1
+eta 1.8
+roughness 0.0001
+sheen 0
+
+Sphere
+id 85
+position 0 -8 -20
+material 85
+radius 1
+medium false
+
+Material[Dielectric]
+id 86
+albedo 1 1 1
+eta 1.8
+roughness 0.05
+sheen 0
+
+Sphere
+id 86
+position 0 -8 -18
+material 86
+radius 1
+medium false
+
+Material[Dielectric]
+id 87
+albedo 1 1 1
+eta 1.8
+roughness 0.1
+sheen 0
+
+Sphere
+id 87
+position 0 -8 -16
+material 87
+radius 1
+medium false
+
+Material[Dielectric]
+id 88
+albedo 1 1 1
+eta 1.8
+roughness 0.15
+sheen 0
+
+Sphere
+id 88
+position 0 -8 -14
+material 88
+radius 1
+medium false
+
+Material[Dielectric]
+id 89
+albedo 1 1 1
+eta 1.8
+roughness 0.2
+sheen 0
+
+Sphere
+id 89
+position 0 -8 -12
+material 89
+radius 1
+medium false
+
+Material[Dielectric]
+id 90
+albedo 1 1 1
+eta 1.8
+roughness 0.25
+sheen 0
+
+Sphere
+id 90
+position 0 -8 -10
+material 90
+radius 1
+medium false
+
+Material[Dielectric]
+id 91
+albedo 1 1 1
+eta 1.8
+roughness 0.3
+sheen 0
+
+Sphere
+id 91
+position 0 -8 -8
+material 91
+radius 1
+medium false
+
+Material[Dielectric]
+id 92
+albedo 1 1 1
+eta 1.8
+roughness 0.35
+sheen 0
+
+Sphere
+id 92
+position 0 -8 -6
+material 92
+radius 1
+medium false
+
+Material[Dielectric]
+id 93
+albedo 1 1 1
+eta 1.8
+roughness 0.4
+sheen 0
+
+Sphere
+id 93
+position 0 -8 -4
+material 93
+radius 1
+medium false
+
+Material[Dielectric]
+id 94
+albedo 1 1 1
+eta 1.8
+roughness 0.45
+sheen 0
+
+Sphere
+id 94
+position 0 -8 -2
+material 94
+radius 1
+medium false
+
+Material[Dielectric]
+id 95
+albedo 1 1 1
+eta 1.8
+roughness 0.5
+sheen 0
+
+Sphere
+id 95
+position 0 -8 0
+material 95
+radius 1
+medium false
+
+Material[Dielectric]
+id 96
+albedo 1 1 1
+eta 1.8
+roughness 0.55
+sheen 0
+
+Sphere
+id 96
+position 0 -8 2
+material 96
+radius 1
+medium false
+
+Material[Dielectric]
+id 97
+albedo 1 1 1
+eta 1.8
+roughness 0.6
+sheen 0
+
+Sphere
+id 97
+position 0 -8 4
+material 97
+radius 1
+medium false
+
+Material[Dielectric]
+id 98
+albedo 1 1 1
+eta 1.8
+roughness 0.65
+sheen 0
+
+Sphere
+id 98
+position 0 -8 6
+material 98
+radius 1
+medium false
+
+Material[Dielectric]
+id 99
+albedo 1 1 1
+eta 1.8
+roughness 0.7
+sheen 0
+
+Sphere
+id 99
+position 0 -8 8
+material 99
+radius 1
+medium false
+
+Material[Dielectric]
+id 100
+albedo 1 1 1
+eta 1.8
+roughness 0.75
+sheen 0
+
+Sphere
+id 100
+position 0 -8 10
+material 100
+radius 1
+medium false
+
+Material[Dielectric]
+id 101
+albedo 1 1 1
+eta 1.8
+roughness 0.8
+sheen 0
+
+Sphere
+id 101
+position 0 -8 12
+material 101
+radius 1
+medium false
+
+Material[Dielectric]
+id 102
+albedo 1 1 1
+eta 1.8
+roughness 0.85
+sheen 0
+
+Sphere
+id 102
+position 0 -8 14
+material 102
+radius 1
+medium false
+
+Material[Dielectric]
+id 103
+albedo 1 1 1
+eta 1.8
+roughness 0.9
+sheen 0
+
+Sphere
+id 103
+position 0 -8 16
+material 103
+radius 1
+medium false
+
+Material[Dielectric]
+id 104
+albedo 1 1 1
+eta 1.8
+roughness 0.95
+sheen 0
+
+Sphere
+id 104
+position 0 -8 18
+material 104
+radius 1
+medium false
+
+Material[Dielectric]
+id 105
+albedo 1 1 1
+eta 1.8
+roughness 1
+sheen 0
+
+Sphere
+id 105
+position 0 -8 20
+material 105
+radius 1
+medium false
+
+Material[Dielectric]
+id 106
+albedo 1 1 1
+eta 2
+roughness 0.0001
+sheen 0
+
+Sphere
+id 106
+position 0 -10 -20
+material 106
+radius 1
+medium false
+
+Material[Dielectric]
+id 107
+albedo 1 1 1
+eta 2
+roughness 0.05
+sheen 0
+
+Sphere
+id 107
+position 0 -10 -18
+material 107
+radius 1
+medium false
+
+Material[Dielectric]
+id 108
+albedo 1 1 1
+eta 2
+roughness 0.1
+sheen 0
+
+Sphere
+id 108
+position 0 -10 -16
+material 108
+radius 1
+medium false
+
+Material[Dielectric]
+id 109
+albedo 1 1 1
+eta 2
+roughness 0.15
+sheen 0
+
+Sphere
+id 109
+position 0 -10 -14
+material 109
+radius 1
+medium false
+
+Material[Dielectric]
+id 110
+albedo 1 1 1
+eta 2
+roughness 0.2
+sheen 0
+
+Sphere
+id 110
+position 0 -10 -12
+material 110
+radius 1
+medium false
+
+Material[Dielectric]
+id 111
+albedo 1 1 1
+eta 2
+roughness 0.25
+sheen 0
+
+Sphere
+id 111
+position 0 -10 -10
+material 111
+radius 1
+medium false
+
+Material[Dielectric]
+id 112
+albedo 1 1 1
+eta 2
+roughness 0.3
+sheen 0
+
+Sphere
+id 112
+position 0 -10 -8
+material 112
+radius 1
+medium false
+
+Material[Dielectric]
+id 113
+albedo 1 1 1
+eta 2
+roughness 0.35
+sheen 0
+
+Sphere
+id 113
+position 0 -10 -6
+material 113
+radius 1
+medium false
+
+Material[Dielectric]
+id 114
+albedo 1 1 1
+eta 2
+roughness 0.4
+sheen 0
+
+Sphere
+id 114
+position 0 -10 -4
+material 114
+radius 1
+medium false
+
+Material[Dielectric]
+id 115
+albedo 1 1 1
+eta 2
+roughness 0.45
+sheen 0
+
+Sphere
+id 115
+position 0 -10 -2
+material 115
+radius 1
+medium false
+
+Material[Dielectric]
+id 116
+albedo 1 1 1
+eta 2
+roughness 0.5
+sheen 0
+
+Sphere
+id 116
+position 0 -10 0
+material 116
+radius 1
+medium false
+
+Material[Dielectric]
+id 117
+albedo 1 1 1
+eta 2
+roughness 0.55
+sheen 0
+
+Sphere
+id 117
+position 0 -10 2
+material 117
+radius 1
+medium false
+
+Material[Dielectric]
+id 118
+albedo 1 1 1
+eta 2
+roughness 0.6
+sheen 0
+
+Sphere
+id 118
+position 0 -10 4
+material 118
+radius 1
+medium false
+
+Material[Dielectric]
+id 119
+albedo 1 1 1
+eta 2
+roughness 0.65
+sheen 0
+
+Sphere
+id 119
+position 0 -10 6
+material 119
+radius 1
+medium false
+
+Material[Dielectric]
+id 120
+albedo 1 1 1
+eta 2
+roughness 0.7
+sheen 0
+
+Sphere
+id 120
+position 0 -10 8
+material 120
+radius 1
+medium false
+
+Material[Dielectric]
+id 121
+albedo 1 1 1
+eta 2
+roughness 0.75
+sheen 0
+
+Sphere
+id 121
+position 0 -10 10
+material 121
+radius 1
+medium false
+
+Material[Dielectric]
+id 122
+albedo 1 1 1
+eta 2
+roughness 0.8
+sheen 0
+
+Sphere
+id 122
+position 0 -10 12
+material 122
+radius 1
+medium false
+
+Material[Dielectric]
+id 123
+albedo 1 1 1
+eta 2
+roughness 0.85
+sheen 0
+
+Sphere
+id 123
+position 0 -10 14
+material 123
+radius 1
+medium false
+
+Material[Dielectric]
+id 124
+albedo 1 1 1
+eta 2
+roughness 0.9
+sheen 0
+
+Sphere
+id 124
+position 0 -10 16
+material 124
+radius 1
+medium false
+
+Material[Dielectric]
+id 125
+albedo 1 1 1
+eta 2
+roughness 0.95
+sheen 0
+
+Sphere
+id 125
+position 0 -10 18
+material 125
+radius 1
+medium false
+
+Material[Dielectric]
+id 126
+albedo 1 1 1
+eta 2
+roughness 1
+sheen 0
+
+Sphere
+id 126
+position 0 -10 20
+material 126
+radius 1
+medium false
+
+Material[Dielectric]
+id 127
+albedo 1 1 1
+eta 2.2
+roughness 0.0001
+sheen 0
+
+Sphere
+id 127
+position 0 -12 -20
+material 127
+radius 1
+medium false
+
+Material[Dielectric]
+id 128
+albedo 1 1 1
+eta 2.2
+roughness 0.05
+sheen 0
+
+Sphere
+id 128
+position 0 -12 -18
+material 128
+radius 1
+medium false
+
+Material[Dielectric]
+id 129
+albedo 1 1 1
+eta 2.2
+roughness 0.1
+sheen 0
+
+Sphere
+id 129
+position 0 -12 -16
+material 129
+radius 1
+medium false
+
+Material[Dielectric]
+id 130
+albedo 1 1 1
+eta 2.2
+roughness 0.15
+sheen 0
+
+Sphere
+id 130
+position 0 -12 -14
+material 130
+radius 1
+medium false
+
+Material[Dielectric]
+id 131
+albedo 1 1 1
+eta 2.2
+roughness 0.2
+sheen 0
+
+Sphere
+id 131
+position 0 -12 -12
+material 131
+radius 1
+medium false
+
+Material[Dielectric]
+id 132
+albedo 1 1 1
+eta 2.2
+roughness 0.25
+sheen 0
+
+Sphere
+id 132
+position 0 -12 -10
+material 132
+radius 1
+medium false
+
+Material[Dielectric]
+id 133
+albedo 1 1 1
+eta 2.2
+roughness 0.3
+sheen 0
+
+Sphere
+id 133
+position 0 -12 -8
+material 133
+radius 1
+medium false
+
+Material[Dielectric]
+id 134
+albedo 1 1 1
+eta 2.2
+roughness 0.35
+sheen 0
+
+Sphere
+id 134
+position 0 -12 -6
+material 134
+radius 1
+medium false
+
+Material[Dielectric]
+id 135
+albedo 1 1 1
+eta 2.2
+roughness 0.4
+sheen 0
+
+Sphere
+id 135
+position 0 -12 -4
+material 135
+radius 1
+medium false
+
+Material[Dielectric]
+id 136
+albedo 1 1 1
+eta 2.2
+roughness 0.45
+sheen 0
+
+Sphere
+id 136
+position 0 -12 -2
+material 136
+radius 1
+medium false
+
+Material[Dielectric]
+id 137
+albedo 1 1 1
+eta 2.2
+roughness 0.5
+sheen 0
+
+Sphere
+id 137
+position 0 -12 0
+material 137
+radius 1
+medium false
+
+Material[Dielectric]
+id 138
+albedo 1 1 1
+eta 2.2
+roughness 0.55
+sheen 0
+
+Sphere
+id 138
+position 0 -12 2
+material 138
+radius 1
+medium false
+
+Material[Dielectric]
+id 139
+albedo 1 1 1
+eta 2.2
+roughness 0.6
+sheen 0
+
+Sphere
+id 139
+position 0 -12 4
+material 139
+radius 1
+medium false
+
+Material[Dielectric]
+id 140
+albedo 1 1 1
+eta 2.2
+roughness 0.65
+sheen 0
+
+Sphere
+id 140
+position 0 -12 6
+material 140
+radius 1
+medium false
+
+Material[Dielectric]
+id 141
+albedo 1 1 1
+eta 2.2
+roughness 0.7
+sheen 0
+
+Sphere
+id 141
+position 0 -12 8
+material 141
+radius 1
+medium false
+
+Material[Dielectric]
+id 142
+albedo 1 1 1
+eta 2.2
+roughness 0.75
+sheen 0
+
+Sphere
+id 142
+position 0 -12 10
+material 142
+radius 1
+medium false
+
+Material[Dielectric]
+id 143
+albedo 1 1 1
+eta 2.2
+roughness 0.8
+sheen 0
+
+Sphere
+id 143
+position 0 -12 12
+material 143
+radius 1
+medium false
+
+Material[Dielectric]
+id 144
+albedo 1 1 1
+eta 2.2
+roughness 0.85
+sheen 0
+
+Sphere
+id 144
+position 0 -12 14
+material 144
+radius 1
+medium false
+
+Material[Dielectric]
+id 145
+albedo 1 1 1
+eta 2.2
+roughness 0.9
+sheen 0
+
+Sphere
+id 145
+position 0 -12 16
+material 145
+radius 1
+medium false
+
+Material[Dielectric]
+id 146
+albedo 1 1 1
+eta 2.2
+roughness 0.95
+sheen 0
+
+Sphere
+id 146
+position 0 -12 18
+material 146
+radius 1
+medium false
+
+Material[Dielectric]
+id 147
+albedo 1 1 1
+eta 2.2
+roughness 1
+sheen 0
+
+Sphere
+id 147
+position 0 -12 20
+material 147
+radius 1
+medium false
+
+Material[Dielectric]
+id 148
+albedo 1 1 1
+eta 2.4
+roughness 0.0001
+sheen 0
+
+Sphere
+id 148
+position 0 -14 -20
+material 148
+radius 1
+medium false
+
+Material[Dielectric]
+id 149
+albedo 1 1 1
+eta 2.4
+roughness 0.05
+sheen 0
+
+Sphere
+id 149
+position 0 -14 -18
+material 149
+radius 1
+medium false
+
+Material[Dielectric]
+id 150
+albedo 1 1 1
+eta 2.4
+roughness 0.1
+sheen 0
+
+Sphere
+id 150
+position 0 -14 -16
+material 150
+radius 1
+medium false
+
+Material[Dielectric]
+id 151
+albedo 1 1 1
+eta 2.4
+roughness 0.15
+sheen 0
+
+Sphere
+id 151
+position 0 -14 -14
+material 151
+radius 1
+medium false
+
+Material[Dielectric]
+id 152
+albedo 1 1 1
+eta 2.4
+roughness 0.2
+sheen 0
+
+Sphere
+id 152
+position 0 -14 -12
+material 152
+radius 1
+medium false
+
+Material[Dielectric]
+id 153
+albedo 1 1 1
+eta 2.4
+roughness 0.25
+sheen 0
+
+Sphere
+id 153
+position 0 -14 -10
+material 153
+radius 1
+medium false
+
+Material[Dielectric]
+id 154
+albedo 1 1 1
+eta 2.4
+roughness 0.3
+sheen 0
+
+Sphere
+id 154
+position 0 -14 -8
+material 154
+radius 1
+medium false
+
+Material[Dielectric]
+id 155
+albedo 1 1 1
+eta 2.4
+roughness 0.35
+sheen 0
+
+Sphere
+id 155
+position 0 -14 -6
+material 155
+radius 1
+medium false
+
+Material[Dielectric]
+id 156
+albedo 1 1 1
+eta 2.4
+roughness 0.4
+sheen 0
+
+Sphere
+id 156
+position 0 -14 -4
+material 156
+radius 1
+medium false
+
+Material[Dielectric]
+id 157
+albedo 1 1 1
+eta 2.4
+roughness 0.45
+sheen 0
+
+Sphere
+id 157
+position 0 -14 -2
+material 157
+radius 1
+medium false
+
+Material[Dielectric]
+id 158
+albedo 1 1 1
+eta 2.4
+roughness 0.5
+sheen 0
+
+Sphere
+id 158
+position 0 -14 0
+material 158
+radius 1
+medium false
+
+Material[Dielectric]
+id 159
+albedo 1 1 1
+eta 2.4
+roughness 0.55
+sheen 0
+
+Sphere
+id 159
+position 0 -14 2
+material 159
+radius 1
+medium false
+
+Material[Dielectric]
+id 160
+albedo 1 1 1
+eta 2.4
+roughness 0.6
+sheen 0
+
+Sphere
+id 160
+position 0 -14 4
+material 160
+radius 1
+medium false
+
+Material[Dielectric]
+id 161
+albedo 1 1 1
+eta 2.4
+roughness 0.65
+sheen 0
+
+Sphere
+id 161
+position 0 -14 6
+material 161
+radius 1
+medium false
+
+Material[Dielectric]
+id 162
+albedo 1 1 1
+eta 2.4
+roughness 0.7
+sheen 0
+
+Sphere
+id 162
+position 0 -14 8
+material 162
+radius 1
+medium false
+
+Material[Dielectric]
+id 163
+albedo 1 1 1
+eta 2.4
+roughness 0.75
+sheen 0
+
+Sphere
+id 163
+position 0 -14 10
+material 163
+radius 1
+medium false
+
+Material[Dielectric]
+id 164
+albedo 1 1 1
+eta 2.4
+roughness 0.8
+sheen 0
+
+Sphere
+id 164
+position 0 -14 12
+material 164
+radius 1
+medium false
+
+Material[Dielectric]
+id 165
+albedo 1 1 1
+eta 2.4
+roughness 0.85
+sheen 0
+
+Sphere
+id 165
+position 0 -14 14
+material 165
+radius 1
+medium false
+
+Material[Dielectric]
+id 166
+albedo 1 1 1
+eta 2.4
+roughness 0.9
+sheen 0
+
+Sphere
+id 166
+position 0 -14 16
+material 166
+radius 1
+medium false
+
+Material[Dielectric]
+id 167
+albedo 1 1 1
+eta 2.4
+roughness 0.95
+sheen 0
+
+Sphere
+id 167
+position 0 -14 18
+material 167
+radius 1
+medium false
+
+Material[Dielectric]
+id 168
+albedo 1 1 1
+eta 2.4
+roughness 1
+sheen 0
+
+Sphere
+id 168
+position 0 -14 20
+material 168
+radius 1
+medium false
+
+Material[Dielectric]
+id 169
+albedo 1 1 1
+eta 2.6
+roughness 0.0001
+sheen 0
+
+Sphere
+id 169
+position 0 -16 -20
+material 169
+radius 1
+medium false
+
+Material[Dielectric]
+id 170
+albedo 1 1 1
+eta 2.6
+roughness 0.05
+sheen 0
+
+Sphere
+id 170
+position 0 -16 -18
+material 170
+radius 1
+medium false
+
+Material[Dielectric]
+id 171
+albedo 1 1 1
+eta 2.6
+roughness 0.1
+sheen 0
+
+Sphere
+id 171
+position 0 -16 -16
+material 171
+radius 1
+medium false
+
+Material[Dielectric]
+id 172
+albedo 1 1 1
+eta 2.6
+roughness 0.15
+sheen 0
+
+Sphere
+id 172
+position 0 -16 -14
+material 172
+radius 1
+medium false
+
+Material[Dielectric]
+id 173
+albedo 1 1 1
+eta 2.6
+roughness 0.2
+sheen 0
+
+Sphere
+id 173
+position 0 -16 -12
+material 173
+radius 1
+medium false
+
+Material[Dielectric]
+id 174
+albedo 1 1 1
+eta 2.6
+roughness 0.25
+sheen 0
+
+Sphere
+id 174
+position 0 -16 -10
+material 174
+radius 1
+medium false
+
+Material[Dielectric]
+id 175
+albedo 1 1 1
+eta 2.6
+roughness 0.3
+sheen 0
+
+Sphere
+id 175
+position 0 -16 -8
+material 175
+radius 1
+medium false
+
+Material[Dielectric]
+id 176
+albedo 1 1 1
+eta 2.6
+roughness 0.35
+sheen 0
+
+Sphere
+id 176
+position 0 -16 -6
+material 176
+radius 1
+medium false
+
+Material[Dielectric]
+id 177
+albedo 1 1 1
+eta 2.6
+roughness 0.4
+sheen 0
+
+Sphere
+id 177
+position 0 -16 -4
+material 177
+radius 1
+medium false
+
+Material[Dielectric]
+id 178
+albedo 1 1 1
+eta 2.6
+roughness 0.45
+sheen 0
+
+Sphere
+id 178
+position 0 -16 -2
+material 178
+radius 1
+medium false
+
+Material[Dielectric]
+id 179
+albedo 1 1 1
+eta 2.6
+roughness 0.5
+sheen 0
+
+Sphere
+id 179
+position 0 -16 0
+material 179
+radius 1
+medium false
+
+Material[Dielectric]
+id 180
+albedo 1 1 1
+eta 2.6
+roughness 0.55
+sheen 0
+
+Sphere
+id 180
+position 0 -16 2
+material 180
+radius 1
+medium false
+
+Material[Dielectric]
+id 181
+albedo 1 1 1
+eta 2.6
+roughness 0.6
+sheen 0
+
+Sphere
+id 181
+position 0 -16 4
+material 181
+radius 1
+medium false
+
+Material[Dielectric]
+id 182
+albedo 1 1 1
+eta 2.6
+roughness 0.65
+sheen 0
+
+Sphere
+id 182
+position 0 -16 6
+material 182
+radius 1
+medium false
+
+Material[Dielectric]
+id 183
+albedo 1 1 1
+eta 2.6
+roughness 0.7
+sheen 0
+
+Sphere
+id 183
+position 0 -16 8
+material 183
+radius 1
+medium false
+
+Material[Dielectric]
+id 184
+albedo 1 1 1
+eta 2.6
+roughness 0.75
+sheen 0
+
+Sphere
+id 184
+position 0 -16 10
+material 184
+radius 1
+medium false
+
+Material[Dielectric]
+id 185
+albedo 1 1 1
+eta 2.6
+roughness 0.8
+sheen 0
+
+Sphere
+id 185
+position 0 -16 12
+material 185
+radius 1
+medium false
+
+Material[Dielectric]
+id 186
+albedo 1 1 1
+eta 2.6
+roughness 0.85
+sheen 0
+
+Sphere
+id 186
+position 0 -16 14
+material 186
+radius 1
+medium false
+
+Material[Dielectric]
+id 187
+albedo 1 1 1
+eta 2.6
+roughness 0.9
+sheen 0
+
+Sphere
+id 187
+position 0 -16 16
+material 187
+radius 1
+medium false
+
+Material[Dielectric]
+id 188
+albedo 1 1 1
+eta 2.6
+roughness 0.95
+sheen 0
+
+Sphere
+id 188
+position 0 -16 18
+material 188
+radius 1
+medium false
+
+Material[Dielectric]
+id 189
+albedo 1 1 1
+eta 2.6
+roughness 1
+sheen 0
+
+Sphere
+id 189
+position 0 -16 20
+material 189
+radius 1
+medium false
+
+Material[Dielectric]
+id 190
+albedo 1 1 1
+eta 2.8
+roughness 0.0001
+sheen 0
+
+Sphere
+id 190
+position 0 -18 -20
+material 190
+radius 1
+medium false
+
+Material[Dielectric]
+id 191
+albedo 1 1 1
+eta 2.8
+roughness 0.05
+sheen 0
+
+Sphere
+id 191
+position 0 -18 -18
+material 191
+radius 1
+medium false
+
+Material[Dielectric]
+id 192
+albedo 1 1 1
+eta 2.8
+roughness 0.1
+sheen 0
+
+Sphere
+id 192
+position 0 -18 -16
+material 192
+radius 1
+medium false
+
+Material[Dielectric]
+id 193
+albedo 1 1 1
+eta 2.8
+roughness 0.15
+sheen 0
+
+Sphere
+id 193
+position 0 -18 -14
+material 193
+radius 1
+medium false
+
+Material[Dielectric]
+id 194
+albedo 1 1 1
+eta 2.8
+roughness 0.2
+sheen 0
+
+Sphere
+id 194
+position 0 -18 -12
+material 194
+radius 1
+medium false
+
+Material[Dielectric]
+id 195
+albedo 1 1 1
+eta 2.8
+roughness 0.25
+sheen 0
+
+Sphere
+id 195
+position 0 -18 -10
+material 195
+radius 1
+medium false
+
+Material[Dielectric]
+id 196
+albedo 1 1 1
+eta 2.8
+roughness 0.3
+sheen 0
+
+Sphere
+id 196
+position 0 -18 -8
+material 196
+radius 1
+medium false
+
+Material[Dielectric]
+id 197
+albedo 1 1 1
+eta 2.8
+roughness 0.35
+sheen 0
+
+Sphere
+id 197
+position 0 -18 -6
+material 197
+radius 1
+medium false
+
+Material[Dielectric]
+id 198
+albedo 1 1 1
+eta 2.8
+roughness 0.4
+sheen 0
+
+Sphere
+id 198
+position 0 -18 -4
+material 198
+radius 1
+medium false
+
+Material[Dielectric]
+id 199
+albedo 1 1 1
+eta 2.8
+roughness 0.45
+sheen 0
+
+Sphere
+id 199
+position 0 -18 -2
+material 199
+radius 1
+medium false
+
+Material[Dielectric]
+id 200
+albedo 1 1 1
+eta 2.8
+roughness 0.5
+sheen 0
+
+Sphere
+id 200
+position 0 -18 0
+material 200
+radius 1
+medium false
+
+Material[Dielectric]
+id 201
+albedo 1 1 1
+eta 2.8
+roughness 0.55
+sheen 0
+
+Sphere
+id 201
+position 0 -18 2
+material 201
+radius 1
+medium false
+
+Material[Dielectric]
+id 202
+albedo 1 1 1
+eta 2.8
+roughness 0.6
+sheen 0
+
+Sphere
+id 202
+position 0 -18 4
+material 202
+radius 1
+medium false
+
+Material[Dielectric]
+id 203
+albedo 1 1 1
+eta 2.8
+roughness 0.65
+sheen 0
+
+Sphere
+id 203
+position 0 -18 6
+material 203
+radius 1
+medium false
+
+Material[Dielectric]
+id 204
+albedo 1 1 1
+eta 2.8
+roughness 0.7
+sheen 0
+
+Sphere
+id 204
+position 0 -18 8
+material 204
+radius 1
+medium false
+
+Material[Dielectric]
+id 205
+albedo 1 1 1
+eta 2.8
+roughness 0.75
+sheen 0
+
+Sphere
+id 205
+position 0 -18 10
+material 205
+radius 1
+medium false
+
+Material[Dielectric]
+id 206
+albedo 1 1 1
+eta 2.8
+roughness 0.8
+sheen 0
+
+Sphere
+id 206
+position 0 -18 12
+material 206
+radius 1
+medium false
+
+Material[Dielectric]
+id 207
+albedo 1 1 1
+eta 2.8
+roughness 0.85
+sheen 0
+
+Sphere
+id 207
+position 0 -18 14
+material 207
+radius 1
+medium false
+
+Material[Dielectric]
+id 208
+albedo 1 1 1
+eta 2.8
+roughness 0.9
+sheen 0
+
+Sphere
+id 208
+position 0 -18 16
+material 208
+radius 1
+medium false
+
+Material[Dielectric]
+id 209
+albedo 1 1 1
+eta 2.8
+roughness 0.95
+sheen 0
+
+Sphere
+id 209
+position 0 -18 18
+material 209
+radius 1
+medium false
+
+Material[Dielectric]
+id 210
+albedo 1 1 1
+eta 2.8
+roughness 1
+sheen 0
+
+Sphere
+id 210
+position 0 -18 20
+material 210
+radius 1
+medium false
+
+Material[Dielectric]
+id 211
+albedo 1 1 1
+eta 3
+roughness 0.0001
+sheen 0
+
+Sphere
+id 211
+position 0 -20 -20
+material 211
+radius 1
+medium false
+
+Material[Dielectric]
+id 212
+albedo 1 1 1
+eta 3
+roughness 0.05
+sheen 0
+
+Sphere
+id 212
+position 0 -20 -18
+material 212
+radius 1
+medium false
+
+Material[Dielectric]
+id 213
+albedo 1 1 1
+eta 3
+roughness 0.1
+sheen 0
+
+Sphere
+id 213
+position 0 -20 -16
+material 213
+radius 1
+medium false
+
+Material[Dielectric]
+id 214
+albedo 1 1 1
+eta 3
+roughness 0.15
+sheen 0
+
+Sphere
+id 214
+position 0 -20 -14
+material 214
+radius 1
+medium false
+
+Material[Dielectric]
+id 215
+albedo 1 1 1
+eta 3
+roughness 0.2
+sheen 0
+
+Sphere
+id 215
+position 0 -20 -12
+material 215
+radius 1
+medium false
+
+Material[Dielectric]
+id 216
+albedo 1 1 1
+eta 3
+roughness 0.25
+sheen 0
+
+Sphere
+id 216
+position 0 -20 -10
+material 216
+radius 1
+medium false
+
+Material[Dielectric]
+id 217
+albedo 1 1 1
+eta 3
+roughness 0.3
+sheen 0
+
+Sphere
+id 217
+position 0 -20 -8
+material 217
+radius 1
+medium false
+
+Material[Dielectric]
+id 218
+albedo 1 1 1
+eta 3
+roughness 0.35
+sheen 0
+
+Sphere
+id 218
+position 0 -20 -6
+material 218
+radius 1
+medium false
+
+Material[Dielectric]
+id 219
+albedo 1 1 1
+eta 3
+roughness 0.4
+sheen 0
+
+Sphere
+id 219
+position 0 -20 -4
+material 219
+radius 1
+medium false
+
+Material[Dielectric]
+id 220
+albedo 1 1 1
+eta 3
+roughness 0.45
+sheen 0
+
+Sphere
+id 220
+position 0 -20 -2
+material 220
+radius 1
+medium false
+
+Material[Dielectric]
+id 221
+albedo 1 1 1
+eta 3
+roughness 0.5
+sheen 0
+
+Sphere
+id 221
+position 0 -20 0
+material 221
+radius 1
+medium false
+
+Material[Dielectric]
+id 222
+albedo 1 1 1
+eta 3
+roughness 0.55
+sheen 0
+
+Sphere
+id 222
+position 0 -20 2
+material 222
+radius 1
+medium false
+
+Material[Dielectric]
+id 223
+albedo 1 1 1
+eta 3
+roughness 0.6
+sheen 0
+
+Sphere
+id 223
+position 0 -20 4
+material 223
+radius 1
+medium false
+
+Material[Dielectric]
+id 224
+albedo 1 1 1
+eta 3
+roughness 0.65
+sheen 0
+
+Sphere
+id 224
+position 0 -20 6
+material 224
+radius 1
+medium false
+
+Material[Dielectric]
+id 225
+albedo 1 1 1
+eta 3
+roughness 0.7
+sheen 0
+
+Sphere
+id 225
+position 0 -20 8
+material 225
+radius 1
+medium false
+
+Material[Dielectric]
+id 226
+albedo 1 1 1
+eta 3
+roughness 0.75
+sheen 0
+
+Sphere
+id 226
+position 0 -20 10
+material 226
+radius 1
+medium false
+
+Material[Dielectric]
+id 227
+albedo 1 1 1
+eta 3
+roughness 0.8
+sheen 0
+
+Sphere
+id 227
+position 0 -20 12
+material 227
+radius 1
+medium false
+
+Material[Dielectric]
+id 228
+albedo 1 1 1
+eta 3
+roughness 0.85
+sheen 0
+
+Sphere
+id 228
+position 0 -20 14
+material 228
+radius 1
+medium false
+
+Material[Dielectric]
+id 229
+albedo 1 1 1
+eta 3
+roughness 0.9
+sheen 0
+
+Sphere
+id 229
+position 0 -20 16
+material 229
+radius 1
+medium false
+
+Material[Dielectric]
+id 230
+albedo 1 1 1
+eta 3
+roughness 0.95
+sheen 0
+
+Sphere
+id 230
+position 0 -20 18
+material 230
+radius 1
+medium false
+
+Material[Dielectric]
+id 231
+albedo 1 1 1
+eta 3
+roughness 1
+sheen 0
+
+Sphere
+id 231
+position 0 -20 20
+material 231
+radius 1
+medium false
\ No newline at end of file
diff --git a/tests/earth.csr b/tests/earth.csr
index f168610..f989add 100644
--- a/tests/earth.csr
+++ b/tests/earth.csr
@@ -9,12 +9,16 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10.0
 
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
 Texture[Image]
 id earth_texture
 transparency false
 path earthmap.jpg
 
-Material[Lambertian]
+Material[Diffuse]
 id earth_surface
 texture earth_texture
 
diff --git a/tests/instances.csr b/tests/instances.csr
index 4b2a1a7..93ded5f 100644
--- a/tests/instances.csr
+++ b/tests/instances.csr
@@ -10,17 +10,22 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10
 
-# Init a Material object with the Lambertian type.
-Material[Lambertian]
-id red # give it a name
-texture no # we aren't giving it a custom texture
-albedo 1.0 0.0 0.0 # RedGreenBlue code
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
+Material[Diffuse]
+id red
+texture no
+albedo 1.0 0.2 0.2
+roughness 0.0001
 
 Sphere
 id sphere1
 position 0 -1 3
 material red
 radius 1
+medium false
 
 Quad
 id quad1
@@ -36,6 +41,7 @@ a 1 0 0
 b 0 1 0
 c 0 0 1
 material red
+medium false
 
 Instance[SpherePrimitive]
 prim_id sphere1
diff --git a/tests/kylo.csr b/tests/kylo.csr
index 1536223..ab7a189 100644
--- a/tests/kylo.csr
+++ b/tests/kylo.csr
@@ -10,12 +10,16 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10.0
 
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
 Texture[Image]
 id kylo_texture
 transparency true
 path kylo.png
 
-Material[Lambertian]
+Material[Diffuse]
 id kylo_material
 texture kylo_texture
 
diff --git a/tests/quads.csr b/tests/quads.csr
index cbc011e..bebc00f 100644
--- a/tests/quads.csr
+++ b/tests/quads.csr
@@ -9,30 +9,39 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10
 
-Material[Lambertian]
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
+Material[Diffuse]
 id left_red
 texture no
 albedo 1.0 0.2 0.2
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id back_green
 texture no
 albedo 0.2 1.0 0.2 
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id right_blue
 texture no
 albedo 0.2 0.2 1.0
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id upper_orange
 texture no
 albedo 1.0 0.5 0.0
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id lower_teal
 texture no
 albedo 0.2 0.8 0.8
+roughness 0.0001
 
 Quad
 id 1
diff --git a/tests/rough_metals.csr b/tests/rough_metals.csr
new file mode 100644
index 0000000..1a4935e
--- /dev/null
+++ b/tests/rough_metals.csr
@@ -0,0 +1,137 @@
+version 0.1.5
+
+Camera
+lookfrom 20 8 0
+lookat 0 4 0
+vup 0 1 0
+vfov 50
+aspect_ratio 16/9
+aperture 0.0001
+focus_dist 10.0
+
+Sky
+top 0.4 0.4 0.4
+bottom 0.4 0.4 0.4
+
+# Scene Setup with red ground
+Material[Diffuse]
+id gray_ground
+texture no
+albedo 0.8 0.4 0.4
+roughness 0.0001
+
+Sphere
+id ground
+position 0 -5000 0
+material gray_ground
+radius 5000
+medium false
+
+# Emissive
+Material[Emissive]
+id key_light
+rgb 1.0 1.0 1.0
+strength 10
+
+Sphere
+id key_sphere
+position 30 25 -30
+material key_light
+radius 10
+medium false
+
+Material[Metal]
+id m0
+albedo 1.0 1.0 1.0
+roughness 0.0
+
+Sphere
+id sphere_0
+position 0 2 14
+material m0
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.01
+albedo 1.0 1.0 1.0
+roughness 0.01
+
+Sphere
+id sphere_1
+position 0 2 10
+material m0.01
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.05
+albedo 1.0 1.0 1.0
+roughness 0.05
+
+Sphere
+id sphere_2
+position 0 2 6
+material m0.05
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.1
+albedo 1.0 1.0 1.0
+roughness 0.1
+
+Sphere
+id sphere_3
+position 0 2 2
+material m0.1
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.15
+albedo 1.0 1.0 1.0
+roughness 0.15
+
+Sphere
+id sphere_4
+position 0 2 -2
+material m0.15
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.2
+albedo 1.0 1.0 1.0
+roughness 0.2
+
+Sphere
+id sphere_5
+position 0 2 -6
+material m0.2
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.25
+albedo 1.0 1.0 1.0
+roughness 0.25
+
+Sphere
+id sphere_6
+position 0 2 -10
+material m0.25
+radius 1.75
+medium false
+
+Material[Metal]
+id m0.3
+albedo 1.0 1.0 1.0
+roughness 0.3
+
+Sphere
+id sphere_7
+position 0 2 -14
+material m0.3
+radius 1.75
+medium false
\ No newline at end of file
diff --git a/tests/simple_light.csr b/tests/simple_light.csr
index df221c0..3758540 100644
--- a/tests/simple_light.csr
+++ b/tests/simple_light.csr
@@ -9,27 +9,35 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10.0
 
-Material[Lambertian]
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
+Material[Diffuse]
 id red
 texture no
 albedo 1.0 0.2 0.2
+roughness 0.0001
 
-Material[Lambertian]
+Material[Diffuse]
 id green
 texture no
 albedo 0.2 1.0 0.2
+roughness 0.0001
 
 Sphere
 id 1
 position 0 -1000 0
 material red
 radius 1000
+medium false
 
 Sphere
 id 2
 position 0 2 0
 material green
 radius 2
+medium false
 
 Material[Emissive]
 id light_material
@@ -41,6 +49,7 @@ id lightsphere
 position 0 7 0
 material light_material
 radius 2
+medium false
 
 Quad
 id lightquad
diff --git a/tests/two_perlin_spheres.csr b/tests/two_perlin_spheres.csr
index 5bc941b..de92dfa 100644
--- a/tests/two_perlin_spheres.csr
+++ b/tests/two_perlin_spheres.csr
@@ -9,11 +9,15 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10
 
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
 Texture[Noise]
 id noise_texture
 scale 4
 
-Material[Lambertian]
+Material[Diffuse]
 id noise_material
 texture noise_texture
 
@@ -22,9 +26,11 @@ id 1
 position 0 -1000 0
 material noise_material
 radius 1000
+medium false
 
 Sphere
 id 2
 position 0 2 0
 material noise_material
-radius 2
\ No newline at end of file
+radius 2
+medium false
\ No newline at end of file
diff --git a/tests/two_spheres.csr b/tests/two_spheres.csr
index ee185f2..ea9f242 100644
--- a/tests/two_spheres.csr
+++ b/tests/two_spheres.csr
@@ -9,13 +9,17 @@ aspect_ratio 16/9
 aperture 0.0001
 focus_dist 10.0
 
+Sky
+top 0.5 0.7 1.0
+bottom 1.0 1.0 1.0
+
 Texture[Checker]
 id checker
 scale 0.8
 c1 0.2 0.3 0.1
 c2 0.9 0.9 0.9
 
-Material[Lambertian]
+Material[Diffuse]
 id checkered_surface
 texture checker
 
@@ -24,9 +28,11 @@ id top
 position 0 -10 0
 material checkered_surface
 radius 10
+medium false
 
 Sphere
 id bottom
 position 0 10 0
 material checkered_surface
-radius 10
\ No newline at end of file
+radius 10
+medium false
\ No newline at end of file
diff --git a/tools/tests.sh b/tools/tests.sh
old mode 100644
new mode 100755
index a21404a..db52bda
--- a/tools/tests.sh
+++ b/tools/tests.sh
@@ -1,16 +1,25 @@
 #!/bin/bash
 
+# Run this from inside the tools folder, as the following paths are RELATIVE to your
+# current location
+
 # Defaults for if caitlyn is in the build folder. Create a test_outputs folder in build/.
 test_outputs_dir="../build/test_outputs"
 tests_dir="../tests"
 executable="../build/caitlyn"
 
-if [ $1 == 'yes' ]
-then
+if [ "$1" == 'yes' ]; then
+    # Navigate to the build directory
+    pushd ../build > /dev/null
+
     make clean
     make
+    
+    # Return to the original directory
+    popd > /dev/null
 fi
 
 for file in "$tests_dir"/*.csr; do
-    "$executable" -i "$file" -t png -o "$test_outputs_dir/$(basename "$file" .csr).png" -s 2 -d 2 -V
+    echo "$file"
+    "$executable" -i "$file" -t png -o "$test_outputs_dir/$(basename "$file" .csr).png" -s 5 -d 200 -V
 done
\ No newline at end of file