diff --git a/stack/maxima/contrib/vectorcalculus.mac b/stack/maxima/contrib/vectorcalculus.mac
index 9570f3b6e5215800ccee3da107913c59c05610ff..be84021519e6348e047206d2e2ec735483b9131e 100644
--- a/stack/maxima/contrib/vectorcalculus.mac
+++ b/stack/maxima/contrib/vectorcalculus.mac
@@ -1,6 +1,6 @@
/* Author Luke Longworth
University of Canterbury
- Copyright (C) 2023 Luke Longworth
+ Copyright (C) 2024 Luke Longworth
This program is free software: you can redistribute it or modify
it under the terms of the GNU General Public License version two.
@@ -16,16 +16,43 @@
/****************************************************************/
/* Vector calculus functions for STACK */
/* */
-/* V1.0 June 2023 */
+/* V2.0 March 2024 */
/* */
/****************************************************************/
+/* A flag used throughout the file. */
+/* If return_vect is true, then vector answers are returned as an nx1 matrix. */
+/* If return_vect is false, then vector answers are returned as a list. */
+return_vect: true;
+
+/****************************************************************/
+/* Calculate the gradient vector of a multivariate function */
+/****************************************************************/
+grad(f, [vars]):= block([grad_vec],
+ vars: flatten(vars),
+ if emptyp(vars) then vars: listofvars(f),
+ /* TODO: confirm grad should always simplify? */
+ grad_vec: map(lambda([ex], ev(diff(f,vars[ex]), simp)), ev(makelist(ii,ii,1,length(vars)), simp)),
+ if return_vect then return(transpose(matrix(grad_vec))) else return(grad_vec)
+);
+
+s_test_case((return_vect:true, grad(x*y*z,[x,y,z])),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:true, grad(x*y*z,x,y,z)),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:true, grad(x*y*z)),matrix([y*z],[x*z],[x*y]));
+s_test_case((return_vect:false, grad(x*y*z,[x,y,z])),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x*y*z,x,y,z)),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x*y*z)),[y*z,x*z,x*y]);
+s_test_case((return_vect:false, grad(x^2 + x)),[2*x+1]);
+s_test_case((return_vect:true, grad(a+2*b+3*c+4*d+5*p)),matrix([1],[2],[3],[4],[5]));
+s_test_case((return_vect:true, grad(a+2*b+3*c+4*d+5*p,[p,d,c,b,a])),matrix([5],[4],[3],[2],[1]));
+
/****************************************************************/
/* Calculate the divergence of a vector-valued function */
/****************************************************************/
-div(u, vars):= block([div_vec],
- /* TODO: error trapping: if not(listp(vars)) or emptyp(vars) then error("div: the second argument must be a list of variables."), */
+div(u, [vars]):= block([div_vec],
if matrixp(u) then funcs: list_matrix_entries(u) else funcs: flatten(u),
+ vars: flatten(vars),
+ if emptyp(vars) then vars: listofvars(u),
/* TODO: confirm div should always simplify? */
div_vec: map(lambda([ex], ev(diff(funcs[ex],vars[ex]), simp)), ev(makelist(ii,ii,1,length(vars)), simp)),
return(apply("+", div_vec))
@@ -34,15 +61,67 @@ div(u, vars):= block([div_vec],
s_test_case(div([x^2*cos(y),y^3],[x,y]), 2*x*cos(y)+3*y^2);
s_test_case(div(transpose(matrix([x^2*cos(y),y^3])),[x,y]), 2*x*cos(y)+3*y^2);
s_test_case(div(matrix([x^2*cos(y),y^3]),[x,y]), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3],[y,x]), -x^2*sin(y));
+s_test_case(div([y^3,x^2*cos(y)],[y,x]), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3]), 2*x*cos(y)+3*y^2);
+s_test_case(div(transpose(matrix([x^2*cos(y),y^3]))), 2*x*cos(y)+3*y^2);
+s_test_case(div(matrix([x^2*cos(y),y^3])), 2*x*cos(y)+3*y^2);
+s_test_case(div([x^2*cos(y),y^3],x,y), 2*x*cos(y)+3*y^2);
+s_test_case(div(transpose(matrix([x^2*cos(y),y^3])),x,y), 2*x*cos(y)+3*y^2);
+s_test_case(div(matrix([x^2*cos(y),y^3]),x,y), 2*x*cos(y)+3*y^2);
/****************************************************************/
/* Calculate the curl of a vector-valued function */
/****************************************************************/
-curl(u,vars):= block([cux, cuy, cuz],
- /* TODO: error trapping: if not(listp(vars)) or emptyp(vars) then error("curl: the second argument must be a list of 3 variables."), */
+curl(u, [vars]):= block([cux, cuy, cuz],
if matrixp(u) then [ux,uy,uz]: list_matrix_entries(u) else [ux,uy,uz]: flatten(u),
+ vars: flatten(vars),
+ if emptyp(vars) then vars: listofvars(u),
cux: diff(uz,vars[2]) - diff(uy,vars[3]),
cuy: diff(ux,vars[3]) - diff(uz,vars[1]),
cuz: diff(uy,vars[1]) - diff(ux,vars[2]),
- return(transpose(matrix([cux,cuy,cuz])))
+ if return_vect then return(transpose(matrix([cux,cuy,cuz]))) else return([cux,cuy,cuz])
);
+
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z],[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: false, curl([x*y*z,x*y*z,x*y*z],[x,y,z])),[x*z-x*y,x*y-y*z,y*z-x*z]);
+s_test_case((return_vect: false, curl([x*y*z,x*y*z,x*y*z])),[x*z-x*y,x*y-y*z,y*z-x*z]);
+s_test_case((return_vect: true, curl([x*y*z,x*y*z,x*y*z],[y,z,x])),matrix([x*y-y*z],[y*z-x*z],[x*z-x*y]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]),[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]),x,y,z)),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z,x*y*z,x*y*z]))),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]),[x,y,z])),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]),x,y,z)),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+s_test_case((return_vect: true, curl(matrix([x*y*z],[x*y*z],[x*y*z]))),matrix([x*z-x*y],[x*y-y*z],[y*z-x*z]));
+
+/*******************************************************************/
+/* Calculate the directional derivative of a multivariate function */
+/*******************************************************************/
+dir_deriv(f, u, [vars]):= block([unit_u, der],
+ if matrixp(u) then u: list_matrix_entries(u),
+ vars: flatten(vars),
+ if emptyp(vars) then vars: listofvars(f),
+ unit_u: u/sqrt(u . u),
+ der: ev(flatten(args(grad(f, vars))) . unit_u,simp),
+ return(der)
+);
+
+s_test_case((return_vect: false, dir_deriv(x*y*z,[1,2,2],[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,[1,2,2])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,[1,2,2],x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,matrix([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,matrix([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,matrix([1,2,2]),x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,transpose([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),[x,y,z])),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: false, dir_deriv(x*y*z,transpose([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]))),(y*z)/3+(2*x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),[y,z,x])),(2*y*z)/3+(x*z)/3+(2*x*y)/3);
+s_test_case((return_vect: true, dir_deriv(x*y*z,transpose([1,2,2]),x,y,z)),(y*z)/3+(2*x*z)/3+(2*x*y)/3);