Cálculo Vectorial en Scilab 1

Uso de Scilab para el Cálculo Vectorial

Si buscas en Google las palabras "vectores" y "Scilab", o "cálculo vectorial" y "Scilab"; te aparecen diferentes artículos donde se explica el concepto de los arrays y de los vectores en Scilab. Brevemente: arrays es una matrices de datos, agrupación de datos numéricos, o de otro tipo, agrupados en forma de matriz; y los vectores serían arrays de nx1: vector columna (sólo una columna), o 1xn vector fila (sólo una fila). Voy a utilizar la pablara array en vez de vector fila o vector columna, para no confundir con los Vectores (magnitud física).

Si lo que buscabas es el implementar el Vector (magnitud física con dirección y sentido) con Scilab, tendrá que buscar con más detalle. En este post voy a tratar de este tipo de vectores; y de como he utilizado Scilab para resolver problemas con vectores. No voy a realizar una introducción a Scilab (hay muchas y muy buenas) ni al Cálculo Vectorial (lo mismo se podría decir); lo que quiero hacer es ir directamente a presentar como he utilizado yo Scilab para resolver problemas de Cálculo Vectorial (creo que alguna de las funciones que he desarrollado pueden ser útiles para otros).

Un Vector es una magnitud  física con dirección y sentido (p.e. aceleración, velocidad, fuerza), mientras que un Escalar es una magnitud física , que sólo tiene magnitud (p.e. temperatura, presión). De manera que un Escalar lo podemos representar con un simple número, y para un Vector vamos necesitar tres números para cada una de sus coordenadas. Normalmente estamos acostumbrados a trabajar con la coordenadas cartesianas.

Entonces el Vector tendría tres coordenadas, una por caja eje: X,Y,Z.

De igual manera, si queremos representar un Punto determinado en el espacio, también necesitaremos dar sus tres coordenadas.

Primera decisión como hacemos esto en Scilab; lo más sencillo es utilizar un vector (array 1x3) para representar estos conceptos; pero tiene el inconveniente que no diferenciamos ente vector y punto. Una idea, muy recomendable, es utilizar una nomenclatura especial: por ejemplo disponer una letra "p" delante de las variables utilizadas para contener puntos y una "v" para las variables que contienen vectores. Aún podemos ir un poco más allá, y utilizar arrays 1x3 para almacenar los puntos y arrays 3x1 para almacenar los vectores; es como lo he hecho yo (después de descubrir alguna función de Scilab que requería el vector columna para su ejecución, y tener que reescribir las funciones que ya había creado). 

Por ejemplo definamos un punto (0,0,0) en Scilab (-->  es el indicador de la consola, donde se debe introducir el código que queremos ejecutar; las filas posteriores es lo que nos devuelve la consola de Scilab).

-->  pInicial=[0,0,0]
 pInicial  = 
   0.   0.   0.

Explicación: las "comas" (se podrían omitir dejando el espacio en blanco, pero no queda tan claro) delimitan cada elemento de las filas del array 1x3 del punto definido.

Introduzcamos otro punto 

--> pFinal=[1,1,1];

El "punto y coma"  final hace que la consola no nos devuelva el valor introducido, el valor se queda almacenado en la memoria, pero no nos lo muestra. Si queremos ver el valor, lo tenemos en la ventana "Explorador de variables"; o podemos consultarlo introducción su nombre en la consola:

--> pFinal
     pFinal  = 
       1.   1.   1.

Ahora creemos un Vector:

--> vVector=[1;1;1]
     vVector  = 
       1.
       1.
       1.

Queda claro la diferencia entre el Vector y el Punto (array 3x1 y 1x3), y nos puede evitar errores a la hora de usar las variables. 

Otra forma de crear este Vector sería mediante:

--> vVector2=(pFinal-pInicial)'
     vVector2  = 
       1.
       1.
       1.

Aquí he definido el Vector2 (que es el mismo que el anterior), mediante la resta al punto Final del punto Inicial y su trasposición a columna vertical (de eso se encarga el apóstrofe existente al final de la sentencia) .

Después descubrí que Scilab aceptaba muchos más caracteres ASCII, y los he utilizado ampliamente consiguiendo un resultado muy visual:

