aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2001-06-14 11:53:21 +0000
committerjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2001-06-14 11:53:21 +0000
commit6a0395d7e63372be04d20ee447b5d88af4f09f69 (patch)
tree925ce0fd4046ff4bf7b05f7f1c7672be858836fe
parent54ea67b814e9e5e8a00c797b320c0ad2b02bf4f0 (diff)
This commit was generated by cvs2svn to compensate for changes in r2,
which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@3 f88db872-0e4f-0410-b76b-b9085cfa78c5
-rw-r--r--README51
-rw-r--r--interface.ccl2
-rw-r--r--param.ccl2
-rw-r--r--schedule.ccl2
-rw-r--r--src/gr/H.c159
-rw-r--r--src/gr/Makefile29
-rw-r--r--src/gr/gr_gfas.minc25
-rw-r--r--src/gr/horizon_function.maple127
-rw-r--r--src/gr/setup_gfas.maple49
-rw-r--r--src/make.code.defn9
-rw-r--r--src/maple/Diff.maple260
-rw-r--r--src/maple/Makefile29
-rw-r--r--src/maple/coeffs.maple10
-rw-r--r--src/maple/coords.maple68
-rw-r--r--src/maple/coords.minc18
-rw-r--r--src/maple/gfa.maple315
-rw-r--r--src/maple/setup.maple24
-rw-r--r--src/maple/util.maple234
18 files changed, 1413 insertions, 0 deletions
diff --git a/README b/README
new file mode 100644
index 0000000..c5badbf
--- /dev/null
+++ b/README
@@ -0,0 +1,51 @@
+Cactus Code Thorn AHFinderDirect
+Authors : Jonathan Thornburg <jthorn@aei.mpg.de>
+CVS info : $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn
+====================
+
+This thorn finds an apparent horizon given the 3D xyz-grid metric and
+extrinsic curvature. It uses a direct method, writing the apparent
+horizon equation as an elliptic PDE on angular-coordinate space. This
+is fast, but does require an intitial guess for the apparent horizon
+position.
+
+
+Documentation
+=============
+
+The doc/ directory contains a user's guide and a programmer's guide
+for this thorn.
+
+
+Distribution
+============
+
+This thorn should eventually become part of the CactusEinstein arrangement.
+
+This thorn is distributed under terms of the GNU public license; see
+the file COPYING in this directory for details.
+
+
+Other Software Required
+=======================
+
+* This thorn requires the Einstein thorn.
+* This thorn is written in C++, so you'll need a C++ compiler to compile
+ this thorn. It calls various C and Fortran numerical libraries; if
+ these aren't already installed on your system, then you'll need C
+ and Fortran compilers to compile them.
+* Most of this thorn's relativity code is machine-generated using Maple
+ (version 6), but you don't need Maple unless you want to modify the
+ relativity code.
+
+
+Portability
+===========
+
+This thorn is known to compile and run successfully on the following
+systems:
+* Red Hat GNU/Linux, gcc 2.95.3
+* Alpha OSF1 V4.0, gcc 2.95.3
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..3905af4
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,2 @@
+# Interface definition for thorn AHFinderDirect
+# $Header$
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..9e673cc
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,2 @@
+# Parameter definitions for thorn AHFinderDirect
+# $Header$
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..815801e
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,2 @@
+# Schedule definitions for thorn AHFinderDirect
+# $Header$
diff --git a/src/gr/H.c b/src/gr/H.c
new file mode 100644
index 0000000..c345f78
--- /dev/null
+++ b/src/gr/H.c
@@ -0,0 +1,159 @@
+ t1 = g_uu[2][2];
+ t2 = 1/r;
+ t4 = X_ud[0][2];
+ t5 = Diff(h,rho);
+ t7 = X_ud[1][2];
+ t8 = Diff(h,sigma);
+ t10 = zz*t2-t4*t5-t7*t8;
+ t11 = t1*t10;
+ t12 = g_uu[1][1];
+ t14 = X_ud[0][1];
+ t16 = X_ud[1][1];
+ t18 = yy*t2-t14*t5-t16*t8;
+ t19 = t12*t18;
+ t21 = r*r;
+ t23 = 1/r/t21;
+ t30 = Diff(h,rho,rho);
+ t33 = Diff(h,rho,sigma);
+ t38 = Diff(h,sigma,sigma);
+ t40 = -yy*zz*t23-X_udd[0][1][2]*t5-X_udd[1][1][2]*t8-t14*t4*t30-t16*t4*
+t33-t14*t7*t33-t16*t7*t38;
+ t44 = g_uu[0][2];
+ t45 = t44*t44;
+ t46 = t10*t10;
+ t48 = yy*yy;
+ t49 = zz*zz;
+ t56 = X_ud[0][0];
+ t57 = t56*t56;
+ t59 = X_ud[1][0];
+ t63 = t59*t59;
+ t65 = (t48+t49)*t23-X_udd[0][0][0]*t5-X_udd[1][0][0]*t8-t57*t30-2.0*t59*
+t56*t33-t63*t38;
+ t70 = xx*t2-t56*t5-t59*t8;
+ t71 = t70*t70;
+ t72 = t70*t71;
+ t74 = Diff(g_uu[0][0],zz);
+ t77 = g_uu[1][2];
+ t78 = t77*t10;
+ t79 = xx*xx;
+ t86 = t14*t14;
+ t91 = t16*t16;
+ t93 = (t79+t49)*t23-X_udd[0][1][1]*t5-X_udd[1][1][1]*t8-t86*t30-2.0*t16*
+t14*t33-t91*t38;
+ t97 = g_uu[0][1];
+ t99 = Diff(g_uu[0][0],yy);
+ t102 = t77*t18;
+ t109 = t4*t4;
+ t114 = t7*t7;
+ t116 = (t79+t48)*t23-X_udd[0][2][2]*t5-X_udd[1][2][2]*t8-t109*t30-2.0*t7*
+t4*t33-t114*t38;
+ t120 = t44*t10;
+ t135 = -xx*zz*t23-X_udd[0][0][2]*t5-X_udd[1][0][2]*t8-t56*t4*t30-t59*t4*
+t33-t56*t7*t33-t59*t7*t38;
+ t139 = g_uu[0][0];
+ t140 = t139*t70;
+ t141 = Diff(g_uu[1][2],xx);
+ t142 = t141*t10;
+ t145 = t44*t70;
+ t146 = Diff(g_uu[1][2],zz);
+ t147 = t146*t10;
+ t151 = Diff(g_uu[0][0],xx);
+ t154 = t10*t46;
+ t156 = Diff(g_uu[2][2],zz);
+ t159 = t97*t70;
+ t160 = Diff(g_uu[1][2],yy);
+ t161 = t160*t10;
+ t164 = -2.0*t11*t19*t40-t45*t46*t65-t44*t72*t74/2.0-2.0*t78*t19*t93-t97*
+t72*t99/2.0-2.0*t11*t102*t116-2.0*t120*t102*t135-t140*t142*t18-t145*t147*t18-
+t139*t72*t151/2.0-t1*t154*t156/2.0-t159*t161*t18;
+ t166 = Diff(g_uu[2][2],yy);
+ t170 = Diff(g_uu[2][2],xx);
+ t173 = t77*t46;
+ t188 = -xx*yy*t23-X_udd[0][0][1]*t5-X_udd[1][0][1]*t8-t56*t14*t30-t59*t14
+*t33-t56*t16*t33-t59*t16*t38;
+ t192 = t18*t18;
+ t193 = t18*t192;
+ t195 = Diff(g_uu[1][1],zz);
+ t199 = Diff(g_uu[1][1],xx);
+ t202 = t195*t192;
+ t205 = t199*t192;
+ t209 = Diff(g_uu[1][1],yy);
+ t212 = t1*t1;
+ t215 = t1*t46;
+ t216 = t77*t40;
+ t219 = t97*t192;
+ t221 = t12*t192;
+ t223 = -t77*t154*t166/2.0-t44*t154*t170/2.0-2.0*t173*t44*t188-t77*t193*
+t195/2.0-t97*t193*t199/2.0-t145*t202/2.0-t140*t205/2.0-t12*t193*t209/2.0-t212*
+t46*t116-2.0*t215*t216-t219*t142-t221*t161;
+ t225 = t74*t71;
+ t228 = t99*t71;
+ t235 = t151*t71;
+ t238 = t156*t46;
+ t241 = t77*t192;
+ t243 = t209*t192;
+ t246 = t77*t77;
+ t249 = t97*t18;
+ t252 = t170*t46;
+ t255 = t166*t46;
+ t258 = -t102*t225/2.0-t19*t228/2.0-t11*t225/2.0-t78*t228/2.0-t120*t235/
+2.0-t145*t238/2.0-t241*t147-t159*t243/2.0-t246*t46*t93-t249*t235/2.0-t140*t252/
+2.0-t159*t255/2.0;
+ t263 = Diff(g_uu[0][1],zz);
+ t266 = Diff(g_uu[0][1],yy);
+ t267 = t266*t70;
+ t269 = t139*t71;
+ t270 = Diff(g_uu[0][1],xx);
+ t271 = t270*t18;
+ t273 = t97*t71;
+ t274 = Diff(g_uu[0][2],yy);
+ t275 = t274*t10;
+ t277 = Diff(g_uu[0][2],xx);
+ t278 = t277*t10;
+ t282 = t266*t18;
+ t285 = t263*t18;
+ t290 = t44*t71;
+ t291 = Diff(g_uu[0][2],zz);
+ t292 = t291*t10;
+ t297 = -t249*t252/2.0-t102*t238/2.0-t241*t263*t70-t221*t267-t269*t271-
+t273*t275-t269*t278-t19*t255/2.0-t78*t282*t70-t11*t285*t70-t120*t271*t70-t290*
+t292-2.0*t273*t139*t188;
+ t303 = t97*t97;
+ t313 = t97*t188;
+ t316 = t12*t12;
+ t332 = -2.0*t290*t97*t40-t303*t71*t93-2.0*t290*t139*t135-t45*t71*t116-
+t303*t192*t65-2.0*t221*t313-t316*t192*t93-2.0*t241*t97*t135-2.0*t241*t12*t40
+-2.0*t45*t10*t70*t135-t246*t192*t116-t273*t282;
+ t333 = t139*t139;
+ t336 = t145*t40;
+ t351 = t44*t46;
+ t357 = t291*t70;
+ t362 = t159*t40;
+ t365 = -t333*t71*t65-2.0*t78*t336-t249*t278*t70-t102*t292*t70-t19*t275*
+t70-2.0*t11*t249*t135-2.0*t78*t249*t188-t351*t277*t70-t173*t274*t70-t290*t285-
+t215*t357-2.0*t120*t19*t188-2.0*t11*t362;
+ t371 = t140*t65;
+ t374 = t140*t135;
+ t379 = t140*t188;
+ t382 = t159*t93;
+ t401 = -2.0*t246*t10*t18*t40-2.0*t120*t371-2.0*t11*t374-t219*t270*t70-2.0
+*t78*t379-2.0*t78*t382-2.0*t120*t159*t188-2.0*t19*t382-2.0*t303*t18*t70*t188
+-2.0*t19*t379-2.0*t249*t371-2.0*t249*t145*t135;
+ t406 = t145*t116;
+ t422 = t44*t135;
+ t427 = t146*t18;
+ t431 = -2.0*t102*t362-2.0*t102*t374-2.0*t11*t406-2.0*t102*t406-2.0*t19*
+t336-t11*t202/2.0-t78*t243/2.0-2.0*t120*t249*t65-t120*t205/2.0-2.0*t215*t422-
+t351*t141*t18-t215*t427-t173*t160*t18;
+ t441 = t269+2.0*t249*t70+2.0*t120*t70+t221+2.0*t78*t18+t215;
+ t442 = sqrt(t441);
+ t452 = t151*t70+t267+t357+t271+t209*t18+t427+t278+t161+t156*t10+t139*t65+
+2.0*t313+2.0*t422;
+ t456 = partial_d_ln_sqrt_g[0];
+ t459 = partial_d_ln_sqrt_g[1];
+ t462 = partial_d_ln_sqrt_g[2];
+ t477 = t12*t93+2.0*t216+t1*t116+t139*t456*t70+t97*t459*t70+t44*t462*t70+
+t97*t456*t18+t12*t459*t18+t77*t462*t18+t44*t456*t10+t77*t459*t10+t1*t462*t10;
+ H = (t164+t223+t258+t297+t332+t365+t401+t431)/t442/t441+(t452+t477)/t442+
+(K_uu[0][0]*t71+2.0*K_uu[0][1]*t18*t70+2.0*K_uu[0][2]*t10*t70+K_uu[1][1]*t192+
+2.0*K_uu[1][2]*t10*t18+K_uu[2][2]*t46)/t441-K;
diff --git a/src/gr/Makefile b/src/gr/Makefile
new file mode 100644
index 0000000..892450c
--- /dev/null
+++ b/src/gr/Makefile
@@ -0,0 +1,29 @@
+# Makefile for GR code
+# $Id: Makefile,v 1.1.1.1 2001-06-14 11:53:21 jthorn Exp $
+
+#
+# Environment Variables:
+# MAPLE_VERSION used via @ifdef for version control in Maple code;
+# typically set to something like MAPLE_V_RELEASE_4
+#
+# Targets:
+# mm ==> preprocess all *.maple files to produce *.mm files
+# clean ==> delete *.mm (mpp Maple preprocessor output)
+# run ==> load all the code into Maple
+#
+
+###############################################################################
+
+.PHONY : mm
+mm : $(patsubst %.maple, %.mm, $(wildcard *.maple))
+
+%.mm : %.maple $(wildcard *.minc $(gfa_dir)/*.minc $(gfa_dir)/*.maple)
+ mpp -D$(MAPLE_VERSION) <$< >$@
+
+.PHONY : clean
+clean :
+ -rm *.mm
+
+.PHONY : run
+run : mm
+ maple -f <setup.mm >maple.log
diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
new file mode 100644
index 0000000..92e781e
--- /dev/null
+++ b/src/gr/gr_gfas.minc
@@ -0,0 +1,25 @@
+# Maple include file listing all gfa functional-dependence forms
+
+g_dd,
+K_dd,
+
+g_uu, g_uu__fnd,
+K_uu, K_uu__fnd,
+
+K, K__fnd,
+
+partial_d_ln_sqrt_g, partial_d_ln_sqrt_g__fnd,
+
+h, h__fnd,
+
+s_d, s_d__fnd,
+partial_d_s_d, partial_d_s_d__fnd,
+
+n_u, n_u__fnd,
+
+H, H__fnd,
+
+HA, HA__fnd,
+HB, HB__fnd,
+HC, HC__fnd,
+HD, HD__fnd # no comma
diff --git a/src/gr/horizon_function.maple b/src/gr/horizon_function.maple
new file mode 100644
index 0000000..c1750ea
--- /dev/null
+++ b/src/gr/horizon_function.maple
@@ -0,0 +1,127 @@
+# horizon.maple - evaluate LHS of horizon function
+#
+# non_unit_normal - compute s_d
+# non_unit_normal_gradient - compute partial_d_s_d
+# horizon_function - compute HA,HB,HC,HD,H = fn(s_d)
+#
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function computes the non-unit outward-pointing normal vector
+# (field) to the horizon, using the equation on p.7.1 of my apparent
+# horizon finding notes.
+#
+non_unit_normal :=
+proc()
+global
+ @include "../maple/coords.minc",
+ @include "../gr/gr_gfas.minc";
+local i, u;
+
+assert_fnd_exists(h);
+assert_fnd_exists(X_ud);
+assert_fnd_exists(s_d, fnd);
+
+ for i from 1 to N
+ do
+ s_d__fnd[i] := x_xyz[i]/r
+ - sum('X_ud[u,i]*Diff(h,x_rrs[u])', 'u'=A_min..N);
+ end do;
+
+NULL;
+end proc;
+
+
+################################################################################
+
+#
+# This function computes the gradient of the non-unit outward-pointing
+# normal vector (field) to the horizon, using the equation on p.7.2 of
+# my apparent horizon finding notes.
+#
+non_unit_normal_gradient :=
+proc()
+global
+ @include "../maple/coords.minc",
+ @include "../gr/gr_gfas.minc";
+local temp, i, j, k, u, v;
+
+assert_fnd_exists(h);
+assert_fnd_exists(X_ud);
+assert_fnd_exists(X_udd);
+assert_fnd_exists(partial_d_s_d, fnd);
+
+ for i from 1 to N
+ do
+ for j from 1 to N
+ do
+ temp := `if`(i=j, sum('x_xyz[k]^2','k'=1..3), 0);
+ # simplify likes to put things over a common denominator,
+ # so we only apply it to the individual terms, not to
+ # the whole experssion
+ partial_d_s_d__fnd[i,j]
+ := + simplify( (temp-x_xyz[i]*x_xyz[j])/r^3 )
+ - simplify( sum('X_udd[u,i,j]*Diff(h,x_rrs[u])', 'u'=A_min..N) )
+ - simplify( msum('X_ud[u,i]*X_ud[v,j]*Diff(h,x_rrs[u],x_rrs[v])',
+ 'u'=A_min..N, 'v'=A_min..N) );
+ end do;
+ end do;
+
+NULL;
+end proc;
+
+################################################################################
+
+#
+# This function computes the LHS of the apparent horizon equation
+# as a function of 1st and 2nd angular partial derivatives of the
+# apparent horizon radius $h$, in the Maple gridfn array H.fnd ,
+# using equation (15) in my 1996 apparent horizon finding paper.
+#
+horizon_function :=
+proc()
+global
+ @include "../maple/coords.minc",
+ @include "../gr/gr_gfas.minc";
+local i,j,k,l;
+
+assert_fnd_exists(g_uu);
+assert_fnd_exists(K_uu);
+assert_fnd_exists(partial_d_ln_sqrt_g);
+assert_fnd_exists(s_d, fnd);
+assert_fnd_exists(partial_d_s_d, fnd);
+
+#
+# n.b. no simplify() in these subexpressions; empirically it makes
+# them much *more* complicated (measured by codegen(...,optimized))
+#
+
+# (15a) in my paper
+HA__fnd := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l]
+ * partial_d_s_d__fnd[i,j]',
+ 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N)
+ - (1/2)*msum('g_uu[i,j]*s_d__fnd[j]
+ * Diff(g_uu[k,l],x_xyz[i]) * s_d__fnd[k]*s_d__fnd[l]',
+ 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
+
+# (15b) in my paper
+HB__fnd := + msum('Diff(g_uu[i,j],x_xyz[i])*s_d__fnd[j]', 'i'=1..N, 'j'=1..N)
+ + msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N)
+ + msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
+ 'i'=1..N, 'j'=1..N);
+
+# (15c) in my paper
+HC__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+
+# (15d) in my paper
+HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+
+# n.b. no simplify() here, it would try to put things over a common
+# denominator, which would make the equation much more complicated
+H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K;
+
+NULL;
+end proc;
diff --git a/src/gr/setup_gfas.maple b/src/gr/setup_gfas.maple
new file mode 100644
index 0000000..5da990b
--- /dev/null
+++ b/src/gr/setup_gfas.maple
@@ -0,0 +1,49 @@
+# Maple code to set up all gridfn arrays
+# $Id$
+
+#
+# setup_gr_gfas - setup all the GR gfas
+#
+
+################################################################################
+
+setup_gr_gfas :=
+proc()
+global
+ @include "../gr/gr_gfas.minc";
+
+# basic metric & extrinsic curvature fields
+make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric);
+make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric);
+
+make_gfa('g_uu', {inert,fnd}, [1..N, 1..N], symmetric);
+make_gfa('K_uu', {inert,fnd}, [1..N, 1..N], symmetric);
+
+make_gfa('K', {inert, fnd}, [], none);
+
+# partial derivatives of metric determinant
+make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none);
+
+# radius of horizon
+make_gfa('h', {inert, fnd}, [], none);
+
+# outward-pointing *non*-unit normal to horizon
+# and it's xyz-coordinate partial derivatives
+make_gfa('s_d', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none);
+
+# outward-pointing unit normal to horizon
+make_gfa('n_u', {inert,fnd}, [1..N], none);
+
+# LHS of apparent horizon equation
+make_gfa('H', {inert, fnd}, [], none);
+
+# subexpressions for computing LHS of apparent horizon equation
+# ... these are defined by (15) in my 1996 apparent horizon finding paper
+make_gfa('HA', {inert, fnd}, [], none);
+make_gfa('HB', {inert, fnd}, [], none);
+make_gfa('HC', {inert, fnd}, [], none);
+make_gfa('HD', {inert, fnd}, [], none);
+
+NULL;
+end proc;
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..a3d158e
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn AHFinderDirect
+# $Header$
+
+# Source files in this directory
+SRCS =
+
+# Subdirectories containing source files
+SUBDIRS =
+
diff --git a/src/maple/Diff.maple b/src/maple/Diff.maple
new file mode 100644
index 0000000..c21d79b
--- /dev/null
+++ b/src/maple/Diff.maple
@@ -0,0 +1,260 @@
+# Diff.maple -- Diff() and friends
+# $Id$
+
+#
+# simplify/Diff - simplify(...) function for expressions containing Diff()
+# Diff - "inert" derivative function (low-level simplifications)
+# Diff/gridfn - simplifications that know about gridfns
+#
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This function is called automatically by simplify() whenever
+# simplify()'s argument contains a Diff() call. This interface
+# to simplify() is described in
+# Michael B. Monagan
+# "Tips for Maple Users"
+# The Maple Technical Newsletter
+# Issue 10 (Fall 1993), p.10-18, especially p.17-18
+# but not in the Maple library manual or the online help, which both
+# imply that simplify(..., Diff) is needed for this procedure to be
+# invoked).
+#
+# This function recurses over the expression structure. For each
+# Diff() call, it calls Diff() itself to do "basic" derivative
+# simplifications, then calls `Diff/gridfn`() to do further simplifications
+# based on gridfn properties.
+#
+# expr = (in) An expression to be simplified.
+#
+# Results:
+# The function returns the "simplified" expression.
+#
+`simplify/Diff` :=
+proc(expr)
+option remember; # performance optimization
+local temp;
+
+# recurse on tables, equations, lists, sets, sums, products, powers
+if type(expr, {table, `=`, list, set, `+`, `*`, `^`})
+ then return map(simplify, expr);
+end if;
+
+# return unchanged on numbers, names, procedures
+if type(expr, {numeric, name, procedure})
+ then return expr;
+end if;
+
+# function call ==> simplify Diff() calls, recurse on others
+if type(expr, function)
+ then
+ if (op(0, expr) = 'Diff')
+ then return Diff(op(expr)); # "basic" derivative simplifications
+ else
+ temp := map(simplify, [op(expr)]);
+ op(0,expr); return '%'(op(temp));
+ end if;
+end if;
+
+# unknown type
+ERROR(`expr has unknown type`, `whattype(expr)=`, whattype(expr));
+
+end proc;
+
+###############################################################################
+
+#
+# This function implements some of the more useful simplification
+# rules of algebra and calculus for Maple's "inert" differentiation
+# function Diff() . It then calls `Diff/gridfn`() to do any
+# further simplifications based on gridfn properties.
+#
+# In particular:
+# - It distributes over sums.
+# - It distributes over leading numeric factors. (Maple's built-in
+# simplifier will always make numeric factors leading.)
+# - It knows that derivatives of numeric constants are zero.
+# - It knows that the derivative of a name with respect to itself is one.
+# - It knows the product rule
+# Diff(f*g,x) = Diff(f,x)*g + f*Diff(g,x)
+# - It knows the power rule
+# Diff(f^n,x) = n * f^(n-1) * Diff(f,x)
+# - It flattens multi-level "Diff() trees", i.e. it makes the transformation
+# Diff(Diff(f,x),y) = Diff(f,x,y);
+# - It canonicalizes the order of the things we're differentiating with
+# respect to, so that simplify() will recognize that derivatives
+# commute, i.e. that
+# Diff(f,x,y) = Diff(f,y,x)
+#
+# Arguments:
+# operand = (in) The thing to be differentiated.
+# var_seq = (in) (varargs) An expression sequence of the variables to
+# differentiate with respect to.
+#
+unprotect('Diff');
+Diff :=
+proc(operand) # varargs
+option remember; # performance optimization
+local var_list,
+ nn, nderiv,
+ op_cdr,
+ abs_operand,
+ f, g, x, x_car, x_cdr, temp,
+ n,
+ inner_operand, inner_var_list, k,
+ sorted_var_list;
+
+var_list := [args[2..nargs]];
+
+nn := nops(operand);
+nderiv := nops(var_list);
+
+##print(`Diff(): nn=`,nn,` nderiv=`,nderiv,
+## ` operand=`,operand,` var_list=`,var_list);
+
+# distribute (recurse) over sums
+if (type(operand, `+`))
+ then return simplify(
+ sum('Diff(op(k, operand), op(var_list))', 'k'=1..nn)
+ );
+end if;
+
+# distribute (recurse) over leading numeric factors
+if (type(operand, `*`) and (nn >= 1) and type(op(1,operand), numeric))
+ then
+ op_cdr := product('op(k,operand)', 'k'=2..nn);
+ return simplify( op(1,operand) * Diff(op_cdr, op(var_list)) );
+end if;
+
+# derivatives of numbers are zero
+if (type(operand, numeric))
+ then return 0;
+end if;
+
+# derivative of a name with respect to itself is one
+if (type(operand, name) and (var_list = [operand]))
+ then return 1;
+end if;
+
+# implement the product rule
+# ... we actually implement it only for a 1st derivative of the
+# product of two factors; for fancier cases we recurse
+if (type(operand, `*`))
+ then
+ if (nn = 0)
+ then # empty product = 1
+ return 0; # ==> Diff(1) = 0
+
+ elif (nn = 1)
+ then # singleton product is equal to the single factor
+ return simplify( Diff(op(1,operand), op(var_list)) );
+
+ elif (nn >= 2)
+ then # nontrivial product rule case
+
+ # note that if we have more than 2 factors, g will
+ # be the product of the cdr of the factor list, so
+ # our later Diff(g,...) calls will recurse as
+ # necessary to handle this case
+
+ f := op(1,operand);
+ g := product('op(k,operand)', 'k'=2..nn);
+
+ if (nderiv = 1)
+ then # product rule for 1st derivatives:
+ # Diff(f*g,x) = Diff(f,x)*g + f*Diff(g,x)
+ x := var_list[1];
+
+ return simplify( Diff(f,x)*g + f*Diff(g,x) );
+
+ else # product rule for 2nd (and higher) derivatives
+ # ==> recurse on derivatives
+ x_car := var_list[1];
+ x_cdr := var_list[2..nderiv];
+ temp := simplify( Diff(f*g, x_car) );
+ return simplify( Diff(temp, op(x_cdr)) );
+ end if;
+ else ERROR(`impossible value of nn in product rule!`,
+ `(this should never happen!)`,
+ ` operand=`,operand,` nn=`,nn);
+ end if;
+end if;
+
+# implement the power rule
+# ... we actually implement it only for 1st derivatives; for higher
+# derivatives we recurse
+if (type(operand, `^`))
+ then
+ f := op(1,operand);
+ n := op(2,operand);
+
+ if (nderiv = 1)
+ then # power rule for 1st derivatives:
+ # Diff(f^n,x) = n * f^(n-1) * Diff(f,x)
+ x := var_list[1];
+
+ return simplify( n * f^(n-1) * Diff(f, x) );
+
+ else # power rule for 2nd (and higher) derivatives:
+ # ==> recurse on derivatives
+ x_car := var_list[1];
+ x_cdr := var_list[2..nderiv];
+ temp := simplify( Diff(f^n, x_car) );
+
+ return simplify( Diff(temp, op(x_cdr)) );
+ end if;
+end if;
+
+# flatten (recurse on) multi-level "Diff() trees"
+if (type(operand, function) and (op(0, operand) = 'Diff'))
+ then
+ inner_operand := op(1, operand);
+ inner_var_list := [op(2..nn, operand)];
+
+ return simplify(
+ Diff(inner_operand, op(inner_var_list), op(var_list))
+ );
+end if;
+
+# get to here
+# ==> no other transformations to make
+# ==> canonicalize ordering of variables and return reconstructed Diff() call
+sorted_var_list := sort_var_list(var_list);
+return `Diff/gridfn`(operand, op(sorted_var_list));
+
+end proc;
+
+###############################################################################
+
+#
+# This function implements further simplification rules for Diff()
+# based on gridfn properties.
+#
+# It currently knows about the following simplifications:
+# - Diff(X_ud[u,i], x_xyz[j]) --> X_udd[u,i,j]
+#
+# Arguments:
+# operand = (in) The thing to be differentiated.
+# var_seq = (in) (varargs) An expression sequence of the variables to
+# differentiate with respect to.
+#
+`Diff/gridfn` :=
+proc(operand) # varargs
+option remember; # performance optimization
+global
+ @include "../maple/coords.minc";
+local var_list, posn;
+
+var_list := [args[2..nargs]];
+
+if ( type(operand, indexed) and (op(0,operand) = 'X_ud')
+ and (nops(var_list) = 1) and member(var_list[1],x_xyz_list,'posn') )
+ then return X_udd[op(operand), posn];
+ else return 'Diff'(operand, op(var_list)); # return is *unevaluated*
+ # to avoid infinite recursion!
+end if;
+
+end proc;
diff --git a/src/maple/Makefile b/src/maple/Makefile
new file mode 100644
index 0000000..c46eea4
--- /dev/null
+++ b/src/maple/Makefile
@@ -0,0 +1,29 @@
+# Makefile for maple code
+# $Id: Makefile,v 1.1.1.1 2001-06-14 11:53:21 jthorn Exp $
+
+#
+# Environment Variables:
+# MAPLE_VERSION used via @ifdef for version control in Maple code;
+# typically set to something like MAPLE_V_RELEASE_4
+#
+# Targets:
+# mm ==> preprocess all *.maple files to produce *.mm files
+# clean ==> delete *.mm (mpp Maple preprocessor output)
+# run ==> load all the code into Maple
+#
+
+###############################################################################
+
+.PHONY : mm
+mm : $(patsubst %.maple, %.mm, $(wildcard *.maple))
+
+%.mm : %.maple $(wildcard *.minc $(gfa_dir)/*.minc $(gfa_dir)/*.maple)
+ mpp -D$(MAPLE_VERSION) <$< >$@
+
+.PHONY : clean
+clean :
+ -rm *.mm
+
+.PHONY : run
+run : mm
+ maple -f <setup.mm >maple.log
diff --git a/src/maple/coeffs.maple b/src/maple/coeffs.maple
new file mode 100644
index 0000000..b45114b
--- /dev/null
+++ b/src/maple/coeffs.maple
@@ -0,0 +1,10 @@
+# coeffs.maple - set up coordinate coefficients
+# $Id$
+#
+# setup_X_ud
+# setup_X_udd
+#
+
+################################################################################
+
+
diff --git a/src/maple/coords.maple b/src/maple/coords.maple
new file mode 100644
index 0000000..4e8953d
--- /dev/null
+++ b/src/maple/coords.maple
@@ -0,0 +1,68 @@
+# coords.maple -- set up and covert between coordinates, gridfn arrays, etc
+# $Id$
+
+#
+# setup_coords - set up coordinates
+# setup_coord_gfas - setup up coordinate-transformation gfas
+#
+
+###############################################################################
+
+#
+# This function sets up the coordinates.
+# It also sets up the Kronecker delta
+#
+setup_coords :=
+proc()
+global
+@include "coords.minc";
+
+N := 3; # number of spatial coordinates,
+ # i.e. spatial coordinate indices
+ # are 1..N
+N_ang := 2; # number of angular spatial coordinates
+A_min := 2; # min angular-spatial-coordinate
+ # spatial coordinate index,
+ # i.e. angular spatial coordinate
+ # indices are A_min..N
+
+delta := array(1..N, 1..N, identity); # Kronecker delta
+
+# (x,y,z) coordinates
+x_xyz := array(1..N);
+x_xyz[1] := xx;
+x_xyz[2] := yy;
+x_xyz[3] := zz;
+x_xyz_list := map(op, [entries(x_xyz)]);
+
+# (r,rho,sigma) coordinates
+# we hard-wire that r is the first coordinate and remaining ones are angular
+x_rrs := array(1..N);
+x_rrs[1] := r;
+x_rrs[2] := rho;
+x_rrs[3] := sigma;
+x_rrs_list := map(op, [entries(x_rrs)]);
+
+# list of all coordinates (for sorting things into "nice" order)
+x_all_list := [op(x_xyz_list), op(x_rrs_list)];
+
+NULL;
+end proc;
+
+###############################################################################
+
+#
+# This function sets up the coordinate-transformation gfas
+#
+setup_coord_gfas :=
+proc()
+global
+@include "coords.minc";
+
+# 1st ($X^u_i$) and 2nd ($X^u_{ij}$) partial derivs
+# of (r,rho,sigma) wrt (xx,yy,zz)
+make_gfa('X_ud', {inert,fnd}, [A_min..N, 1..N], none);
+make_gfa('X_udd', {inert,fnd}, [A_min..N, 1..N, 1..N], symmetric3_23);
+
+NULL;
+end proc;
diff --git a/src/maple/coords.minc b/src/maple/coords.minc
new file mode 100644
index 0000000..532bfa0
--- /dev/null
+++ b/src/maple/coords.minc
@@ -0,0 +1,18 @@
+# coords.minc -- global variables for coordinates
+# $Id$
+
+# coordinates
+N, # number of spatial coordinates
+N_ang, # number of angular spatial coordinates
+A_min, # min angular-spatial-coordinate coordinate index
+delta, # Kronecker delta
+x_xyz, x_xyz_list, # array and list of (xx,yy,zz) coordinates
+x_rrs, x_rrs_list, # array and list of (r,rho,sigma) coordinates
+r, rho, sigma, # individual patch coordinates
+xx, yy, zz, # Cartesian coordinates
+x_all_list, # list of all coordinates
+
+# coordinate-transformation gfas
+X_ud, X_ud__fnd,
+X_udd, X_udd_fnd
+ # no comma
diff --git a/src/maple/gfa.maple b/src/maple/gfa.maple
new file mode 100644
index 0000000..e0733b6
--- /dev/null
+++ b/src/maple/gfa.maple
@@ -0,0 +1,315 @@
+# gfa.maple -- low-level gridfn array "database"
+# $Id$
+#
+# ***** gridfn arrays *****
+# ***** functional dependences *****
+# ***** gridfn-array sets *****
+# make_gfa - construct a new gridfn array
+# assert_fnd_exists - assert (verify) that a given gridfn array fnd form exists
+# fnd_name - Maple name of functional-dependence form
+#
+# index/symmetric3_23 - Maple indexing fn for rank 3 array sym on axes 23
+#
+# print_symmetric3_23 - prettyprint a symmetric3_23 array
+#
+
+###############################################################################
+
+#
+# ***** gridfn arrays *****
+#
+
+#
+# The basic data objects for the Maple code are grid function arrays,
+# a.k.a. gridfn arrays, such as the metric g_dd, the Christoffel symbols
+# Gamma_udd, etc etc.
+#
+# By convention, the name of a nonscalar geometric object ends with an
+# underscore ("_") followed by a sequence of "d" or "u" specifying whether
+# the indices are down (covariant) or up (contravariant). Thus, for
+# example, g_dd is the covariant metric, while g_uu is the inverse metric.
+#
+# Scalar geometric objects are represented by Maple scalars, while
+# nonscalar objects (often but not necessarily tensors) are represented
+# by Maple arrays, using Maple indexing functions for any symmetries
+# or sparsity. For example, the metric g_dd is represented by an
+# array
+# array(1..N, 1..N, symmetric)
+# By convention, for a symmetric array of this type, we only compute
+# the upper triangle.
+#
+
+###############################################################################
+
+#
+# ***** functional dependences *****
+#
+
+#
+# Within the Maple code, we distinguish between two distinct forms
+# that a gridfn array may take:
+# - The "inert" form of a gridfn array represents (only) the symmetries,
+# zero (or nonzero but constant) elements, and indexing semantics of
+# the gridfn array.
+# - The "functional dependence" form of a gridfn array additionally
+# represents how that gridfn array is computed from other gridfn
+# arrays.
+#
+# For example, consider the inverse metric g_uu . It's inert form
+# is a 3*3 symmetric array, with all elements being unassigned symbols
+# (i.e. evaluating to themselves, as is usual in Maple). It's functional
+# dependence form, in contrast, is also a 3*3 symmetric array, but its
+# elements are expressions giving the inverse-metric elements in terms
+# of components of the metric and/or the metric determinant.
+#
+# Any given gridfn array may exist in either or (the usual case) both
+# an inert form and a functional dependence forms. There may also be
+# multiple functional dependence forms, corresponding to alternative
+# formulas for its computation.)
+#
+# The inert form of a gridfn array is represented by a Maple variable
+# or array (see below) of the same name (eg g_uu ), unassigned except
+# for any symmetries and/or zero components the gridfn array has.
+#
+# The functional dependence form of a gridfn array is represented by
+# a Maple array (see below) with a name formed by appending "__fnd"
+# to the gridfn array's name (eg g_uu__fnd ); this array has the same
+# symmetries as the corresponding inert Maple array, but only those
+# elements needed will be assigned. If there are multiple functional
+# dependence forms, by convention each uses a name beginning in "__"
+# and ending in "fnd", eg H__ps_fnd .
+#
+# The PDE compiler normally generates C code referring to the (inert)
+# gridfn array names, so these should be valid C variable names.
+#
+
+#
+# To allow some run-time error checking, we keep a global table
+# of all the gridfn array functional-dependence forms:
+# `gfa/fnd_table`
+# ... table index is gfa__fnd_name
+# ... table value is true
+# assert_fnd_exists() uses this table to verify that a given gridfn
+# array functional-dependence form does indeed exist.
+#
+
+#
+# Functional-dependence forms may be multilevel, as in (eg)
+# g_dd__Kerr_fnd__ps_fnd . cat() may be used to combine names in
+# this manner.
+#
+
+###############################################################################
+
+#
+# This function makes the Maple representation of a single gridfn array.
+#
+# Arguments:
+# gfa_name = (out*) The Maple variable to be created for the gridfn
+# array. This argument should be passed as an
+# unevaluated (quoted) name.
+# fnd_name_set = (in) A set of functional dependence types the gridfn
+# array is to have. The set's members should be
+# drawn from from {inert, fnd, ...} as appropriate,
+# where the special value inert refers to the
+# inert form.
+# index_bounds = (in) The list of index bounds for the gridfn array,
+# index_fn = (in) The name of the array indexing function, if any, or
+# 'none' if the gridfn array has no symmetries and thus
+# should be represented by a Maple array with the standard
+# "ordinary array" indexing semantics.
+#
+# Thus, for example, the calls
+# make_gfa('g', {inert, fnd}, [], none);
+# make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric);
+# make_gfa('R_dd', {inert, fnd}, [1..N, 1..N], symmetric);
+# make_gfa('Gamma_udd', {inert, fnd}, [1..N, 1..N, 1..N], symmetric3_23);
+# would set up the Maple representations of the metric determinant,
+# the metric, the Ricci tensor, and the Christoffel symbols respectively.
+#
+make_gfa :=
+proc(gfa_name::name,
+ fnd_name_set::set(name),
+ index_bounds::list(integer..integer),
+ index_fn::{identical(none), name}) # note "name" here, not "procedure",
+ # since we'll get the name of the
+ # indexing fn, which may not have
+ # been assigned its procedure yet
+global N, `gfa/fnd_table`;
+local fnd_name, var_name, index_fn_seq;
+
+# firewall checks
+if (not type(N, integer))
+ then ERROR("must set up coordinates first!");
+end if;
+if ((index_bounds = []) and (index_fn <> 'none'))
+ then ERROR("not meaningful to specify a symmetry function for a scalar!",
+ " gfa_name=",gfa_name);
+end if;
+
+# create the Maple variable(s)
+ for fnd_name in fnd_name_set
+ do
+ var_name := Maple_name(gfa_name, fnd_name);
+ if (assigned(`gfa/fnd_table`[eval(var_name,1)]))
+ then ERROR("duplicate gfa/fnd definition!",
+ " gfa_name=", gfa_name,
+ " fnd_name_set=", fnd_name_set);
+ end if;
+ `gfa/fnd_table`[eval(var_name,1)] := true;
+ if (index_bounds = [])
+ then # scalar
+ unassign( eval(var_name,1) );
+ else # array
+ if (index_fn = 'none')
+ then index_fn_seq := NULL;
+ else index_fn_seq := index_fn;
+ end if;
+ assign(
+ eval(var_name,1) = array(index_fn_seq, op(index_bounds))
+ );
+ end if;
+ end do;
+
+NULL;
+end proc;
+
+###############################################################################
+
+#
+# This function checks that a given gridfn array functional-dependence
+# form exists. If so, this function is a no-op; if not, this function
+# calls ERROR(). Hence after calling this function, the caller can assume
+# that the given functional-dependence form does indeed exist.
+#
+# This function may be called with either 1 or 2 arguments (the two
+# types of calls are equivalent).
+#
+# Arguments (1-argument form):
+# gfa_name = The name of the fnd-form to verify, eg 'g_dd.fnd'.
+# If this isn't inert, it should be quoted so this function
+# gets the name, not whatever value may have been assigned to it.
+#
+# Arguments (2-argument form):
+# gfa_name = The name of the gridfn array itself, eg `g_dd`.
+# fnd_name = The name of the fnd itself, eg `fnd`.
+#
+assert_fnd_exists :=
+proc(gfa_name::name, fnd_name::name)
+global `gfa/fnd_table`;
+local var_name;
+
+var_name := Maple_name(gfa_name, args[2..nargs]);
+
+if (not assigned(`gfa/fnd_table`[eval(var_name,1)]))
+ then ERROR("functional-dependence form doesn't exist!",
+ "var_name=",var_name);
+end if;
+end proc;
+
+###############################################################################
+
+#
+# This function computes the Maple name corresponding to a given
+# gridfn array and functional dependence form.
+#
+# Arguments:
+# gfa_name = The name of the gridfn array itself, eg `g_dd`.
+# fnd_name = (optional)
+# The name of the functional-dependence form, eg `fnd`,
+# or the special name 'inert'. If this argument is omitted
+# it defaults to 'inert'.
+#
+Maple_name :=
+proc(gfa_name::name, fnd_name::name)
+if ((nargs = 1) or (fnd_name = 'inert'))
+ then return gfa_name;
+ else return cat(gfa_name, "__", fnd_name);
+end if;
+end proc;
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This function is a Maple indexing function for a rank 3 table/array
+# which is symmetric in the 2nd and 3rd indices. This function is
+# called automagically by Maple whenever a component of such an array
+# is referenced.
+#
+# The ordering is done using the Maple sort() function, which sorts
+# numbers numerically, strings lexicographically, and "other things"
+# by machine address.
+#
+# Arguments:
+# ilist = (in) The indices to be used. There must be exactly 3 of these.
+# tab = (in out) The table/array to be indexed. This uses any remaining
+# indexing methods.
+# vlist = (optional)
+# If present, this is a list of the values to be assigned to
+# the (lvalue) table entry.
+#
+`index/symmetric3_23` :=
+proc(ilist::list, tab::table, vlist::list)
+local k, i, j, index_seq;
+
+if (not (nops(ilist) = 3))
+ then ERROR(`must have exactly 3 indices!`);
+end if;
+
+# ... note we need eval() here to get actual integers,
+# rather than variable names whose values are integers
+k := eval(ilist[1]);
+i := eval(ilist[2]);
+j := eval(ilist[3]);
+
+index_seq := k, op(sort([i,j]));
+
+if (nargs = 2)
+ then return tab[index_seq]; # rvalue
+ else tab[index_seq] := op(vlist); # lvalue
+end if;
+end proc;
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This function prettyprints a symmetric3_23 array
+#
+# Arguments:
+# A = (in) The array to be prettyprinted
+#
+# Results:
+# There is no explicit result; the array is prettyprinted as a side effect.
+#
+print_symmetric3_23 :=
+proc(A::array)
+local bounds, k, i, j, M23;
+
+if (op(1, eval(A)) <> 'symmetric3_23')
+ then ERROR(`can only print symmetric3_23 arrays`);
+end if;
+
+bounds := op(2, eval(A));
+
+M23 := array(bounds[2..3], symmetric);
+
+ for k in $(bounds[1]) do
+
+ for i in $(bounds[2]) do
+ for j in $(bounds[3]) do
+ if (i <= j) # only need upper triangle in (i,j)
+ then M23[i,j] := A[k,i,j];
+ end if;
+ end do;
+ end do;
+
+ printf("[%d] = \n", k);
+ print(M23);
+ end do;
+
+NULL;
+end proc;
diff --git a/src/maple/setup.maple b/src/maple/setup.maple
new file mode 100644
index 0000000..8602d4e
--- /dev/null
+++ b/src/maple/setup.maple
@@ -0,0 +1,24 @@
+# setup.maple -- top-level Maple input file to set up PDE compiler system
+# $Id$
+
+#
+# Note not all files here actually *use* the preprocessor, but for
+# convenience we preprocess everything (this way we don't have to keep
+# track of which files should be used as .maple and which as .mm).
+#
+
+#
+# Note paths are evaluated relative to where read command gets executed,
+# not relative to this directory, so we use ../maple/ to try and avoid
+# ambiguities...
+#
+
+read "../maple/util.mm";
+
+read "../maple/coords.mm";
+setup_coords();
+
+read "../maple/Diff.mm";
+
+read "../maple/gfa.mm";
+setup_coord_gfas();
diff --git a/src/maple/util.maple b/src/maple/util.maple
new file mode 100644
index 0000000..2e574df
--- /dev/null
+++ b/src/maple/util.maple
@@ -0,0 +1,234 @@
+# util.maple -- Maple utility functions
+# $Id$
+
+#
+# msum - multiple-index sum()
+# arctan_xy - 2-argument arctan() function
+#
+# sort_var_list - sort a list of variables into a canonical order
+# lexorder_vars - lexorder() for "interesting" variables
+#
+# gensym - generate a "new" symbol name
+# gensym/init - reinitialize the gensym() global naming counter
+#
+# saveit - optionally save variables for debugging
+#
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This function is identical to the standard Maple library function sum()
+# except that it allows multiple summation indices. It's usage is best
+# illustrated by a series of examples:
+#
+# msum('f', 'k'=m..n)
+# = sum('f', 'k'=m..n)
+#
+# msum('f', 'k1'=m1..n1, 'k2'=m2..n2)
+# = sum( 'sum('f', 'k1'=m1..n1)' , 'k2'=m2..n2 )
+#
+# msum('f', 'k1'=m1..n1, 'k2'=m2..n2, 'k3'=m3..n3)
+# = sum(
+# 'sum( 'sum('f', 'k1'=m1..n1)' , 'k2'=m2..n2 )'
+# ,
+# 'k3'=m3..n3
+# )
+#
+# Note, in particular, that the summation indices are nested with the
+# innermost summation appearing first. This is the same nesting convention
+# as Fortran I/O statement's implicit do loops use.
+#
+# Bugs:
+# - The quoting in the code is pretty obscure...
+#
+msum :=
+proc(fn::algebraic) # varargs
+local expr, sum_index;
+
+if (nargs < 2)
+ then ERROR("must have two or more arguments")
+end if;
+
+expr := fn;
+
+ for sum_index in [args[2..nargs]]
+ do
+ # loop invariant:
+ # msum(fn, <args processed so far>) = sum(expr, sum_index)
+ expr; sum_index; expr:= 'sum'(''%%'', ''%'');
+ end do;
+
+return eval(expr);
+end proc;
+
+###############################################################################
+
+#
+# This function computes the 4-quadrant arctangent of arctan(y/x) .
+# It's "just" a wrapper around the Maple 2-argument arctan function,
+# providing the what-I-find-natural (x,y) argument
+#
+# Arguments:
+# (x,y) = (In) The x,y coordinates.
+#
+arctan_xy :=
+proc(x::algebraic, y::algebraic)
+arctan(y,x); # !!! reversedargument order !!!
+end proc;
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function sorts a list of names (in practice, a list of coordinates)
+# into a canonical order.
+#
+sort_var_list :=
+proc(var_list::list(name))
+option remember; # performance optimization
+
+# only get to here the first time we see a given variable list
+# (i.e. if it's not in the remember table)
+RETURN( sort(var_list, lexorder_vars) );
+end;
+
+################################################################################
+
+#
+# This function defines a lexical ordering on variable names.
+#
+# Arguments:
+# x, y = Two names to be compared.
+#
+# Results:
+# This function returns true iff a < b, false otherwise.
+#
+lexorder_vars :=
+proc(x::name, y::name)
+option remember; # performance optimization
+global
+ @include "coords.minc";
+
+if (member(x, x_all_list, 'xposn') and member(y, x_all_list, 'yposn'))
+ then return evalb(xposn < yposn);
+ else return lexorder(x, y);
+fi;
+
+end;
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This procedure generates a (presumably) new symbol name of the
+# form
+# cat(base_name,count)
+# where base_name is passed as an (optional, with a default value
+# supplied if it's omitted) argument, and count is the current
+# value of the global variable
+# `gensym/counter`
+# which is postincremented after being used. If the counter is
+# unassigned on entry to this function, `gensym/init`() is first
+# called to initialize it.
+#
+# `gensym/init`() may be used to (re)initialize the counter as
+# desired, or `gensym/counter` may be examined or altered directly
+# to otherwise modify the name sequence generated by this function.
+#
+# Arguments:
+# base_name = (in) (optional) The base name for the name generation.
+# This argument defaults to `temp_` ifo
+# omitted.
+#
+# Bugs:
+# No checks are made for whether the names in question are indeed
+# unassigned.
+#
+gensym :=
+proc(opt_base_name::string)
+global `gensym/counter`; # in out
+local base_name, tn;
+
+if (nargs >= 1)
+ then base_name := opt_base_name;
+ else base_name := '`temp_`';
+end if;
+
+# force an initialization if it hasn't been done before
+if (not assigned(`gensym/counter`))
+ then `gensym/init`();
+end if;
+
+tn := cat(base_name, `gensym/counter`);
+`gensym/counter` := `gensym/counter` + 1;
+
+tn; return '%';
+end proc;
+
+###############################################################################
+
+#
+# This procedure (re)initializes the gensym() global naming counter.
+#
+# Arguments:
+# initial_counter = (in) (optional)
+# The naming counter value to (re)initialize with.
+# This argument defaults to 1 if omitted.
+#
+`gensym/init` :=
+proc(opt_initial_counter::integer)
+global `gensym/counter`; # in out
+local initial_counter;
+
+if (nargs >= 1)
+ then initial_counter := opt_initial_counter;
+ else initial_counter := 1;
+end if;
+
+`gensym/counter` := initial_counter;
+NULL;
+end proc;
+
+###############################################################################
+###############################################################################
+###############################################################################
+
+#
+# This function optionally saves a local variable to a global variable
+# for later examination. This is useful for debugging purposes, since
+# Maple doesn't have a proper debugger :-( .
+#
+# Arguments (explicit):
+# fn_name = (in) The calling function or subsystem's name.
+# expr = (in) The temporary result to save.
+# name_postfix... = (in) (varargs)
+# One or more arguments which this function
+# concatenates to form the name postfix for
+# the global variable.
+#
+# Arguments (implicit, as global variables)
+# `saveit/level` = (in) (optional) If this name is assigned (as determined
+# by the assigned() function), this
+# function does the saving. Otherwise,
+# this function is a nop.
+#
+saveit :=
+proc(n::integer, fn_name::string, label::string, expr::anything)
+global `saveit/level`;
+local save_name;
+
+if ( assigned(`saveit/level`)
+ and type(`saveit/level`, integer)
+ and (`saveit/level` >= n) )
+ then
+ save_name := cat(fn_name,`/`,label);
+ lprint(printf(` --> \"%s\"`, save_name));
+ assign( eval(save_name,1) = expr );
+end if;
+
+NULL;
+end proc;