PL/pgSQL vector and matrix functions

vector.sql

CREATE OR REPLACE FUNCTION get_vector_from_angles(azimuth double precision, altitude double precision) RETURNS double precision[] AS
$$
DECLARE
    vector double precision[];
BEGIN
    vector := ARRAY [ sin(azimuth) * cos(altitude),cos(azimuth) * cos(altitude), sin(altitude) ];

    RETURN vector;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION named_vector(direction text) RETURNS double precision[] AS
$$
BEGIN
    CASE direction
        WHEN 'up' THEN RETURN ARRAY [0, 0, 1];
        WHEN 'down' THEN RETURN ARRAY [0, 0, -1];
        WHEN 'left' THEN RETURN ARRAY [-1, 0, 0];
        WHEN 'right' THEN RETURN ARRAY [1, 0, 0];
        WHEN 'forward' THEN RETURN ARRAY [0, 1, 0];
        WHEN 'backward' THEN RETURN ARRAY [0, -1, 0];
        ELSE RAISE EXCEPTION 'Invalid direction parameter.'; END CASE;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION rotate(axis character, theta float) RETURNS double precision[][] AS
$$
DECLARE
    matrix double precision[][];
BEGIN
    CASE axis
        WHEN 'x' THEN matrix :=
                ARRAY [ ARRAY [1, 0, 0], ARRAY [0, cos(theta), -sin(theta)], ARRAY [0, sin(theta), cos(theta)] ];
        WHEN 'y' THEN matrix :=
                ARRAY [ ARRAY [cos(theta), 0, sin(theta)], ARRAY [0, 1, 0], ARRAY [-sin(theta), 0, cos(theta)] ];
        WHEN 'z' THEN matrix :=
                ARRAY [ ARRAY [cos(theta), -sin(theta), 0], ARRAY [sin(theta), cos(theta), 0], ARRAY [0, 0, 1] ];
        ELSE RAISE EXCEPTION 'Invalid axis parameter. Please specify "x", "y", or "z".'; END CASE;

    RETURN matrix;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION matrix_vector_product(vector double precision[], matrix double precision[][]) RETURNS double precision[] AS
$$
DECLARE
    result      double precision[];
    num_rows    integer := array_length(matrix, 2);
    num_columns integer := array_length(vector, 1);
    i           integer;
    j           integer;
BEGIN
    IF array_length(vector, 1) <> num_columns THEN
        RAISE EXCEPTION 'Matrix and vector dimensions are not compatible.';
    END IF;

    result := array_fill(0.0, ARRAY [num_rows]);

    FOR i IN 1..num_rows
        LOOP
            FOR j IN 1..num_columns
                LOOP
                    result[i] := result[i] + matrix[i][j] * vector[j];
                END LOOP;
        END LOOP;

    RETURN result;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION vector_multiply(vector1 double precision[], vector2 double precision[]) RETURNS double precision[] AS
$$
DECLARE
    result double precision[];
    i      integer;
BEGIN
    IF array_length(vector1, 1) <> array_length(vector2, 1) THEN
        RAISE EXCEPTION 'Vector dimensions are not compatible for multiplication.';
    END IF;

    result := array_fill(0.0, ARRAY [array_length(vector1, 1)]);

    FOR i IN 1..array_length(vector1, 1)
        LOOP
            result[i] := vector1[i] * vector2[i];
        END LOOP;

    RETURN result;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION vector_scalar_product(vector1 double precision[], vector2 double precision[]) RETURNS double precision AS
$$
DECLARE
    result double precision := 0.0;
    i      integer;
BEGIN
    IF array_length(vector1, 1) <> array_length(vector2, 1) THEN
        RAISE EXCEPTION 'Vector dimensions are not compatible for scalar product.';
    END IF;

    FOR i IN 1..array_length(vector1, 1)
        LOOP
            result := result + (vector1[i] * vector2[i]);
        END LOOP;

    RETURN result;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

CREATE OR REPLACE FUNCTION angle_of_incidence_from_az_alt(tilt double precision, rotation double precision,
                                                          azimuth double precision,
                                                          altitude double precision) RETURNS double precision AS
$$
DECLARE
    rotated_plane      double precision[];
    vector_from_angles double precision[];
    result double precision;
BEGIN
    vector_from_angles := get_vector_from_angles(pi() + azimuth, altitude);
    rotated_plane := matrix_vector_product(matrix_vector_product(named_vector('up'), rotate('x', radians(tilt))),
                                           rotate('z', radians(rotation)));
    result := degrees(acos(vector_scalar_product(vector_from_angles, rotated_plane)));
    return result;
END;
$$ LANGUAGE plpgsql  IMMUTABLE;