--> ↗VA=vVector
     ↗VA  = 
       1.
       1.
       1.

La flechita (y muchos más caracteres ASCII) se pueden utilizar como nombres de las variables y de las funciones que se utilicen en SCILAB; con lo que el código queda mucho más explicito.

Presentación legible por humanos

Sin embargo, esta forma de presentación no me acomodaba completamente a lo que quería realizar (solucionar problemas de Cálculo Vectorial y presentarlos de una forma legible); por lo que lo primero que hice fue crear unas funciones para presentar los valores de las variables en la consola, de una forma controlada.

Las funciones en Scilab se almacena en archivos de extensión .sci (el mismo archivo puede contener más de una función, yo las agrupo por uso, y las cargo y ejecutó antes de utilizar mis programas; para los programas se utiliza la extensión .sce).

El programa Scilab tiene un editor propio, en el menú "Aplicaciones" la opción: "SciNotes"; que nos permite crear estas funciones (.sci) (o programas: .sce). Tanto SciNotes como la consola de Scilab permiten configurar su apariencia, colores de fondo y de los distintos elementos que componen los programas (para facilitar la lectura no incluiré los colores en los códigos que utilizó, pues podría no ser legibles en otras pantallas; ya que tengo configurado el fondo negro y las letras en blanco y otros colores; en función del tipo).

Veamos mi primera función, para representar valores escalares:

function[]=ver_escalar(Titol,Escalar,Unid,Separador,Resolucio,Ocultar)
    //Titol: string con nombre magnitud
    //Escalar: Magnitud a presentar
    //Unid: string a presentar después de la magnitud
    //Resolucio: número de dígitos que tendrá el escalar (6 por defecto)
    //Ocultar: si el número es más pequeño que 1/(10^Ocultar) ⇒ 0
    if ~exists("Unid") then  Unid=""; end
    if ~exists("Separador") then  Separador=""; end
    if ~exists("Resolucio") then  Resolucio=6; end
    if ~exists("Ocultar") then  Ocultar=100; end
    if type(Escalar)==2 then // escalar con polinomios
        mprintf(Titol); mprintf('%s',pol2str(Escalar)+Unid); mprintf(Separador);
    else
        v=format(); format(Resolucio);
        if abs(Escalar)<1/(10^Ocultar) then Escalar=0; end 
        mprintf(Titol); mprintf('%s',string(Escalar)+Unid); mprintf(Separador);
        format(v(2));
    end
endfunction

Las dos barras laterales "//" indican comentario, y los utilizó para documentar mis funciones. Las variables de tipo "string" son una cadena de caracteres de longitud variable; la función "mprintf" (nativa de Scilab) es la que presenta los strings en pantalla. Veámosla en uso, para presentar un valor (una magnitud de 25,5 m):

--> ver_escalar(" Valor Escalar = ",25.5," m","\n")
     Valor Escalar = 25.5 m

En la llamada a la función hay campos que no he cumplimentado, pero como se puede comprobar en el código esto ya está previsto; y en el caso de que no se indiquen, el mismo código los detecta "~exists" y toma la decisión adecuada (para que no genere error y presente la salida deseada -por mi-).

En el campo "Separador" al llamar la función se ha introducido "\n", esto provoca un salto de línea.

Veamos otra función, esta para presentar los puntos en pantalla:

function[]=ver_punto(Titol,Punto,F,Unid,Separador,Resolucio)
//Titol: string con nombre magnitud
//Punto: punto en forma [x,y,z];[rho,phi,z];[rho,theta,phi]
//F: formato: Rectangular, Cilíndrico, Esférico
//Unid: string con unidades de la magnitud
//Resolucio: número de digitos que tendrán los componentes (6 por defecto)
if ~exists("Unid") then  Unid=""; end
if ~exists("Separador") then  Separador=""; end
if ~exists("Resolucio") then  Resolucio=6; end
if type(Punto)==2 then // vector con polinomios
    Comp_1s=pol2str(Punto(1)); 
    Comp_2s=pol2str(Punto(2)); 
    Comp_3s=pol2str(Punto(3));
    if (F=='R')|(F=='r')|(F=='0') then  //Rectangular
        text1= Comp_1s+", ";
        text2=" "+Comp_2s+", ";
        text3=" "+Comp_3s+" ";
        text="( "+text1+text2+text3+")"+Unid;
    elseif (F=='C')|(F=='c')|(F=='1') then  //Cilíndricas
        text1="ρ: "+Comp_1s+Unid;
        text2="φ: "+Comp_2s+"° ";
        text3="z: "+Comp_3s+Unid; 
        text="( "+text1+", "+text2+", "+text3+" )";
    elseif (F=='E')|(F=='e')|(F=='S')|(F=='s')|(F=='2') then  //Esféricas
        text1="r: "+Comp_1s+Unid;
        text2="θ: "+Comp_2s+"° ";
        text3="φ: "+Comp_3s+"°";
        text="( "+text1+", "+text2+", "+text3+" )";
    else
        mprintf(" Tipo de formato no reconocido \n"); text="";
    end
    mprintf('%s\n',text);
