/*****************************************************************************
*                                                                            *
*   VECCROSS.SPL  Copyright (C) 2022 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Computes the vector cross product                          *
*                                                                            *
*   Revisions:    17 May 2022  RRR  Creation                                 *
*                                                                            *
*****************************************************************************/


#if @HELP_VECCROSS

    VECCROSS

    Purpose: Computes the vector cross product.

    Syntax:  VECCROSS(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 series, the coordinates of one or more a vectors that
             are perpendicular to the input vectors.

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

             c = veccross(u, v);

             c == {-12,
                    -1,
                    39}

             The input and output coordinates are column-wise.

             To verify that vector C is perpendicular to both U and V:

             d1 = vecdot(c, u);
             d2 = vecdot(c, v);

             d1 == d2 == 0.0.

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

             c = veccross(u, v);

             c == {{-12, -1, 39}}

             Same as above except the coordinates are row-wise.

    Example:
             W1: ravel(1..9, 3)
             W2: ravel(3..11, 3)
             W3: veccross(w1, w2)
             W4: veccross(w1, w2, 2)

             W3 == {{-6, 12, -6},
                    {-6, 12, -6},
                    {-6, 12, -6}}

             W4 == {{-6, 12, -6},
                    {-6, 12, -6},
                    {-6, 12, -6}}

            Both inputs are 3x3 arrays. Both cross products are computed with
            the coordinates going row-wise.

    Example:
             W1: ravel(1..9, 3)
             W2: ravel(3..11, 3)
             W3: veccross(w1, w2, 1)

             W3 == {{-2, -2, -2},
                    { 4,  4,  4},
                    {-2, -2, -2}}

             Same as above except the coordinates go column-wise.

    Remarks:
             The cross product of two vectors A and B is defined as:
                                            
             a x b = ||a|| * ||b|| * sin(theta) n

             where:
             
             theta is the angle between A and B in the plane that
             contains them.

             ||a|| and ||b|| are the magnitudes of A and B.

             n is the unit vector perpendicular to the plane that contains
             A and B.

             The cross product of two 3x1 vectors A and B can be computed with:

             a x b = {a[2]*b[3] - a[3] * b[2],
                      a[1]*b[3] - a[3] * b[1],                        
                      a[1]*b[2] - a[2] * b[1]}

             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.

             VECCROSS handles XY and XYZ series as vector inputs.

             See VECDOT to compute the vector cross product.

    See Also:
             Norm
             Sum
             Vecdot
#endif


/* dot product of two vectors */
veccross(a, b, dim)
{
        local vcross, nr, nc;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("veccross - 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) == 3 || numcols(b) == 3) ? 2 : 1;
        }
                
        /* dot product */
        vcross = veccross_cross(a, b, dim, __FUNC__);

        return(vcross);
}


/* cross product */
veccross_cross(a, b, dim, func)
{
        local nra, nrb, nca, ncb, c;

        (nra, nca) = size(a);
        (nrb, ncb) = size(b);

        if (not(nra == nrb && nca == ncb))
        {
                error(sprintf("%s - input dimensions must be equal", func));
        }

        if (((dim == 1) && not(nra == 3)) || ((dim == 2) && not(nca == 3)))
        {
                error(sprintf("%s - input dimension must be 3 in cross product direction", func));
        }

        if (dim == 1)
        {
                c = {a[2,..] * b[3,..] - a[3,..] * b[2,..],
                         a[3,..] * b[1,..] - a[1,..] * b[3,..],
                         a[1,..] * b[2,..] - a[2,..] * b[1,..]};

                setcomment(c, "i j k", 1, 1);
        }
        else
        {
                c = {a[..,2] * b[..,3] - a[..,3] * b[..,2],
                         a[..,3] * b[..,1] - a[..,1] * b[..,3],
                         a[..,1] * b[..,2] - a[..,2] * b[..,1]};

                /* reshape */
                c = ravel(c, length(a));

                /* attributes */
                setcomment(c, "i", 1, 1);
                setcomment(c, "j", 1, 2);
                setcomment(c, "k", 1, 3);
        }

        return(c);
}