math.go

/*
Most of the code here will be taken from HandmadeMath, 
which is covered under the CC0 1.0 Universal license,
see LICENSE for more information
*/
package math

import omath "core:math"

PI :: omath.PI;
PI2 :: PI * 2;

// Use this to switch between f32 and f64
Float :: f32;

Vec2       :: distinct [2]Float;
Vec3       :: distinct [3]Float;
Vec4       :: distinct [4]Float;
Quaternion :: distinct [4]Float;
Mat4       :: distinct [4][4]Float;

toradians :: proc(angle: Float) -> Float {
	return angle / 180 * PI;
}

lerp :: proc(a: Float, time: Float, b: Float) -> Float {
	return (1.0 - time) * a + time * b;
}

clamp :: proc(min: Float, value: Float, max: Float) -> Float {
	if value < min do value = min;
	if value > max do value = max;
	return value;
}

dot :: proc [
	dot_vec2,
	dot_vec3,
	dot_vec4,
	quaternion_dot,
];

dot_vec2 :: proc(a, b: Vec2) -> Float {
	return (a.y * b.y) + (a.y * b.y);
}

dot_vec3 :: proc(a, b: Vec3) -> Float {
	return (a.x * b.x) + (a.y * b.y) + (a.z * b.z);
}

dot_vec4 :: proc(a, b: Vec4) -> Float {
	return (a.x * b.x) + (a.y * b.y) + (a.z * b.z) + (a.w * b.w);	
}

cross :: proc(a, b: Vec3) -> (result: Vec3) {
	result.x = (a.y * b.z) - (a.z * b.y);
	result.y = (a.z * b.x) - (a.x * b.z);
	result.z = (a.x * b.y) - (a.y * b.x);
	
	return;
}

length_squared :: proc[
	length_squared_vec2,
	length_squared_vec3,
	length_squared_vec4,
];

length_squared_vec2 :: proc(a: Vec2) -> Float {
	return a.x*a.x + a.y*a.y;
}

length_squared_vec3 :: proc(a: Vec3) -> Float {
	return a.x*a.x + a.y*a.y + a.z*a.z;
}

length_squared_vec4 :: proc(a: Vec4) -> Float {
	return a.x*a.x + a.y*a.y + a.z*a.z + a.w*a.w;
}

length :: proc[
	length_vec2,
	length_vec3,
	length_vec4,
];

length_vec2 :: proc(a: Vec2) -> Float {
	return omath.sqrt(length_squared_vec2(a));
}

length_vec3 :: proc(a: Vec3) -> Float {
	return omath.sqrt(length_squared_vec3(a));
}

length_vec4 :: proc(a: Vec4) -> Float {
	return omath.sqrt(length_squared_vec4(a));
}

normalize :: proc[
	normalize_vec2,
	normalize_vec3,
	normalize_vec4,
	quaternion_normalize,
];

normalize_vec2 :: proc(a: Vec2) -> (result: Vec2) {
	length := length(a);
	if length != 0 {
		result.x = a.x / length;
		result.y = a.y / length;
	}
	return;
}

normalize_vec3 :: proc(a: Vec3) -> (result: Vec3) {
	length := length(a);
	if length != 0 {
		result.x = a.x / length;
		result.y = a.y / length;
		result.z = a.z / length;
	}
	return;
}

normalize_vec4 :: proc(a: Vec4) -> (result: Vec4) {
	length := length(a);
	if length != 0 {
		result.x = a.x / length;
		result.y = a.y / length;
		result.z = a.z / length;
		result.w = a.w / length;
	}
	return;
}

mat4 :: proc(diagonal: Float) -> (result: Mat4) {
	result[0][0] = diagonal;
	result[1][1] = diagonal;
	result[2][2] = diagonal;
	result[3][3] = diagonal;
	return;
}

orthographic :: proc(left, right, bottom, top: Float, near, far: Float) -> (result: Mat4) {
	result[0][0] = 2.0 / (right - left);
	result[1][1] = 2.0 / (top - bottom);
	result[2][2] = 2.0 / (near - far);
	result[3][3] = 1.0;
	result[3][0] = (left + right) / (left - right);
	result[3][1] = (bottom + top) / (bottom - top);
	result[3][2] = (far + near) / (near - far);
	return;
}

perspective :: proc(fov, aspect, near, far: Float) -> (result: Mat4) {
	tanthetaover2 := omath.tan(fov * (PI / 360.0));
	
	result[0][0] = 1.0 / tanthetaover2;
	result[1][1] = aspect / tanthetaover2;
	result[2][3] = -1.0;
	result[2][2] = (near + far) / (near - far);
	result[3][2] = (2.0 * near * far) / (near - far);
	result[3][3] = 0.0;
	
	return;
}

translate :: proc(translation: Vec3) -> (result: Mat4) {
	result[3][0] = translation.x;
	result[3][1] = translation.y;
	result[3][2] = translation.z;
	
	return;
}

scale :: proc(scale: Vec3) -> (result: Mat4) {
	result[0][0] = scale.x;
	result[1][1] = scale.y;
	result[2][2] = scale.z;
	
	return;
}