else  // vector con valores numéricos
    v=format(); format(Resolucio);
    // Basado en ver_Vector
    if (abs(Punto(1))+Punto(1))==0 then 
        Comp_1=abs(Punto(1))*(-1);
    else
        Comp_1=abs(Punto(1));
    end
    if (abs(Punto(2))+Punto(2))==0 then
        Comp_2=abs(Punto(2))*(-1);
    else
        Comp_2=abs(Punto(2));
    end
    if (abs(Punto(3))+Punto(3))==0 then
        Comp_3=abs(Punto(3))*(-1);
    else
        Comp_3=abs(Punto(3));
    end
    if (F=='R')|(F=='r')|(F=='0') then  //Rectangular
        text1= string(Comp_1)+", "; text2=string(Comp_2)+", "; 
        text3=string(Comp_3)+" ";
        text="( "+text1+text2+text3+")"+Unid;
    elseif (F=='C')|(F=='c')|(F=='1') then  //Cilíndricas
        text1="ρ: "+string(Comp_1)+" "+Unid; 
        text2="φ: "+string(Comp_2)+"° ";
        text3="z: "+string(Comp_3)+Unid; 
        text="( "+text1+", "+text2+", "+text3+")";
    elseif (F=='E')|(F=='e')|(F=='S')|(F=='s')|(F=='2') then  //Esféricas
        text1="r: "+string(Comp_1)+Unid;
        text2="θ: "+string(Comp_2)+"° ";
        text3="φ: "+string(Comp_3)+"°";
        text="( "+text1+", "+text2+", "+text3+" )";
    else
        mprintf(" Tipo de formato no reconocido \n"); text="";
    end
    mprintf(Titol); mprintf('%s',text); mprintf(Separador);
    format(v(2));
end
endfunction

Veámoslo en uso:

--> ver_punto(" Punto Final = ",pFinal,"R"," y "); ver_punto(" Punto Inicial = ",pInicial,"R");
     Punto Final = ( 1, 1, 1 ) y  Punto Inicial = ( 0, 0, 0 )

El "Formato" "R" es para la presentación de un punto en coordenadas rectangulares, más adelante entraré en otras coordenadas. Ver el uso del "Separador" que permite enlazar varias presentaciones, o saltar líneas como en el ejemplo anterior; y permite ajustar la salida de la forma que deseemos.

Por último veamos la función que cree para presentar vectores:

