View Raw SPL
/*****************************************************************************
*                                                                            *
*   VECDOT.SPL   Copyright (C) 2012-2022 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Computes the vector dot product and optional angle         *
*                                                                            *
*   Revisions:    21 Jun 2012  RRR  Creation                                 *
*                                                                            *
*****************************************************************************/


#if @HELP_VECDOT

    VECDOT

    Purpose: Computes the vector dot product and optional angle.

    Syntax:  VECDOT(a, b, dim)

             (vdot, vang) = VECDOT(a, b, dim)

                  a - A series, XY or XYZ series.

                  b - A series, XY or XYZ series.

                dim - Optional, the scalar direction, defaults
                      to 0, the dimension is determined by the
                      input series.

    Returns: A scalar or series, the dot product of each vector.

             (vdot, vang) = VECDOT(a, b, dim) returns the dot product and
             angle in radians between each vector.

    Example:
             u = {2, -1,  3};
             v = {5,  3, -1};

             d = vecdot(u, v);

             d == 4, the dot product of vector u and v.

    Example:
             u = {2, -1,  3};
             v = {5,  3, -1};

             (d, a) = vecdot(u, v);

             d == 11, the dot product of vector u and v and a == 1.389097
             radians (79.589372 degrees).

    Example:
             u = {2, -1,  3};
             v = {3,  3, -1};

             (d, a) = vecdot(u, v);

             d == 0, the dot product of vector u and v and a == 1.570796
             radians (90 degrees) showing that u and v are perpendicular.

    Example:
             W1: {{2, -1, 3},
                  {3,  3, 2},
                  {1,  2, 4}}

             W2: {{3,  1, 3},
                  {2,  3, 1},
                  {1, -2, 2}}

             W3: vecdot(w1, w2)

             W3 == {14, 17, 5}, the dot product of W1 and W2 where each
             row is considered a vector with X, Y, Z coordinates.

    Example:
             W1: {{2, -1, 3},
                  {3,  3, 2},
                  {1,  2, 4}}

             W2: {{3,  1, 3},
                  {2,  3, 1},
                  {1, -2, 2}}

             W3: vecdot(w1, w2, 1)

             W3 == {{13, 4, 19}}, the dot product of W1 and W2 where each
             column is considered a vector with X, Y, Z coordinates.

    Example:
             W1: {{2, -1, 3},
                  {3,  3, 2},
                  {1,  2, 4}}

             W2: vecdot(w1, w1)

             W3: vecdot(w1, w1, 1)

             W4: rowsum(w1 * w1)

             W5: colsum(w1 * w1)


             W2 == W4 == {14, 22, 21}
             W3 == W5 == {{14, 14, 29}}

             W2 and W4 compute the magnitude squared of the vectors
             row-wise.

             W3 and W5 compute the magnitude squared of the vectors
             colunm-wise.

    Example:
             W1: {{2, -1, 3},
                  {3,  3, 2},
                  {1,  2, 4}}

             W2: {{3,  1, 3},
                  {2,  3, 1},
                  {1, -2, 2}}

             W3: (d, a) = vecdot(w1, w2);ravel(d, a)

             W3 ==  {{14.0, 0.539},
                     {17.0, 0.251},
                     { 5.0, 1.199}}

             where the first column of W1 is the dot product and the second
             column is the angle.

    Example:
             u = xyz({{2, -1, 3}, {3,  3, 2}, {1,  2, 4}});
             v = xyz({{3, 1, 3},  {2, 3, 1},  {1, -2, 2}});
             
             (d, a) = vecdot(u, v);

             d == {14, 17, 5}
             a == {0.539, 0.251, 1.199}

             Same as above except the vectors are in XYZ form and the result
             is returned in two separate variables.

    Remarks:
             The dot product of two vectors A and B can be computed with:

             a . b = sum(a1 * b1 + a2 * b2 + ... + aN * bN)

             Note that a . a = ||a||^2, the magnitude squared of the vector.
             
             The angle between two vectors A and B is computed with:

             theta = acos(a . b / (||a|| * ||b||))

             DIM determines the direction of the coordinates.

             DIM == 0, data dependent
             DIM == 1, column-wise
             DIM == 2, row-wise

             For DIM == 0, if the number of columns is 3, the coordinates
             are row-wise (i, j, k), otherwise the coordinates are column-wise.

             VECDOT handles XY and XYZ series as vector inputs.

             See VECCROSS to compute the vector cross product.

    See Also:
             Norm
             Sum
             Veccross
#endif


/* dot product of two vectors */
vecdot(a, b, dim)
{
        local vdot, vangle, nr, nc;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("vecdot - input series required");

                        /* build axis projection vector */
                        (nr, nc) = size(a);
                        b = ones(nr, nc);

                        if (nc > 1)
                        {
                                b[.., nc] = zeros(nr, 1);
                        }
                }

                /* use inputs for dimension */
                dim = 0;
        }

        /* handle XY and XYZ series */
        if (isxyz(a) || isxy(a))
        {
                a = ravel(xvals(a), yvals(a), zvals(a));
        }

        if (isxyz(b) || isxy(b))
        {
                b = ravel(xvals(b), yvals(b), zvals(b));
        }

        /* dimension */
        if (dim == 0)
        {
                dim = (numcols(a) > 1 || numcols(b) > 1) ? 2 : 1;
        }
                
        /* dot product */
        vdot = vecdot_dot(a, b, dim);

        if (outargc > 1)
        {
                /* divide by product of vector lengths */
                vangle = vdot / (vecdot_mag(a, dim) * vecdot_mag(b, dim));

                /* angle in radians */
                vangle = acos(vangle);

                return(vdot, vangle);
        }
        else
        {
                return(vdot);
        }
}


/* dot product */
vecdot_dot(a, b, dim)
{
        local dot;

        dot = vecdot_sum(conj(a) * b, dim);

        return(dot);
}


/* length of vector */
vecdot_mag(a, dim)
{
        local len;

        len = sqrt(vecdot_sum(conj(a) * a, dim));

        return(len);
}


/* sum vector along dimension */
vecdot_sum(a, dim)
{
        local vsum;

        if (dim > 2)
        {
                vsum = a;
        }
        else if (dim > 1)
        {
                vsum = rowsum(a);
        }
        else
        {
                if (numcols(a) > 1)
                {
                        vsum = colsum(a);
                }
                else
                {
                        vsum = sum(a);
                }
        }

        return(vsum);
}