transpose :: proc(matrix: Mat4) -> (result: Mat4) {
	for col := 0; col < 4; col += 1 {
		for row := 0; row < 4; row += 1 {
			result[row][col] = matrix[col][row];
		}
	}
	return;
}

rotate :: proc(angle: Float, axis: Vec3) -> (result: Mat4) {
	axis = normalize(axis);
	
	sintheta: Float = omath.sin(toradians(angle));
	costheta: Float = omath.cos(toradians(angle));
	cosvalue: Float = 1.0 - costheta;
	
	result[0][0] = (axis.x * axis.x * cosvalue) + costheta;
	result[0][1] = (axis.x * axis.y * cosvalue) + (axis.z * sintheta);
	result[0][2] = (axis.x * axis.z * cosvalue) - (axis.y * sintheta);
	result[1][0] = (axis.y * axis.x * cosvalue) - (axis.z * sintheta);
	result[1][1] = (axis.y * axis.y * cosvalue) + costheta;
	result[1][2] = (axis.y * axis.z * cosvalue) + (axis.x * sintheta);
	result[2][0] = (axis.z * axis.x * cosvalue) + (axis.y * sintheta);
	result[2][1] = (axis.z * axis.y * cosvalue) - (axis.x * sintheta);
	result[2][2] = (axis.z * axis.z * cosvalue) + costheta;
	
	return;
}

look_at :: proc(eye, center, up: Vec3) -> (result: Mat4) {
	f := normalize(center - eye);
	s := normalize(cross(f, up));
	u := cross(s, f);

	result[0][0] = s.x;
	result[0][1] = u.x;
	result[0][2] = -f.x;
	result[0][3] = 0.0;
	result[1][0] = s.y;
	result[1][1] = u.y;
	result[1][2] = -f.y;
	result[1][3] = 0.0;
	result[2][0] = s.z;
	result[2][1] = u.z;
	result[2][2] = -f.z;
	result[2][3] = 0.0;
	result[3][0] = -dot(s, eye);
	result[3][1] = -dot(u, eye);
	result[3][2] = dot(f, eye);
	result[3][3] = 1.0;
	
	return;
}

mul_mat4_vec4 :: proc(matrix: Mat4, vec: Vec4) -> (result: Vec4) {
	for row := 0; row < 4; row += 1 {
		sum: Float = 0.0;
		for col := 0; col < 4; col += 1 {
			sum += matrix[col][row] * vec[col];
		}
		
		result[row] = sum;
	}
	return;
}

quaternion_mul :: proc(left, right: Quaternion) -> (result: Quaternion) {
	result.x = (left.x * right.w) + (left.y * right.z) - (left.z * right.y) + (left.w * right.x);
	result.y = (-left.x * right.z) + (left.y * right.w) + (left.z * right.x) + (left.w * right.y);
	result.z = (left.x * right.y) - (left.y * right.x) + (left.z * right.w) + (left.w * right.z);
	result.w = (-left.x * right.x) - (left.y * right.y) - (left.z * right.z) + (left.w * right.w);
	
	return;
}

quaternion_dot :: proc(left, right: Quaternion) -> Float {
	return (left.x * right.x) + (left.y * right.y) + (left.z * right.z) + (left.w * right.w);
}

quaternion_normalize :: proc(left: Quaternion) -> (result: Quaternion) {
	length := omath.sqrt(quaternion_dot(left, left));
	result = left / length;	
	return;
}

quaternion_inverse :: proc(left: Quaternion) -> (result: Quaternion) {
	conjugate := Quaternion{-left.x, -left.y, -left.z, left.w};
	norm: Float = omath.sqrt(dot(left, left));
	normsquared := norm * norm;

	result.x = conjugate.x / normsquared;
	result.y = conjugate.y / normsquared;
	result.z = conjugate.z / normsquared;
	result.w = conjugate.w / normsquared;
	
	return;
}

nlerp :: proc(left: Quaternion, time: Float, right: Quaternion) -> (result: Quaternion) {
	result.x = lerp(left.x, time, right.x);
	result.y = lerp(left.y, time, right.y);
	result.z = lerp(left.z, time, right.z);
	result.w = lerp(left.w, time, right.w);

	result = normalize(result);
	
	return;
}

slerp :: proc(left: Quaternion, time: Float, right: Quaternion) -> (result: Quaternion) {
	assert(false, "We dont have acos atm. so this is not implemented");
	cos_theta := dot(left, right);
	angle: Float = 0;/*math.acos(cos_theta);*/
	
	s1: Float = omath.sin((1.0 - time) * angle);
	s2: Float = omath.sin(time * angle);
	is: Float = 1.0 / omath.sin(angle);

	quaternionleft := left * s1;
	quaternionright := right * s2;

	result = quaternionleft + quaternionright;
	result = quaternion_mul(result, is);
	
	return;
}

//TODO:
// HMM_QuaternionToMat4
// HMM_QuaternionFromAxisAngle