function[]=ver_vector(Titol,Vector,F,Unid,Separador,Resolucio,Ocultar)
    // Titol: string con nombre magnitud
    // Vector: vector en forma [x;y;z];[rho;phi;z];[rho;theta;phi]
    // F: formato: Rectangular, Cilíndrico, Esférico
    // Unid: string con unidades de la magnitud
    // Separador: string, por defecto espacio en blanco,  \n: salto de línea
    // Resolucio: número de dígitos que tendrán los componentes (6 por defecto)
    // Ocultar: si componente más pequeño que 1/(10^Ocultar) ⇒ 0
    if ~exists("Unid") then  Unid=""; end
    if ~exists("Separador") then  Separador=""; end
    if ~exists("Resolucio") then  Resolucio=6; end
    if ~exists("Ocultar") then  Ocultar=8; end
    v=format(); format(Resolucio);
    if (F=='R')|(F=='r')|(F=='0') then  //Rectangular
       coord1=" ↗i "; coord2=" ↗j "; coord3=" ↗k ";
    elseif (F=='C')|(F=='c')|(F=='1') then  //Cilíndricas 
       coord1=" ↗ρ "; coord2=" ↗φ "; coord3=" ↗z ";
    elseif (F=='E')|(F=='e')|(F=='S')|(F=='s')|(F=='2') then  //Esféricas
       coord1=" ↗r "; coord2=" ↗θ "; coord3=" ↗φ ";
    else
       mprintf(" Tipo de formato no reconocido \n"); text="";
    end
    if type(Vector)==2 then // vector con polinomios
        Comp_1s=pol2str(Vector(1)); Comp_2s=pol2str(Vector(2)); Comp_3s=pol2str(Vector(3));
        text1= Comp_1s+coord1;text2=" "+Comp_2s+coord2;text3=" "+Comp_3s+coord3;
        text="< "+text1+text2+text3+">  "+Unid;
    else  // vector con valores numéricos
        if (abs(Vector(1))+Vector(1))==0 then 
            Comp_1=abs(Vector(1))*(-1);
        else
            Comp_1=abs(Vector(1))
        end
        if (abs(Vector(2))+Vector(2))==0 then
            Comp_2=abs(Vector(2))*(-1);
        else
            Comp_2=abs(Vector(2))
        end
        if (abs(Vector(3))+Vector(3))==0 then
            Comp_3=abs(Vector(3))*(-1);
        else
            Comp_3=abs(Vector(3))
        end
        if abs(Comp_1)<1 abs="" comp_1="" comp_2="" comp_3="" coord1="" coord2="" coord3="" cultar="" else="" elseif="" end="" if="" omp_1="" omp_2="" omp_3="" string="" text1="" text2="" text3="" text="< " then=""> "+Unid;
    end
    mprintf(Titol); mprintf('%s',text); mprintf(Separador);
    format(v(2));
endfunction 

En uso:

--> ver_vector(" ↗V1 = ",vVector,"R"," m","\n"); ver_vector(" ↗V2 = ",↗VA,"R"," m","\n");
     ↗V1 = < +1 ↗i +1 ↗j +1 ↗k >  m
     ↗V2 = < +1 ↗i +1 ↗j +1 ↗k >  m

Para mi gusto, queda mucho más claro. Nota: ↗i, ↗j y ↗k son los vectores unitarios en los ejes X, Y y Z respectivamente.

Graficado de Vectores

Para presentar de manera gráfica vectores en Scilab (a partir de la versión 5.5.0) se dispone del comando xarrows, con la sintaxis:xarrows(nx, ny [, nz] [, arsize, color])

Donde:

nx, ny y nz son vectores o matrices del mismo tamaño, de forma que el elemento nx(i) da la coordenada X de origen y nx(i+1) la coordenada X final de la flecha a representar (idem para ny y nz).
arsize es un número que da el tamaño de la cabeza de la flecha (-1 por defecto)
color es una matriz o escalar que da el color de la flecha: >0 indica el color de la flecha, <0 color actual. Mediante un,Array  color(i), tamaño 1xi, se indica el color de cada una de las flechas.

Para representar las flechas en la escala adecuada, se hace necesario utilizar el comando plot3d, por ejemplo:

clf //Limpiar gráfico anterior
plot3d(xf,yf,zf)

Donde:

xf, yf y zf son matrices de tamaño [nf, n], donde n indica las caras que hay, cada cara queda definida por un polígono de n puntos; los ejes son dados por xf(:,i), yf(:,i), zf(:,i).

Por ejemplo en la ayuda de Scilab ponen una demostración con:

clf
plot3d([-1 1 1 1], [-1 -1 -1 1], [-1 -1 2 2])

Adaptándolo al caso de querer representar el vector ↗V1= [1;1;1] :  ↗V1 = < +1 ↗i +1 ↗j +1 ↗k >:

nx=[0;1];ny=[0;1];nz=[0;1];
clf, plot3d([0 1 1 1],[0 1 1 1],[0 1 1 1],10,60)
xarrows(nx,ny,nz,-1)

Se abrirá una ventana gráfica donde se mostrará el vector:

Dentro de la ventana, podemos escalar la imagen, girarla, y exportar la imagen resultante como gráfico (ente otros formatos: svg, png). 

