Vector.h

#ifndef VECTOR_H
#define VECTOR_H

#include <array>
#include <cmath>
#include <ostream>

typedef struct Vector Vector;
double dot(const Vector& lhs, const Vector& rhs);

struct Vector{
  double x, y, z;
  Vector() : x{0}, y{0}, z{0} {}
  Vector(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} {}
  Vector(const std::array<double, 3>& xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} {}

  double length() const{
    return sqrt(x*x + y*y + z*z);
  }

  Vector normalize(){
    double magnitude = length();
    x /= magnitude;
    y /= magnitude;
    z /= magnitude;
    return *this;
  }

  Vector operator-(const Vector& rhs) const{
    return Vector(x - rhs.x, y - rhs.y, z - rhs.z);
  }

  Vector operator-() const{
    return Vector(-x, -y, -z);
  }

  Vector operator+(const Vector& rhs) const{
    return Vector(x + rhs.x, y + rhs.y, z + rhs.z);
  }

  Vector operator*(const double& rhs) const{
    return Vector(x * rhs, y * rhs, z * rhs);
  }

  //assumign both vectors point away from the reflection point
  Vector reflect(const Vector& normal) const{
    return normal * 2 * dot(*this, normal) - *this;
  }

  // --------- Taken from scratchapixel.com ---------
  // https://www.scratchapixel.com/lessons/3d-basic-rendering/introduction-to-shading/reflection-refraction-fresnel.html
  Vector refract(const Vector &normal, const float &indexOfRefraction) const{
    float cosi = dot(*this, normal);
    float etai = 1, etat = indexOfRefraction;
    Vector n = normal;
    if (cosi < 0) { //comming from outside
      cosi = -cosi; 
    } 
    else { //comming from the inside
      std::swap(etai, etat);
      n = -normal; 
    }
    float eta = etai / etat;
    float k = 1 - eta * eta * (1 - cosi * cosi);
    if (k < 0) { //total internal reflection
      return (-*this).reflect(normal);
    }
    else {
      return *this * eta + n * (eta * cosi - sqrtf(k));
    }
  }
};

inline double dot(const Vector& lhs, const Vector& rhs) {
  return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
}

inline Vector cross(const Vector& lhs, const Vector& rhs) {
  return Vector(
            lhs.y * rhs.z - lhs.z * rhs.y,
            lhs.z * rhs.x - lhs.x * rhs.z,
            lhs.x * rhs.y - lhs.y * rhs.x
  );
}

inline Vector getNormalized(const Vector& v){
  return Vector(v.x, v.y, v.z).normalize();
}

inline std::ostream& operator<< (std::ostream& os, const Vector& v){
  os << "(" << v.x << ", " << v.y << ", " << v.z << ")";
  return os;
}

#endif