Implemented rectangular area lights.

Also added Cornell Box test scene.

.gitignore

@@ -3,3 +3,4 @@ target
 
 .zedstate
 test_renders
+perf.data*

example_scenes/cornell_box.psy

@@ -0,0 +1,114 @@
+Scene $Scene_fr1 {
+    Output {
+        Path ["test_renders/cornell_box.ppm"]
+    }
+    RenderSettings {
+        Resolution [512 512]
+        SamplesPerPixel [16]
+        Seed [1]
+    }
+    Camera {
+        Fov [39.449188]
+        FocalDistance [10.620000]
+        ApertureRadius [0.000000]
+        Transform [1.000000 -0.000000 0.000000 0.000000 -0.000000 0.000000 1.000000 0.000000 0.000000 1.000000 -0.000000 0.000000 -2.779998 -8.000000 2.730010 1.000000]
+    }
+    World {
+        BackgroundShader {
+            Type [Color]
+            Color [0.000000 0.000000 0.000000]
+        }
+    }
+    Assembly {
+        SurfaceShader $Green {
+            Type [Lambert]
+            Color [0.117000 0.412500 0.115000]
+        }
+        SurfaceShader $Red {
+            Type [Lambert]
+            Color [0.611000 0.055500 0.062000]
+        }
+        SurfaceShader $White {
+            Type [Lambert]
+            Color [0.729500 0.735500 0.729000]
+        }
+        RectangleLight $__Area {
+            Color [84.300003 53.800003 18.500000]
+            Dimensions [1.350000 1.100000]
+        }
+        Instance {
+            Data [$__Area]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 2.779475 -2.794788 -5.498045 1.000000]
+        }
+        MeshSurface $__Plane.010_ {
+            Vertices [-2.649998 2.959996 3.299997 -4.229996 2.469997 3.299997 -3.139998 4.559995 3.299997 -4.719996 4.059995 3.299997 -4.719996 4.059996 0.000000 -3.139998 4.559995 0.000000 -4.229996 2.469997 0.000000 -2.649998 2.959997 0.000000 ]
+            FaceVertCounts [4 4 4 4 4 ]
+            FaceVertIndices [0 1 3 2 1 0 7 6 3 1 6 4 2 3 4 5 0 2 5 7 ]
+        }
+        Instance {
+            Data [$__Plane.010_]
+            SurfaceShaderBind [$White]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane.008_ {
+            Vertices [-1.299999 0.649999 1.649998 -0.820000 2.249998 1.649999 -2.899997 1.139998 1.649999 -2.399998 2.719997 1.649999 -1.299999 0.649999 0.000000 -0.820000 2.249998 0.000000 -2.899997 1.139998 0.000000 -2.399998 2.719997 0.000000 ]
+            FaceVertCounts [4 4 4 4 4 ]
+            FaceVertIndices [0 2 3 1 3 2 6 7 1 3 7 5 0 1 5 4 2 0 4 6 ]
+        }
+        Instance {
+            Data [$__Plane.008_]
+            SurfaceShaderBind [$White]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane.006_ {
+            Vertices [-5.495996 5.591994 0.000000 -5.527995 -0.000001 -0.000000 -5.559996 5.591993 5.487995 -5.559995 -0.000001 5.487995 ]
+            FaceVertCounts [4 ]
+            FaceVertIndices [0 1 3 2 ]
+        }
+        Instance {
+            Data [$__Plane.006_]
+            SurfaceShaderBind [$Red]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane.004_ {
+            Vertices [-0.000001 5.591995 0.000000 0.000000 0.000000 0.000000 -0.000001 5.591994 5.487995 0.000000 -0.000000 5.487995 ]
+            FaceVertCounts [4 ]
+            FaceVertIndices [1 0 2 3 ]
+        }
+        Instance {
+            Data [$__Plane.004_]
+            SurfaceShaderBind [$Green]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane.002_ {
+            Vertices [-5.495996 5.591994 0.000000 -0.000001 5.591995 0.000000 -5.559996 5.591993 5.487995 -0.000001 5.591994 5.487995 ]
+            FaceVertCounts [4 ]
+            FaceVertIndices [0 1 3 2 ]
+        }
+        Instance {
+            Data [$__Plane.002_]
+            SurfaceShaderBind [$White]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane.001_ {
+            Vertices [-5.559996 5.591993 5.487995 -0.000001 5.591994 5.487995 -5.559995 -0.000001 5.487995 0.000000 -0.000000 5.487995 -3.429997 3.319996 5.487995 -2.129998 3.319996 5.487995 -3.429997 2.269997 5.487995 -2.129998 2.269997 5.487995 ]
+            FaceVertCounts [4 4 4 4 ]
+            FaceVertIndices [1 5 4 0 0 4 6 2 2 6 7 3 7 5 1 3 ]
+        }
+        Instance {
+            Data [$__Plane.001_]
+            SurfaceShaderBind [$White]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+        MeshSurface $__Plane_ {
+            Vertices [-5.495996 5.591994 0.000000 -0.000001 5.591995 0.000000 -5.527995 -0.000001 -0.000000 0.000000 0.000000 0.000000 ]
+            FaceVertCounts [4 ]
+            FaceVertIndices [0 1 3 2 ]
+        }
+        Instance {
+            Data [$__Plane_]
+            SurfaceShaderBind [$White]
+            Transform [1.000000 -0.000000 0.000000 -0.000000 -0.000000 1.000000 -0.000000 0.000000 0.000000 -0.000000 1.000000 -0.000000 -0.000000 0.000000 -0.000000 1.000000]
+        }
+    }
+}

src/lerp.rs

@@ -67,6 +67,12 @@ impl Lerp for f64 {
     }
 }
 
+impl<T: Lerp> Lerp for (T, T) {
+    fn lerp(self, other: (T, T), alpha: f32) -> (T, T) {
+        (self.0.lerp(other.0, alpha), self.1.lerp(other.1, alpha))
+    }
+}
+
 
 #[cfg(test)]
 mod tests {

src/light/mod.rs

@@ -1,8 +1,10 @@
 mod sphere_light;
+mod rectangle_light;
 
 use std::fmt::Debug;
 
 pub use self::sphere_light::SphereLight;
+pub use self::rectangle_light::RectangleLight;
 
 use math::{Vector, Point, Matrix4x4};
 use color::SpectralSample;

src/light/rectangle_light.rs

@@ -0,0 +1,157 @@
+use math::{Vector, Point, Matrix4x4, cross};
+use bbox::BBox;
+use boundable::Boundable;
+use color::{XYZ, SpectralSample, Color};
+use super::LightSource;
+use lerp::lerp_slice;
+use sampling::{spherical_triangle_solid_angle, uniform_sample_spherical_triangle};
+use std::f64::consts::PI as PI_64;
+use std::f32::consts::PI as PI_32;
+
+#[derive(Debug)]
+pub struct RectangleLight {
+    dimensions: Vec<(f32, f32)>,
+    colors: Vec<XYZ>,
+    bounds_: Vec<BBox>,
+}
+
+impl RectangleLight {
+    pub fn new(dimensions: Vec<(f32, f32)>, colors: Vec<XYZ>) -> RectangleLight {
+        let bbs = dimensions.iter()
+            .map(|d| {
+                BBox {
+                    min: Point::new(d.0 * -0.5, d.1 * -0.5, 0.0),
+                    max: Point::new(d.0 * 0.5, d.1 * 0.5, 0.0),
+                }
+            })
+            .collect();
+        RectangleLight {
+            dimensions: dimensions,
+            colors: colors,
+            bounds_: bbs,
+        }
+    }
+}
+
+impl LightSource for RectangleLight {
+    fn sample(&self,
+              space: &Matrix4x4,
+              arr: Point,
+              u: f32,
+              v: f32,
+              wavelength: f32,
+              time: f32)
+              -> (SpectralSample, Vector, f32) {
+        // Calculate time interpolated values
+        let dim = lerp_slice(&self.dimensions, time);
+        let col = lerp_slice(&self.colors, time);
+
+        // TODO: Is this right?  Do we need to get the surface area post-transform?
+        let surface_area_inv: f64 = 1.0 / (dim.0 as f64 * dim.1 as f64);
+
+        // Get the four corners of the rectangle, transformed into world space
+        let space_inv = space.inverse();
+        let p1 = Point::new(dim.0 * 0.5, dim.1 * 0.5, 0.0) * space_inv;
+        let p2 = Point::new(dim.0 * -0.5, dim.1 * 0.5, 0.0) * space_inv;
+        let p3 = Point::new(dim.0 * -0.5, dim.1 * -0.5, 0.0) * space_inv;
+        let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv;
+
+        // Get the four corners of the rectangle, projected on to the unit
+        // sphere centered around arr.
+        let sp1 = (p1 - arr).normalized();
+        let sp2 = (p2 - arr).normalized();
+        let sp3 = (p3 - arr).normalized();
+        let sp4 = (p4 - arr).normalized();
+
+        // Get the solid angles of the rectangle split into two triangles
+        let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3);
+        let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3);
+
+        // Normalize the solid angles for selection purposes
+        let prob_1 = area_1 / (area_1 + area_2);
+        let prob_2 = 1.0 - prob_1;
+
+        // Select one of the triangles and sample it
+        let shadow_vec = if u < prob_1 {
+            uniform_sample_spherical_triangle(sp2, sp1, sp3, v, u / prob_1)
+        } else {
+            uniform_sample_spherical_triangle(sp4, sp1, sp3, v, 1.0 - ((u - prob_1) / prob_2))
+        };
+
+        // Project shadow_vec back onto the light's surface
+        let arr_local = arr * *space;
+        let shadow_vec_local = shadow_vec * *space;
+        let shadow_vec_local = shadow_vec_local * (-arr_local[2] / shadow_vec_local[2]);
+        let mut sample_point_local = arr_local + shadow_vec_local;
+        sample_point_local[0] = sample_point_local[0].max(dim.0 * -0.5).min(dim.0 * 0.5);
+        sample_point_local[1] = sample_point_local[1].max(dim.1 * -0.5).min(dim.1 * 0.5);
+        sample_point_local[2] = 0.0;
+        let sample_point = sample_point_local * space_inv;
+        let shadow_vec = sample_point - arr;
+
+        // Calculate pdf and light energy
+        let pdf = 1.0 / (area_1 + area_2); // PDF of the ray direction being sampled
+        let spectral_sample = (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength);
+
+        return (spectral_sample, shadow_vec, pdf as f32);
+    }
+
+    fn sample_pdf(&self,
+                  space: &Matrix4x4,
+                  arr: Point,
+                  sample_dir: Vector,
+                  sample_u: f32,
+                  sample_v: f32,
+                  wavelength: f32,
+                  time: f32)
+                  -> f32 {
+        let dim = lerp_slice(&self.dimensions, time);
+
+        // Get the four corners of the rectangle, transformed into world space
+        let space_inv = space.inverse();
+        let p1 = Point::new(dim.0 * 0.5, dim.1 * 0.5, 0.0) * space_inv;
+        let p2 = Point::new(dim.0 * -0.5, dim.1 * 0.5, 0.0) * space_inv;
+        let p3 = Point::new(dim.0 * -0.5, dim.1 * -0.5, 0.0) * space_inv;
+        let p4 = Point::new(dim.0 * 0.5, dim.1 * -0.5, 0.0) * space_inv;
+
+        // Get the four corners of the rectangle, projected on to the unit
+        // sphere centered around arr.
+        let sp1 = (p1 - arr).normalized();
+        let sp2 = (p2 - arr).normalized();
+        let sp3 = (p3 - arr).normalized();
+        let sp4 = (p4 - arr).normalized();
+
+        // Get the solid angles of the rectangle split into two triangles
+        let area_1 = spherical_triangle_solid_angle(sp2, sp1, sp3);
+        let area_2 = spherical_triangle_solid_angle(sp4, sp1, sp3);
+
+        1.0 / (area_1 + area_2)
+    }
+
+    fn outgoing(&self,
+                space: &Matrix4x4,
+                dir: Vector,
+                u: f32,
+                v: f32,
+                wavelength: f32,
+                time: f32)
+                -> SpectralSample {
+        let dim = lerp_slice(&self.dimensions, time);
+        let col = lerp_slice(&self.colors, time);
+
+        // TODO: Is this right?  Do we need to get the surface area post-transform?
+        let surface_area_inv: f64 = 1.0 / (dim.0 as f64 * dim.1 as f64);
+
+        (col * surface_area_inv as f32 * 0.5).to_spectral_sample(wavelength)
+    }
+
+    fn is_delta(&self) -> bool {
+        false
+    }
+}
+
+impl Boundable for RectangleLight {
+    fn bounds<'a>(&'a self) -> &'a [BBox] {
+        &self.bounds_
+    }
+}

src/math/normal.rs

@@ -1,6 +1,6 @@
 #![allow(dead_code)]
 
-use std::ops::{Index, IndexMut, Add, Sub, Mul, Div};
+use std::ops::{Index, IndexMut, Add, Sub, Mul, Div, Neg};
 use std::cmp::PartialEq;
 
 use lerp::Lerp;
@@ -114,6 +114,15 @@ impl Div<f32> for Normal {
 }
 
 
+impl Neg for Normal {
+    type Output = Normal;
+
+    fn neg(self) -> Normal {
+        Normal { co: self.co * -1.0 }
+    }
+}
+
+
 impl Lerp for Normal {
     fn lerp(self, other: Normal, alpha: f32) -> Normal {
         (self * (1.0 - alpha)) + (other * alpha)

src/math/point.rs

@@ -39,6 +39,10 @@ impl Point {
 
         Point { co: n1.co.v_max(n2.co) }
     }
+
+    pub fn into_vector(self) -> Vector {
+        Vector::new(self[0], self[1], self[2])
+    }
 }
 
 

src/math/vector.rs

@@ -1,6 +1,6 @@
 #![allow(dead_code)]
 
-use std::ops::{Index, IndexMut, Add, Sub, Mul, Div};
+use std::ops::{Index, IndexMut, Add, Sub, Mul, Div, Neg};
 use std::cmp::PartialEq;
 
 use lerp::Lerp;
@@ -114,6 +114,15 @@ impl Div<f32> for Vector {
 }
 
 
+impl Neg for Vector {
+    type Output = Vector;
+
+    fn neg(self) -> Vector {
+        Vector { co: self.co * -1.0 }
+    }
+}
+
+
 impl Lerp for Vector {
     fn lerp(self, other: Vector, alpha: f32) -> Vector {
         (self * (1.0 - alpha)) + (other * alpha)

src/parse/psy_assembly.rs

@@ -5,7 +5,7 @@ use std::result::Result;
 use super::DataTree;
 use super::psy::{parse_matrix, PsyParseError};
 use super::psy_mesh_surface::parse_mesh_surface;
-use super::psy_light::parse_sphere_light;
+use super::psy_light::{parse_sphere_light, parse_rectangle_light};
 
 use assembly::{Assembly, AssemblyBuilder, Object};
 
@@ -82,6 +82,19 @@ pub fn parse_assembly(tree: &DataTree) -> Result<Assembly, PsyParseError> {
                     }
                 }
 
+                // Sphere Light
+                "RectangleLight" => {
+                    if let &DataTree::Internal { ident: Some(ident), .. } = child {
+                        builder.add_object(ident,
+                                        Object::Light(Box::new(
+                                            try!(parse_rectangle_light(&child))
+                                        )));
+                    } else {
+                        // TODO: error condition of some kind, because no ident
+                        panic!();
+                    }
+                }
+
                 _ => {
                     // TODO: some kind of error, because not a known type name
                 }

src/parse/psy_light.rs

@@ -8,7 +8,7 @@ use super::DataTree;
 use super::basics::{ws_usize, ws_f32};
 use super::psy::PsyParseError;
 
-use light::SphereLight;
+use light::{SphereLight, RectangleLight};
 use math::Point;
 use color::XYZ;
 
@@ -54,3 +54,47 @@ pub fn parse_sphere_light(tree: &DataTree) -> Result<SphereLight, PsyParseError>
         return Err(PsyParseError::UnknownError);
     }
 }
+
+pub fn parse_rectangle_light(tree: &DataTree) -> Result<RectangleLight, PsyParseError> {
+    if let &DataTree::Internal { ref children, .. } = tree {
+        let mut dimensions = Vec::new();
+        let mut colors = Vec::new();
+
+        // Parse
+        for child in children.iter() {
+            match child {
+                // Dimensions
+                &DataTree::Leaf { type_name, contents } if type_name == "Dimensions" => {
+                    if let IResult::Done(_, radius) =
+                           closure!(tuple!(ws_f32, ws_f32))(contents.as_bytes()) {
+                        dimensions.push(radius);
+                    } else {
+                        // Found dimensions, but its contents is not in the right format
+                        return Err(PsyParseError::UnknownError);
+                    }
+                }
+
+                // Color
+                &DataTree::Leaf { type_name, contents } if type_name == "Color" => {
+                    if let IResult::Done(_, color) = closure!(tuple!(ws_f32,
+                                                                     ws_f32,
+                                                                     ws_f32))(contents.as_bytes()) {
+                        // TODO: handle color space conversions properly.
+                        // Probably will need a special color type with its
+                        // own parser...?
+                        colors.push(XYZ::new(color.0, color.1, color.2));
+                    } else {
+                        // Found color, but its contents is not in the right format
+                        return Err(PsyParseError::UnknownError);
+                    }
+                }
+
+                _ => {}
+            }
+        }
+
+        return Ok(RectangleLight::new(dimensions, colors));
+    } else {
+        return Err(PsyParseError::UnknownError);
+    }
+}

src/renderer.rs

@@ -260,8 +260,12 @@ impl LightPath {
                         let (_, shadow_vec, _) =
                             light.sample(&space, pos, lu, lv, self.wavelength, self.time);
 
-                        let la = dot(nor.normalized().into_vector(), shadow_vec.normalized())
+                        let rnor = if dot(nor.into_vector(), ray.dir) > 0.0 {
-                            .max(0.0);
+                            -nor.into_vector().normalized()
+                        } else {
+                            nor.into_vector().normalized()
+                        };
+                        let la = dot(rnor, shadow_vec.normalized()).max(0.0);
                         self.light_attenuation = XYZ::from_spectral_sample(&XYZ::new(la, la, la)
                             .to_spectral_sample(self.wavelength));
                         *ray = Ray::new(pos + shadow_vec.normalized() * 0.0001,

src/sampling.rs

@@ -1,4 +1,4 @@
-use math::Vector;
+use math::{Vector, dot};
 
 use std::f64::consts::PI as PI_64;
 use std::f32::consts::PI as PI_32;
@@ -69,3 +69,96 @@ pub fn uniform_sample_cone_pdf(cos_theta_max: f64) -> f64 {
     // 1.0 / solid angle
     1.0 / (2.0 * PI_64 * (1.0 - cos_theta_max))
 }
+
+/// Calculates the projected solid angle of a spherical triangle.
+///
+/// A, B, and C are the points of the triangle on a unit sphere.
+pub fn spherical_triangle_solid_angle(va: Vector, vb: Vector, vc: Vector) -> f32 {
+    // Calculate sines and cosines of the spherical triangle's edge lengths
+    let cos_a: f64 = dot(vb, vc).max(-1.0).min(1.0) as f64;
+    let cos_b: f64 = dot(vc, va).max(-1.0).min(1.0) as f64;
+    let cos_c: f64 = dot(va, vb).max(-1.0).min(1.0) as f64;
+    let sin_a: f64 = (1.0 - (cos_a * cos_a)).sqrt();
+    let sin_b: f64 = (1.0 - (cos_b * cos_b)).sqrt();
+    let sin_c: f64 = (1.0 - (cos_c * cos_c)).sqrt();
+
+    // If two of the vertices are coincident, area is zero.
+    // Return early to avoid a divide by zero below.
+    if cos_a == 1.0 || cos_b == 1.0 || cos_c == 1.0 {
+        return 0.0;
+    }
+
+    // Calculate the cosine of the angles at the vertices
+    let cos_va = ((cos_a - (cos_b * cos_c)) / (sin_b * sin_c)).max(-1.0).min(1.0);
+    let cos_vb = ((cos_b - (cos_c * cos_a)) / (sin_c * sin_a)).max(-1.0).min(1.0);
+    let cos_vc = ((cos_c - (cos_a * cos_b)) / (sin_a * sin_b)).max(-1.0).min(1.0);
+
+    // Calculate the angles themselves, in radians
+    let ang_va = cos_va.acos();
+    let ang_vb = cos_vb.acos();
+    let ang_vc = cos_vc.acos();
+
+    // Calculate and return the solid angle of the triangle
+    (ang_va + ang_vb + ang_vc - PI_64) as f32
+}
+
+/// Generates a uniform sample on a spherical triangle given two uniform
+/// random variables i and j in [0, 1].
+pub fn uniform_sample_spherical_triangle(va: Vector,
+                                         vb: Vector,
+                                         vc: Vector,
+                                         i: f32,
+                                         j: f32)
+                                         -> Vector {
+    // Calculate sines and cosines of the spherical triangle's edge lengths
+    let cos_a: f64 = dot(vb, vc).max(-1.0).min(1.0) as f64;
+    let cos_b: f64 = dot(vc, va).max(-1.0).min(1.0) as f64;
+    let cos_c: f64 = dot(va, vb).max(-1.0).min(1.0) as f64;
+    let sin_a: f64 = (1.0 - (cos_a * cos_a)).sqrt();
+    let sin_b: f64 = (1.0 - (cos_b * cos_b)).sqrt();
+    let sin_c: f64 = (1.0 - (cos_c * cos_c)).sqrt();
+
+    // If two of the vertices are coincident, area is zero.
+    // Return early to avoid a divide by zero below.
+    if cos_a == 1.0 || cos_b == 1.0 || cos_c == 1.0 {
+        // TODO: do something more intelligent here, in the case that it's
+        // an infinitely thin line.
+        return va;
+    }
+
+    // Calculate the cosine of the angles at the vertices
+    let cos_va = ((cos_a - (cos_b * cos_c)) / (sin_b * sin_c)).max(-1.0).min(1.0);
+    let cos_vb = ((cos_b - (cos_c * cos_a)) / (sin_c * sin_a)).max(-1.0).min(1.0);
+    let cos_vc = ((cos_c - (cos_a * cos_b)) / (sin_a * sin_b)).max(-1.0).min(1.0);
+
+    // Calculate sine for A
+    let sin_va = (1.0 - (cos_va * cos_va)).sqrt();
+
+    // Calculate the angles themselves, in radians
+    let ang_va = cos_va.acos();
+    let ang_vb = cos_vb.acos();
+    let ang_vc = cos_vc.acos();
+
+    // Calculate the area of the spherical triangle
+    let area = ang_va + ang_vb + ang_vc - PI_64;
+
+    // The rest of this is from the paper "Stratified Sampling of Spherical
+    // Triangles" by James Arvo.
+    let area_2 = area * i as f64;
+
+    let s = (area_2 - ang_va).sin();
+    let t = (area_2 - ang_va).cos();
+    let u = t - cos_va;
+    let v = s + (sin_va * cos_c);
+
+    let q_top = (((v * t) - (u * s)) * cos_va) - v;
+    let q_bottom = ((v * s) + (u * t)) * sin_va;
+    let q = q_top / q_bottom;
+
+    let vc_2 = (va * q as f32) +
+               ((vc - (va * dot(vc, va))).normalized() * (1.0 - (q * q)).sqrt() as f32);
+
+    let z = 1.0 - (j * (1.0 - dot(vc_2, vb)));
+
+    (vb * z) + ((vc_2 - (vb * dot(vc_2, vb))).normalized() * (1.0 - (z * z)).sqrt())
+}