El formato svg es un formato gráfico vectorial (una nueva utilización de la palabra vector), esto significa que le gráfico no esta guardado como puntos, sino como líneas y objetos escalables; por lo que al cambiarle la escala no pierde calidad. El programa Inkscape permite crear y editar este tipo de gráficos; y programas como TeX, LaTeX y LyX permiten utilizar este formato de imágenes.

El formato png no es un formato gráfico vectorial, sino un gráfico de puntos pixeles, es uno de los que se puede utilizar para  subir imágenes a Blogger (aún no he investigado si se pueden subir, de alguna forma imágenes svg a Blogger).

Si quisiéramos representar dos vectores en el mismo gráfico, por ejemplo añadir un vector ↗V3=[1;1;0] :  ↗V3 = < +1 ↗i +1 ↗j >  al ↗V1 antes indicado; deberemos proporcionar las coordenadas (punto inicial y punto final) de cada uno de los vectores, siendo necesario añadir una nueva columna en las variables (arrays matrices) nx, ny y nz

nx=[0 0;1 1];ny=[0 0;1 1];nz=[0 0;1 0];
clf, plot3d([0 1 1 1],[0 1 1 1],[0 1 1 1],10,60)
xarrows(nx,ny,nz,-1)

El asunto se va complicando; por lo cual, para facilitar el graficado de estas flechas he elaborado una pequeña función. Esta función utiliza unas variables globales (que no se definen en el interior  de la función, sino que deben definirse en el programa que las llame; y tienen el mismo nombre tanto dentro como fuera de la función): nx=[];ny=[];nz=[];colores=[];

function[auxx,auxy,auxz,auxc]=AddArrow(pF,pI,c)
   // adaptada para vectores verticales array 3x1  (v0=zeros(3,1))
   // pF: coordenadas del punto donde acaba la flecha
   // pI: coordenadas del punto donde comienza la flecha 
   // c: número entero: color de la flecha
   ax=[pI(1,1);pF(1,1)];ay=[pI(1,2);pF(1,2)];az=[pI(1,3);pF(1,3)];
   auxx=[nx, ax];auxy=[ny, ay];auxz=[nz, az];auxc=[colores, c];
endfunction 

Veamosla en uso en un programa, que representa la suma y resta de dos vectores:

//Previamente cargar funciones predefinidas
// exec('/ruta/AddArrow.sci',-1); //Línea anulada con comentario
printf(" Problema 1.7.3 Libro Matlab \n")
↗A=[ 3; 4; 5]; ver_vector(" ↗A = ",↗A,"R","","\n");
↗B=[-5; 4;-3]; ver_vector(" ↗B = ",↗B,"R","","\n");
↗C=↗A+↗B; ver_vector(" ↗A + ↗B = ",↗C,"R","","\n");
↗D=↗A-↗B; ver_vector(" ↗A - ↗B = ",↗D,"R","","\n");
//Graficar vectores
// Puntos iniciales y finales de las flechas
p0=zeros(1,3); pA= ↗A'; pB= ↗B'; pC= ↗C'; pD= ↗D';
// variables globales necesarias para la función AddArrow()
nx=[];ny=[];nz=[];colores=[];
// con cada llamada de la función AddArrow() se añade una flecha al gráfico
[nx,ny,nz,colores]= AddArrow(pA,p0,1);
[nx,ny,nz,colores]= AddArrow(pC,pA,1);
[nx,ny,nz,colores]= AddArrow(pC,p0,2);
[nx,ny,nz,colores]= AddArrow(pD,pA,3);
[nx,ny,nz,colores]= AddArrow(pD,p0,4);
// limpia y fija las coordenadas x,y,z y los ángulos de visión del gráfico
clf; plot3d([0 1 1 1],[0 1 1 1],[0 1 1 1],15,75);
// presenta las flechas de la gráfica;
xarrows(nx,ny,nz,-1, colores);

Obteniéndose en la consola el resultado de las operaciones; y abriéndose una ventana gráfica con la imagen de dichos vectores:

 Problema 1.7.3 Libro Matlab 
 ↗A = < +3 ↗i +4 ↗j +5 ↗k > 
 ↗B = < -5 ↗i +4 ↗j -3 ↗k > 
 ↗A + ↗B = < -2 ↗i +8 ↗j +2 ↗k > 
 ↗A - ↗B = < +8 ↗i +8 ↗k > 

