DISTANCE BETWEEN POINT AND RAY IN 3D

// ## task DISTANCE BETWEEN POINT AND RAY IN 3D ###

/* ### Description ###

There are a point and a ray in 3D space. Find the distance between them.

Input:
Nine reals, the 3D coordinates of:

1. The given point.

2. The initial point of the ray.

3. The direction vector of the ray.

The length of the direction vector is greater than 1e-8.

Output:
A real, the distance between the given point and the ray.

*/

#include <iostream>

#include <iomanip>

#include <cmath>

#include <vector>

using namespace std;

//vector input
template <typename T> void read_vector(vector<T>& vec)
{
    T in;

    for (int i = 0; i < 3; i++)
    {
        cin >> in;

        vec.push_back(in);
    }

}

//difference of two vectors
template <typename T> vector<T> vector_diff(const vector<T>& vec_a, const vector<T>& vec_b)
{
    vector<T> result;

    for (int i = 0; i < 3; i++)
    {
        result.push_back(vec_a[i] - vec_b[i]);
    }

    return result;
}

//sum of two vectors
template <typename T> vector<T> vector_sum(const vector<T>& vec_a, const vector<T>& vec_b)
{
    vector<T> result;

    for (int i = 0; i < 3; i++)
    {
        result.push_back(vec_a[i] + vec_b[i]);
    }

    return result;
}

//mult of number and vector
template <typename T> vector<T> vector_mult(const T a, const vector<T>& vec_b)
{
    vector<T> result;

    for (int i = 0; i < 3; i++)
    {
        result.push_back(a * vec_b[i]);
    }

    return result;
}

//scalar product of two vectors
template <typename T> T scalar_product(const vector<T>& vec_a, const vector<T>& vec_b)
{
    T res = 0.0;

    for (int i = 0; i < 3; i++)
    {
        res += vec_a[i] * vec_b[i];
    }

    return res;
}

//module or length of vector
template <typename T> double vector_module(const vector<T>& vec)
{
    double square_sum = 0.0;

    for (int i = 0; i < 3; i++)
    {
        square_sum += pow(vec[i], 2);
    }

    double res = sqrt(square_sum);

    return res;
}

int main()
{
    //point, point on vector and vector
    vector<double> p, vec_p, vec;

    //reading point
    read_vector<double>(p);

    //reading point on vector
    read_vector<double>(vec_p);

    //reading vector
    read_vector<double>(vec);

    //calculations
    double t0 = scalar_product<double>(vec, vector_diff<double>(p, vec_p)) / scalar_product<double>(vec, vec);

    double d;

    if (t0 <= 0)
    {
        d = vector_module<double>(vector_diff<double>(p, vec_p));
    }
    else
    {
        d = vector_module<double>(vector_diff<double>(p, vector_sum<double>(vec_p, vector_mult<double>(t0, vec))));
    }

    cout << fixed << setprecision(9) << d;

    return 0;
}