Edito el post: He aprendido una forma de incluir gráficos svg en Blogguer. Estoy probándolo. Las subo a Drive, copio el enlace al archivo en Drive (poner que el enlace lo pueda ver todo el mundo) y la inserto como imagen desde URL (con tamaño original, muy grande o el que elija).

Operaciones Básicas

Suma y resta de Vectores

Las operaciones de suma y resta de vectores se realizan directamente con las variables que tenemos definidas en Scilab; ya en el post anterior realizamos una suma y resta de vectores (antes hay que cargar y ejecutar la función ver_vector).

↗A=[ 3; 4; 5]; ver_vector(" ↗A = ",↗A,"R","","\n");
↗B=[-5; 4;-3]; ver_vector(" ↗B = ",↗B,"R","","\n");
//disp(↗B,"↗B=",↗A,"↗A=")
↗C=↗A+↗B; ver_vector(" ↗A + ↗B = ",↗C,"R","","\n");
↗D=↗A-↗B; ver_vector(" ↗A - ↗B = ",↗D,"R","","\n");

Obteniendo:

 ↗A = < +3 ↗i +4 ↗j +5 ↗k > 
 ↗B = < -5 ↗i +4 ↗j -3 ↗k > 
 ↗A + ↗B = < -2 ↗i +8 ↗j +2 ↗k > 
 ↗A - ↗B = < +8 ↗i +8 ↗k >

Producto de vectores

El producto escalar de dos vectores, se realiza mediante la multiplicación ( operador *) de un array 1x3 por un array 3x1. Como los vectores los definimos como arrays 3x1, debo transponer el primer vector (array 3x1) a un array 1x3 ( mediante el operador apostrofe ').

-> ↗A'*↗B
 ans  =
  -14.

El producto vectorial de dos vectores, se realiza mediante la función cross

--> cross(↗A,↗B)
 ans  =
  -32.
  -16.
   32.

Otras operaciones con vectores

El modulo de un vector (magnitud, longitud de la flecha) se puede realizar con la función norm

--> ↗A
 ↗A  = 
   3.
   4.
   5.
--> norm(↗A)
 ans  =
   7.0710678

El vector unitario, ↗uA,  de un vector ↗A (yo lo estudié como cosenos directores) es el vector de módulo unidad y con la dirección del vector ↗A; se obtiene dividiendo (operador /)  el vector por su modulo.

--> ↗uA =↗A/norm(↗A)
 ↗uA  = 
   0.4242641
   0.5656854
   0.7071068

Dado que cada vez se iba complicando más la expresión, se me ocurrió crear una función vector que agrupase todas las operaciones a realizar con vectores. Esta función la fui desarrollando poco a poco, añadiendo las operaciones que necesité al resolver diversos problemas.

Dado el tamaño de la función, listado de código extenso, voy a presentar una versión reducida con las operaciones antes enumeradas; y después conforme vaya avanzando en el desarrollo iré analizando el código que se añade.

La función debe recibir un primer parámetro de tipo string con el que se indica la operación a realizar, y tres parámetros más con los que opera (no siempre son necesarios los tres parámetros, y en algunas operaciones el segundo parámetro es un punto).

El código de la función empieza con unos comentarios donde se intenta documentar todas las opciones y modos de uso que tiene.

Después viene una operación lógica de selección: if… ifelse…end , donde se determina que operación se quiere realizar y se ejecuta el código correspondiente (he marcado donde se insertarían el resto de operaciones, que iremos añadiendo en entradas posteriores). No obstante, en el archivo están todas las funciones que he generado (a día de hoy) y utilizó para resolver estos problemas.

function[d]=vector(ope,a_,b_,c_)
//Función que realiza diversas operaciones con vectores
//  a_, b_, c_ : vectores (de 3 componentes) (en cambio coordenadas b_ es punto evaluación)
//  ope : string descriptivo de la operación a realizar
// Operaciones con vectores (válida para todas las coord.)
//    Mod (mag) : modulo del vector a_ (válida para todas las coord.)
//    u_ (cosDirect) : vector unitario (o cosenos directores) del vector a_
// Operaciones con vectores en coord. rectangulares (salvo cambio coord.)
//    . (ProdEscalar / *e) : producto escalar entre los vectoras a_ y b_  (a_ . b_)
//    x (ProdVectorial / *v) : producto vectorial entre vectores a_ y b_  (a_ x b_)
// AQUÍ se incluirían más descripciones de las operaciones
  if ~exists("b_") then b_=[0;0;0]; end
//Operaciones
    if (ope=="Mod") | (ope=="mod")|(ope=="Mag") | (ope=="mag")  then
        d=norm(a_);
    elseif (ope=="CosDirect") | (ope=="cosDirect")| (ope=="u_") | (ope=="↗u") then
        temp=(norm(a_));
        temp=temp +((temp==0)*%eps);// +((temp==0)*%eps)) // evita singularidad 
        d=a_./temp;
    elseif (ope=="ProdEscalar") | (ope=="*e") | (ope==".") then
        d=a_'*b_; //cambio vertical
    elseif (ope=="ProdVectorial") | (ope=="*v") | (ope=="x") then
        d=cross(a_,b_);
// AQUÍ se incluirían los códigos adicionales
    else
        mprintf(" f vector : Operación no definida \n");//c=0;
    end
endfunction  // function[d]=vector(ope,a_,b_,c_)
//

Por ejemplo para el cálculo del ángulo entre dos vectores el código a incluir sería:

//    angulo : ángulo (grados decimales)entre los vectores a_ y b_
//    angulor : ángulo (radianes)entre los vectores a_ y b_

    elseif (ope=="angulo") | (ope=="alpha") then
        temp=(norm(a_)*norm(b_));
        temp=temp +((temp==0)*%eps);// +((temp==0)*%eps)) // evita singularidad 
        d=acosd((a_'*b_)/temp);
    elseif (ope=="angulor") | (ope=="alphar") then
        temp=(norm(a_)*norm(b_));
        temp=temp +((temp==0)*%eps);// +((temp==0)*%eps)) // evita singularidad 
        d=acos((a_'*b_)/temp);

Veamos la función en uso en un programa, se usan las funciones predefinidas: ver_escalar, ver_vector, vector y AddArrow. Siendo su listado el siguiente:

//Previamente cargar funciones predefinidas
↗A=[ 3; 4;  5]; ver_vector(" ↗A = ",↗A,"R","","\n");
↗B=[-5; 4; -3]; ver_vector(" ↗B = ",↗B,"R","","\n");
AB=vector(".",↗A,↗B); ver_escalar(" ↗A ● ↗B = ",AB,"","\n");
↗AxB=vector("x",↗A,↗B); ver_vector(" ↗A ∧ ↗B = ",↗AxB,"R","","\n");
Angulo=vector("angulo",↗A,↗B); ver_escalar(" Ángulo ↗A ∡ ↗B = ",Angulo,"º","\n");
AngRecto=vector("angulo",↗A,↗AxB)
ver_escalar(" Ángulo ↗A ∡ (↗A ∧ ↗B) = ",AngRecto,"º","\n")
//
//Gráficar vectores
p0=zeros(1,3); pA= ↗A'; pB= ↗B'; pAxB=↗AxB';
nx=[];ny=[];nz=[];colores=[];
[nx,ny,nz,colores]= AddArrow(pA,p0,1);
[nx,ny,nz,colores]= AddArrow(pB,p0,1);
[nx,ny,nz,colores]= AddArrow(pAxB,p0,1);
clf; plot3d([0 1 1 1],[0 1 1 1],[0 1 1 1],15,75)
xarrows(nx,ny,nz,-1, colores)

Y obteniendo en consola el resultado y una ventana gráfica con la imagen (producto vectorial):

 ↗A = < +3 ↗i +4 ↗j +5 ↗k > 
 ↗B = < -5 ↗i +4 ↗j -3 ↗k > 
 ↗A ● ↗B = -14
 ↗A ∧ ↗B = < -32 ↗i -16 ↗j +32 ↗k > 
 Ángulo ↗A ∡ ↗B = 106.3º
 Ángulo ↗A ∡ (↗A ∧ ↗B) = 90º