* maple-generated code wants (xx,yy,zz) = *local* xyz coords

  (compute this on the fly at each grid point)
* remove X_ud and X_udd deriv coeffs from grid -- we have to compute
  them anyway at each grid point, so no need to store them in the grid

*** with these changes, code now correctly evaluates H(h)
*** for Schw/EF with the patch center at (0.5,0.5,0.0)


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@584 f88db872-0e4f-0410-b76b-b9085cfa78c5

6 files changed, 109 insertions, 636 deletions

diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index d2dd4a7..b6174eb 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -13,6 +13,9 @@
 // The machine-generated code uses the following variables:
 //	patch& p		// current patch
 //	int irho,isigma		// current generic integer coords within patch
+//	fp xx, yy, zz;		// local (x,y,z) coords of this grid point
+//	fp X_ud_[12][123]			// 1st derivative coefficients
+//	fp X_udd_[12]{11,12,13,22,23,33}	// 2nd derivative coefficients
 //
 
 //
@@ -44,30 +47,6 @@
 //
 // nominal-grid gridfns
 //
-#define xx	p.gridfn(nominal_gfns::gfn__xx, irho,isigma)
-#define yy	p.gridfn(nominal_gfns::gfn__yy, irho,isigma)
-#define zz	p.gridfn(nominal_gfns::gfn__zz, irho,isigma)
-
-#define X_ud_11	p.gridfn(nominal_gfns::gfn__X_ud_11, irho,isigma)
-#define X_ud_12	p.gridfn(nominal_gfns::gfn__X_ud_12, irho,isigma)
-#define X_ud_13	p.gridfn(nominal_gfns::gfn__X_ud_13, irho,isigma)
-#define X_ud_21	p.gridfn(nominal_gfns::gfn__X_ud_21, irho,isigma)
-#define X_ud_22	p.gridfn(nominal_gfns::gfn__X_ud_22, irho,isigma)
-#define X_ud_23	p.gridfn(nominal_gfns::gfn__X_ud_23, irho,isigma)
-
-#define X_udd_111	p.gridfn(nominal_gfns::gfn__X_udd_111, irho,isigma)
-#define X_udd_112	p.gridfn(nominal_gfns::gfn__X_udd_112, irho,isigma)
-#define X_udd_113	p.gridfn(nominal_gfns::gfn__X_udd_113, irho,isigma)
-#define X_udd_122	p.gridfn(nominal_gfns::gfn__X_udd_122, irho,isigma)
-#define X_udd_123	p.gridfn(nominal_gfns::gfn__X_udd_123, irho,isigma)
-#define X_udd_133	p.gridfn(nominal_gfns::gfn__X_udd_133, irho,isigma)
-#define X_udd_211	p.gridfn(nominal_gfns::gfn__X_udd_211, irho,isigma)
-#define X_udd_212	p.gridfn(nominal_gfns::gfn__X_udd_212, irho,isigma)
-#define X_udd_213	p.gridfn(nominal_gfns::gfn__X_udd_213, irho,isigma)
-#define X_udd_222	p.gridfn(nominal_gfns::gfn__X_udd_222, irho,isigma)
-#define X_udd_223	p.gridfn(nominal_gfns::gfn__X_udd_223, irho,isigma)
-#define X_udd_233	p.gridfn(nominal_gfns::gfn__X_udd_233, irho,isigma)
-
 #define g_dd_11	p.gridfn(nominal_gfns::gfn__g_dd_11, irho,isigma)
 #define g_dd_12	p.gridfn(nominal_gfns::gfn__g_dd_12, irho,isigma)
 #define g_dd_13	p.gridfn(nominal_gfns::gfn__g_dd_13, irho,isigma)
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 0abc9f7..fe64dce 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -273,7 +273,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
 		fp r_plus, r_minus;
 		#include "ellipsoid.c"
 
-		// exactly one of the solutions (=def r) should be positive
+		// exactly one of the solutions (call it r) should be positive
 		fp r;
 		if      ((r_plus > 0.0) && (r_minus < 0.0))
 		   then r = r_plus;
diff --git a/src/gr/ellipsoid.out b/src/gr/ellipsoid.out
index 2e77483..4a0ee01 100644
--- a/src/gr/ellipsoid.out
+++ b/src/gr/ellipsoid.out
@@ -4,7 +4,7 @@
  <____ ____>  Waterloo Maple Inc.
       |       Type ? for help.
 # ellipsoid.maple -- compute equations for offset ellipsoid setup
-# $Header: /numrelcvs/ThornburgCVS/AHFinderDirect/src/gr/ellipsoid.maple,v 1.1 2002/04/22 11:02:27 jthorn Exp $
+# $Header: /numrelcvs/ThornburgCVS/AHFinderDirect/src/gr/ellipsoid.maple,v 1.2 2002/04/22 15:10:32 jthorn Exp $
 > 
 #
 # ellipsoid has center (A,B,C), radius (a,b,c)
@@ -33,492 +33,6 @@
 
 > 
 > read "../maple/util.mm";
-msum := proc(fn::algebraic)
-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
-        expr; sum_index; expr := 'sum'(''%%'', ''%'')
-    end do;
-    return eval(expr)
-end proc
-
-      arctan_xy := proc(x::algebraic, y::algebraic) arctan(y, x) end proc
-
-            ssqrt := proc(x::algebraic) sqrt(x, 'symbolic') end proc
-
-indices_in_order :=
-
-    proc(T::table) return sort([indices(T)], lexorder_integer_list) end proc
-
-lexorder_integer_list := proc(list1::list(numeric), list2::list(numeric))
-local len1, len2, k;
-    len1 := nops(list1);
-    len2 := nops(list2);
-    for k to min(len1, len2) do
-        if list1[k] < list2[k] then return true
-        elif list2[k] < list1[k] then return false
-        end if
-    end do;
-    return evalb(len1 < len2)
-end proc
-
-sort_var_list := proc(var_list::list(name))
-global lexorder_vars;
-option remember;
-    return sort(var_list, lexorder_vars)
-end proc
-
-lexorder_vars := proc(x::name, y::name)
-local xposn, yposn;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
-rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set;
-option remember;
-    if member(x, xy_all_list, 'xposn') and member(y, xy_all_list, 'yposn')
-    then return evalb(xposn < yposn)
-    else return lexorder(x, y)
-    end if
-end proc
-
-gensym := proc(opt_base_name::string)
-local base_name, tn;
-global `gensym/counter`;
-    if 1 <= nargs then base_name := opt_base_name
-    else base_name := 'temp_'
-    end if;
-    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
-
-gensym/init := proc(opt_initial_counter::integer)
-local initial_counter;
-global `gensym/counter`;
-    if 1 <= nargs then initial_counter := opt_initial_counter
-    else initial_counter := 1
-    end if;
-    `gensym/counter` := initial_counter;
-    NULL
-end proc
-
-saveit := proc(
-n::integer, fn::{procedure, string}, label::string, expr::anything)
-local save_name;
-global `saveit/level`;
-    if assigned(`saveit/level`) and type(`saveit/level`, integer) and
-    n <= `saveit/level` then
-        save_name := cat(convert(fn, string), "/", label);
-        printf("      --> `%s`\n", save_name);
-        assign(convert(eval(save_name, 1), name) = expr)
-    end if;
-    NULL
-end proc
-
-> read "../maple/codegen2.mm";
-codegen2 := proc(expr_in::{algebraic, list(algebraic)},
-lhs_name::{name, list(name)}, output_file_name::string)
-local expr, expr_temps, input_set, output_set, expr_cost;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
-rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
-fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd;
-    printf("codegen2(%a) --> \"%s\"\n", lhs_name, output_file_name);
-    expr := expr_in;
-    saveit(10, procname, "input", expr);
-    printf("   convert --> equation list\n");
-    expr := cvt_to_eqnlist(expr, lhs_name);
-    saveit(10, procname, "eqnlist", expr);
-    printf("   optimizing computation sequence\n");
-    expr := [codegen[optimize](expr)];
-    saveit(10, procname, "optimize", expr);
-    printf("   find temporary variables\n");
-    expr_temps := temps_in_eqnlist(expr, lhs_name);
-    saveit(10, procname, "temps", expr_temps);
-    printf(
-        "   convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc\n")
-        ;
-    expr := fix_Diff(expr);
-    saveit(10, procname, "fix_Diff", expr);
-    input_set := deindex_names(
-        (indets(map(rhs, expr), name) minus {op(expr_temps)}) minus
-        xy_all_set);
-    output_set :=
-        deindex_names({op(map(lhs, expr))} minus {op(expr_temps)});
-    printf("   convert R_dd[2,3] --> R_dd_23 etc\n");
-    expr := unindex_names(expr);
-    saveit(10, procname, "unindex", expr);
-    expr_cost := codegen[cost](expr);
-    printf("   convert p/q --> RATIONAL(p/q)\n");
-    expr := fix_rationals(expr);
-    saveit(10, procname, "fix_rationals", expr);
-    printf("   writing C code\n");
-    ftruncate(output_file_name);
-    fprintf(output_file_name, "/*\n");
-    fprintf(output_file_name, " * inputs = %a\n", input_set);
-    fprintf(output_file_name, " * outputs = %a\n", output_set);
-    fprintf(output_file_name, " * cost = %a\n", expr_cost);
-    fprintf(output_file_name, " */\n");
-    print_name_list_dcl(expr_temps, "fp", output_file_name);
-    codegen[C](expr, filename = output_file_name);
-    NULL
-end proc
-
-cvt_to_eqnlist := proc(expr::{algebraic, array, list({algebraic, array})},
-lhs_name::{name, list(name)})
-    if type(expr, array) and type(lhs_name, name) then return map(
-        proc(ii) return lhs_name[op(ii)] = expr[op(ii)] end proc,
-        indices_in_order(expr))
-    end if;
-    if type(expr, algebraic) and type(lhs_name, name) then
-        return [lhs_name = expr]
-    end if;
-    if type(expr, list({algebraic, array})) and type(lhs_name, list(name))
-    then return zip(op@cvt_to_eqnlist, expr, lhs_name)
-    end if;
-    error
-        "unknown type for expression!\n   expr=%1\n   whattype(expr)=%2\n",
-        expr, whattype(expr)
-end proc
-
-fix_Diff := proc(
-expr::{algebraic, name = algebraic, list({algebraic, name = algebraic})})
-local nn, k, base, power, fn, fn_args_list, Darg, Dvars;
-global `fix_Diff/remap_table`;
-    if type(expr, list) then return map(fix_Diff, expr) end if;
-    if type(expr, name = algebraic) then
-        return lhs(expr) = fix_Diff(rhs(expr))
-    end if;
-    nn := nops(expr);
-    if type(expr, `+`) then
-        return sum('fix_Diff(op(k, expr))', 'k' = 1 .. nn)
-    end if;
-    if type(expr, `*`) then
-        return product('fix_Diff(op(k, expr))', 'k' = 1 .. nn)
-    end if;
-    if type(expr, `^`) then
-        base := op(1, expr);
-        power := op(2, expr);
-        return fix_Diff(base)^power
-    end if;
-    if type(expr, function) and op(0, expr) <> 'Diff' then
-        fn := op(0, expr);
-        fn_args_list := [op(expr)];
-        fn;
-        return '%'(op(map(fix_Diff, fn_args_list)))
-    end if;
-    if type(expr, function) and op(0, expr) = 'Diff' then
-        Darg := op(1, expr);
-        Dvars := [op(2 .. nn, expr)];
-        if assigned(`fix_Diff/remap_table`[op(Dvars)]) then
-            `fix_Diff/remap_table`[op(Dvars)]; return '%'(Darg)
-        else error "don't know how to remap Diff() call!\n   Darg = %1\n\
-               Dvars = %2\n", Darg, Dvars
-        end if
-    end if;
-    return expr
-end proc
-
-                    fix_Diff/remap_table[rho] := PARTIAL_RHO
-
-                  fix_Diff/remap_table[sigma] := PARTIAL_SIGMA
-
-               fix_Diff/remap_table[rho, rho] := PARTIAL_RHO_RHO
-
-             fix_Diff/remap_table[rho, sigma] := PARTIAL_RHO_SIGMA
-
-           fix_Diff/remap_table[sigma, sigma] := PARTIAL_SIGMA_SIGMA
-
-                     fix_Diff/remap_table[xx] := PARTIAL_X
-
-                     fix_Diff/remap_table[yy] := PARTIAL_Y
-
-                     fix_Diff/remap_table[zz] := PARTIAL_Z
-
-                   fix_Diff/remap_table[xx, xx] := PARTIAL_XX
-
-                   fix_Diff/remap_table[xx, yy] := PARTIAL_XY
-
-                   fix_Diff/remap_table[xx, zz] := PARTIAL_XZ
-
-                   fix_Diff/remap_table[yy, yy] := PARTIAL_YY
-
-                   fix_Diff/remap_table[yy, zz] := PARTIAL_YZ
-
-                   fix_Diff/remap_table[zz, zz] := PARTIAL_ZZ
-
-temps_in_eqnlist := proc(
-expr::list(name = algebraic), result_name::{name, list(name), set(name)})
-    return remove(is_result, map(lhs, expr), result_name)
-end proc
-
-is_result := proc(try_name::name, result_name_in::{name, list(name), set(name)}
-)
-local result_name, rn;
-    if type(result_name_in, name) then result_name := {result_name_in}
-    else result_name := result_name_in
-    end if;
-    for rn in result_name do
-        if try_name = rn then return true
-        elif type(try_name, indexed) and op(0, try_name) = rn then
-            return true
-        end if
-    end do;
-    return false
-end proc
-
-deindex_names := proc(
-expr::{function, name, list({function, name}), set({function, name})})
-local fn, fn_args_list;
-    if type(expr, {set, list}) then return map(deindex_names, expr) end if;
-    if type(expr, function) then
-        fn := op(0, expr);
-        fn_args_list := [op(expr)];
-        fn;
-        return '%'(op(map(deindex_names, fn_args_list)))
-    end if;
-    if type(expr, indexed) then return op(0, expr) end if;
-    if type(expr, {name, numeric}) then return expr end if;
-    error "expr has unknown type!\nwhattype(expr)=%1\n", whattype(expr)
-end proc
-
-unindex_names := proc(expr::{algebraic, name = algebraic,
-list({algebraic, name = algebraic}), set({algebraic, name = algebraic})})
-local nn, k, base, power, fn, fn_args_list, base_name, index_seq;
-    if type(expr, {set, list}) then return map(unindex_names, expr) end if;
-    if type(expr, `=`) then
-        return unindex_names(lhs(expr)) = unindex_names(rhs(expr))
-    end if;
-    nn := nops(expr);
-    if type(expr, `+`) then
-        return sum('unindex_names(op(k, expr))', 'k' = 1 .. nn)
-    end if;
-    if type(expr, `*`) then
-        return product('unindex_names(op(k, expr))', 'k' = 1 .. nn)
-    end if;
-    if type(expr, `^`) then
-        base := op(1, expr);
-        power := op(2, expr);
-        return unindex_names(base)^power
-    end if;
-    if type(expr, function) then
-        fn := op(0, expr);
-        fn_args_list := [op(expr)];
-        fn;
-        return '%'(op(map(unindex_names, fn_args_list)))
-    end if;
-    if type(expr, indexed) then
-        base_name := op(0, expr);
-        index_seq := op(expr);
-        return cat(base_name, "_", index_seq)
-    end if;
-    if type(expr, {name, numeric}) then return expr end if;
-    error "expr has unknown type!\nwhattype(expr)=%1\n", whattype(expr)
-end proc
-
-fix_rationals := proc(
-expr::{algebraic, name = algebraic, list({algebraic, name = algebraic})})
-local nn, k, expr_sign, expr_abs, base, power, fbase, fpower, fn,
-fn_args_list, int_factors, nonint_factors, num, den, mult;
-    if type(expr, list) then return map(fix_rationals, expr) end if;
-    if type(expr, name = algebraic) then
-        return lhs(expr) = fix_rationals(rhs(expr))
-    end if;
-    if type(expr, function) then
-        fn := op(0, expr);
-        if fn <> 'RATIONAL' then
-            fn_args_list := [op(expr)];
-            fn;
-            return '%'(op(map(fix_rationals, fn_args_list)))
-        end if
-    end if;
-    nn := nops(expr);
-    if type(expr, `+`) then
-        return sum('fix_rationals(op(k, expr))', 'k' = 1 .. nn)
-    end if;
-    if type(expr, `*`) then
-        int_factors, nonint_factors := selectremove(type, expr, integer);
-        if 0 < nops(int_factors) then return op(1, int_factors)*product(
-            'fix_rationals(op(k, nonint_factors))',
-            'k' = 1 .. nops(nonint_factors))
-        else return product('fix_rationals(op(k, expr))', 'k' = 1 .. nn)
-        end if
-    end if;
-    if type(expr, `^`) then
-        base := op(1, expr);
-        power := op(2, expr);
-        fbase := fix_rationals(base);
-        if type(power, integer) then fpower := power
-        else fpower := fix_rationals(power)
-        end if;
-        return fbase^fpower
-    end if;
-    if type(expr, integer) then return 'RATIONAL'(expr, 1) end if;
-    if type(expr, fraction) then
-        num := op(1, expr); den := op(2, expr); return 'RATIONAL'(num, den)
-    end if;
-    if type(expr, float) then
-        mult := op(1, expr);
-        power := op(2, expr);
-        return fix_rationals(mult*10^power)
-    end if;
-    if type(expr, name) then return expr end if;
-    error "expr has unknown type!\nwhattype(expr)=%1\nexpr=%2\n",
-        whattype(expr), expr
-end proc
-
-print_name_list_dcl := proc(
-name_list::list({name, string}), name_type::string, file_name::string)
-local nn;
-    nn := nops(name_list);
-    if nn <= 10 then
-        map(convert, name_list, string);
-        ListTools[Join](%, ", ");
-        cat(op(%));
-        fprintf(file_name, "%s %s;\n", name_type, %);
-        NULL;
-        return
-    end if;
-    print_name_list_dcl([op(1 .. 10, name_list)], name_type, file_name);
-    print_name_list_dcl([op(11 .. nn, name_list)], name_type, file_name)
-end proc
-
-> 
-> [solve(eqn, r)];
-bytes used=1000572, alloc=851812, time=0.07
-         2  2              2  2              2  2
-[1/2 (2 a  b  zcos CW + 2 a  c  ycos BV + 2 b  c  xcos AU + 2 (
-
-       4  2          2              2  4          2
-    2 a  b  zcos CW c  ycos BV + 2 a  b  zcos CW c  xcos AU
-
-          2  4          2            4  2     2  2   2    2  4     2  2   2
-     + 2 a  c  ycos BV b  xcos AU - b  c  xcos  a  CW  - b  c  xcos  a  BV
-
-        4  4     2  2    4  2     2  2   2    2  4     2  2   2
-     + b  c  xcos  a  - a  c  ycos  b  CW  - a  c  ycos  b  AU
-
-        4  4     2  2    4  2     2  2   2    2  4     2  2   2
-     + a  c  ycos  b  - a  b  zcos  c  BV  - a  b  zcos  c  AU
-
-        4  4     2  2 1/2    /   2  2     2    2  2     2    2  2     2
-     + a  b  zcos  c )   )  /  (b  c  xcos  + a  c  ycos  + a  b  zcos ), 1/2 (
-                           /
-
-       2  2              2  2              2  2
-    2 a  b  zcos CW + 2 a  c  ycos BV + 2 b  c  xcos AU - 2 (
-
-       4  2          2              2  4          2
-    2 a  b  zcos CW c  ycos BV + 2 a  b  zcos CW c  xcos AU
-
-          2  4          2            4  2     2  2   2    2  4     2  2   2
-     + 2 a  c  ycos BV b  xcos AU - b  c  xcos  a  CW  - b  c  xcos  a  BV
-
-        4  4     2  2    4  2     2  2   2    2  4     2  2   2
-     + b  c  xcos  a  - a  c  ycos  b  CW  - a  c  ycos  b  AU
-
-        4  4     2  2    4  2     2  2   2    2  4     2  2   2
-     + a  c  ycos  b  - a  b  zcos  c  BV  - a  b  zcos  c  AU
-
-        4  4     2  2 1/2    /   2  2     2    2  2     2    2  2     2
-     + a  b  zcos  c )   )  /  (b  c  xcos  + a  c  ycos  + a  b  zcos )]
-                           /
-
-> map(simplify, %);
-bytes used=2000744, alloc=1376004, time=0.12
-bytes used=3001976, alloc=1769148, time=0.19
-   2  2            2  2            2  2              2  2  2
-[(a  b  zcos CW + a  c  ycos BV + b  c  xcos AU + (-a  b  c  (
-
-        2                      2                      2
-    -2 a  zcos CW ycos BV - 2 b  zcos CW xcos AU - 2 c  ycos BV xcos AU
-
-        2     2   2    2     2   2    2  2     2    2     2   2    2     2   2
-     + b  xcos  CW  + c  xcos  BV  - b  c  xcos  + a  ycos  CW  + c  ycos  AU
-
-        2  2     2    2     2   2    2     2   2    2  2     2  1/2    /
-     - a  c  ycos  + a  zcos  BV  + b  zcos  AU  - a  b  zcos ))   )  /  (
-                                                                     /
-
-     2  2     2    2  2     2    2  2     2     2  2            2  2
-    b  c  xcos  + a  c  ycos  + a  b  zcos ), (a  b  zcos CW + a  c  ycos BV
-
-        2  2              2  2  2      2                      2
-     + b  c  xcos AU - (-a  b  c  (-2 a  zcos CW ycos BV - 2 b  zcos CW xcos AU
-
-          2                    2     2   2    2     2   2    2  2     2
-     - 2 c  ycos BV xcos AU + b  xcos  CW  + c  xcos  BV  - b  c  xcos
-
-        2     2   2    2     2   2    2  2     2    2     2   2    2     2   2
-     + a  ycos  CW  + c  ycos  AU  - a  c  ycos  + a  zcos  BV  + b  zcos  AU
-
-        2  2     2  1/2    /   2  2     2    2  2     2    2  2     2
-     - a  b  zcos ))   )  /  (b  c  xcos  + a  c  ycos  + a  b  zcos )]
-                         /
-
-> [r_plus = %[1], r_minus = %[2]];
-            2  2            2  2            2  2              2  2  2
-[r_plus = (a  b  zcos CW + a  c  ycos BV + b  c  xcos AU + (-a  b  c  (
-
-        2                      2                      2
-    -2 a  zcos CW ycos BV - 2 b  zcos CW xcos AU - 2 c  ycos BV xcos AU
-
-        2     2   2    2     2   2    2  2     2    2     2   2    2     2   2
-     + b  xcos  CW  + c  xcos  BV  - b  c  xcos  + a  ycos  CW  + c  ycos  AU
-
-        2  2     2    2     2   2    2     2   2    2  2     2  1/2    /
-     - a  c  ycos  + a  zcos  BV  + b  zcos  AU  - a  b  zcos ))   )  /  (
-                                                                     /
-
-     2  2     2    2  2     2    2  2     2               2  2
-    b  c  xcos  + a  c  ycos  + a  b  zcos ), r_minus = (a  b  zcos CW
-
-        2  2            2  2              2  2  2      2
-     + a  c  ycos BV + b  c  xcos AU - (-a  b  c  (-2 a  zcos CW ycos BV
-
-          2                      2                    2     2   2
-     - 2 b  zcos CW xcos AU - 2 c  ycos BV xcos AU + b  xcos  CW
-
-        2     2   2    2  2     2    2     2   2    2     2   2    2  2     2
-     + c  xcos  BV  - b  c  xcos  + a  ycos  CW  + c  ycos  AU  - a  c  ycos
-
-        2     2   2    2     2   2    2  2     2  1/2    /
-     + a  zcos  BV  + b  zcos  AU  - a  b  zcos ))   )  /  (
-                                                       /
-
-     2  2     2    2  2     2    2  2     2
-    b  c  xcos  + a  c  ycos  + a  b  zcos )]
-
-> solnlist := [codegen[optimize](%)];
-                   2        2                                     2
-solnlist := [t1 = a , t2 = b , t3 = t1 t2, t5 = t3 zcos CW, t6 = c , t7 = t1 t6,
-
-                                                               2          2
-    t9 = t7 ycos BV, t10 = t2 t6, t12 = t10 xcos AU, t28 = xcos , t30 = CW ,
-
-            2                           2          2                          2
-    t33 = BV , t35 = t10 t28, t36 = ycos , t40 = AU , t42 = t7 t36, t43 = zcos ,
-
-    t48 = t3 t43, t49 = -2 t1 zcos CW ycos BV - 2 t2 zcos CW xcos AU
-
-     - 2 t6 ycos BV xcos AU + t2 t28 t30 + t6 t28 t33 - t35 + t1 t36 t30
-
-                                                                           1/2
-     + t6 t36 t40 - t42 + t1 t43 t33 + t2 t43 t40 - t48, t52 = (-t3 t6 t49)   ,
-
-                 1
-    t55 = ---------------, r_plus = (t5 + t9 + t12 + t52) t55,
-          t35 + t42 + t48
-
-    r_minus = (t5 + t9 + t12 - t52) t55]
-
-> 
-> ftruncate("ellipsoid.c");
-                            ftruncate("ellipsoid.c")
-
-> print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]),
-> 		    "fp", "ellipsoid.c");
-> codegen[C](solnlist, filename="ellipsoid.c");
-bytes used=4002168, alloc=1769148, time=0.26
+Error, unable to read `../maple/util.mm`
 > quit
-bytes used=4027884, alloc=1769148, time=0.27
+bytes used=129844, alloc=196572, time=0.03
diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh
index a6c1e11..16af52f 100644
--- a/src/gr/gfn.hh
+++ b/src/gr/gfn.hh
@@ -14,27 +14,9 @@ namespace nominal_gfns
 	  {
 enum	{
 	min_gfn = 1,
-	gfn__xx = min_gfn,
-	gfn__yy,
-	gfn__zz,
-	gfn__X_ud_11,
-	gfn__X_ud_12,
-	gfn__X_ud_13,
-	gfn__X_ud_21,
-	gfn__X_ud_22,
-	gfn__X_ud_23,
-	gfn__X_udd_111,
-	gfn__X_udd_112,
-	gfn__X_udd_113,
-	gfn__X_udd_122,
-	gfn__X_udd_123,
-	gfn__X_udd_133,
-	gfn__X_udd_211,
-	gfn__X_udd_212,
-	gfn__X_udd_213,
-	gfn__X_udd_222,
-	gfn__X_udd_223,
-	gfn__X_udd_233,
+	gfn__global_x = min_gfn,	// n.b. no macro for this in "cg.hh"
+	gfn__global_y,			// (ditto)
+	gfn__global_z,			// (ditto)
 	gfn__g_dd_11,
 	gfn__g_dd_12,
 	gfn__g_dd_13,
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 0cb1f7e..6ea0a8d 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -3,8 +3,9 @@
 //
 // <<<prototypes for functions local to this file>>>
 // horizon_function - top-level driver
-/// xyz_posns_and_deriv_coeffs - init xyz posns of grid points, xyz deriv coeffs
+/// setup_xyz_posns - setup global xyz posns of grid points
 /// interpolate_geometry - interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid
+/// compute_H - compute H(h) given earlier setup
 //
 
 #include <stdio.h>
@@ -44,7 +45,7 @@ using jtutil::error_exit;
 //
 
 namespace {
-void xyz_posns_and_deriv_coeffs(patch_system& ps);
+void setup_xyz_posns(patch_system& ps);
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_interpolator_info& ii);
@@ -67,20 +68,22 @@ void compute_H(patch_system& ps);
 //	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
 //	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
 //
-// Outputs (angular gridfns, all on the nominal grid):
-//	## computed directly from h
+// Outputs (temporaries computed at each grid point)
+//	## computed by hand-written code
 //	xx,yy,zz			# xyz positions of grid points
 //	X_ud_*, X_udd_*			# xyz derivative coefficients
-//	## interpolated from 3-D Cactus grid
-//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
-//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
-//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
 //	## computed by Maple-generated code
 //	g_uu_{11,12,13,22,23,33}	# $g^{ij}$
 //	K				# $K$
 //	K_dd_{11,12,13,22,23,33}	# $K^{ij}$
 //	partial_d_ln_sqrt_g_d		# $\partial_i \ln \sqrt{g}$
 //	partial_d_g_uu_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g^{ij}$
+//
+// Outputs (angular gridfns, all on the nominal grid):
+//	## interpolated from 3-D Cactus grid
+//	g_dd_{11,12,13,22,23,33}			# $g_{ij}$
+//	K_dd_{11,12,13,22,23,33}			# $K_{ij}$
+//	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
 //	H				# $H = H(h)$
 //
 void horizon_function(patch_system& ps,
@@ -92,8 +95,8 @@ CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 // fill in values of ghosted gridfns in ghost zones
 ps.synchronize_ghost_zones(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h);
 
-// compute xyz positions of grid points and xyz derivative coeffs
-xyz_posns_and_deriv_coeffs(ps);
+// set up xyz positions of grid points
+setup_xyz_posns(ps);
 
 // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
 interpolate_geometry(ps, cgi, ii);
@@ -106,13 +109,11 @@ compute_H(ps);
 //******************************************************************************
 
 //
-// This function computes
-// - the xyz positions of the grid points (gridfns xx, yy, and zz)
-// - the xyz derivative coefficients (gridfns X_ud and X_udd).
-// on the nominal grid.
+// This function sets up the global xyz positions of the grid points
+// in the gridfns global_[xyz].  These will be used by interplate_geometry().
 //
 namespace {
-void xyz_posns_and_deriv_coeffs(patch_system& ps)
+void setup_xyz_posns(patch_system& ps)
 {
 CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 
@@ -126,8 +127,8 @@ CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
-		const fp r
-			= p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma);
+		const fp r = p.ghosted_gridfn(ghosted_gfns::gfn__h,
+					      irho,isigma);
 		const fp rho = p.rho_of_irho(irho);
 		const fp sigma = p.sigma_of_isigma(isigma);
 		fp local_x, local_y, local_z;
@@ -136,50 +137,9 @@ CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 		const fp global_y = ps.global_y_of_local_y(local_y);
 		const fp global_z = ps.global_z_of_local_z(local_z);
 
-		// global xyz positions of grid points
-		p.gridfn(nominal_gfns::gfn__xx, irho,isigma) = global_x;
-		p.gridfn(nominal_gfns::gfn__yy, irho,isigma) = global_y;
-		p.gridfn(nominal_gfns::gfn__zz, irho,isigma) = global_z;
-
-		// 1st derivative coefficient gridfns X_ud
-		p.gridfn(nominal_gfns::gfn__X_ud_11, irho,isigma)
-			= p.partial_rho_wrt_x(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_ud_12, irho,isigma)
-			= p.partial_rho_wrt_y(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_ud_13, irho,isigma)
-			= p.partial_rho_wrt_z(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_ud_21, irho,isigma)
-			= p.partial_sigma_wrt_x(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_ud_22, irho,isigma)
-			= p.partial_sigma_wrt_y(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_ud_23, irho,isigma)
-			= p.partial_sigma_wrt_z(local_x,local_y,local_z);
-
-		// 2nd derivative coefficient gridfns X_udd
-		p.gridfn(nominal_gfns::gfn__X_udd_111, irho,isigma)
-			= p.partial2_rho_wrt_xx(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_112, irho,isigma)
-			= p.partial2_rho_wrt_xy(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_113, irho,isigma)
-			= p.partial2_rho_wrt_xz(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_122, irho,isigma)
-			= p.partial2_rho_wrt_yy(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_123, irho,isigma)
-			= p.partial2_rho_wrt_yz(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_133, irho,isigma)
-			= p.partial2_rho_wrt_zz(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_211, irho,isigma)
-			= p.partial2_sigma_wrt_xx(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_212, irho,isigma)
-			= p.partial2_sigma_wrt_xy(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_213, irho,isigma)
-			= p.partial2_sigma_wrt_xz(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_222, irho,isigma)
-			= p.partial2_sigma_wrt_yy(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_223, irho,isigma)
-			= p.partial2_sigma_wrt_yz(local_x,local_y,local_z);
-		p.gridfn(nominal_gfns::gfn__X_udd_233, irho,isigma)
-			= p.partial2_sigma_wrt_zz(local_x,local_y,local_z);
+		p.gridfn(nominal_gfns::gfn__global_x, irho,isigma) = global_x;
+		p.gridfn(nominal_gfns::gfn__global_y, irho,isigma) = global_y;
+		p.gridfn(nominal_gfns::gfn__global_z, irho,isigma) = global_z;
 		}
 		}
 	}
@@ -231,9 +191,9 @@ const int N_interp_points = ps.N_grid_points();
 const int interp_coords_type_code = CCTK_VARIABLE_REAL;
 const void* const interp_coords[N_GRID_DIMS]
   = {
-    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__xx)),
-    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__yy)),
-    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__zz)),
+    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_x)),
+    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_y)),
+    static_cast<const void*>(ps.gridfn_data(nominal_gfns::gfn__global_z)),
     };
 
 const int N_INPUT_ARRAYS = 12;
@@ -410,9 +370,12 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
 
 	assert(out_of_range_pt >= 0);
 	assert(out_of_range_pt < ps.N_grid_points());
-	const double x = ps.gridfn_data(nominal_gfns::gfn__xx)[out_of_range_pt];
-	const double y = ps.gridfn_data(nominal_gfns::gfn__yy)[out_of_range_pt];
-	const double z = ps.gridfn_data(nominal_gfns::gfn__zz)[out_of_range_pt];
+	const double global_x = ps.gridfn_data(nominal_gfns::gfn__global_x)
+				[out_of_range_pt];
+	const double global_y = ps.gridfn_data(nominal_gfns::gfn__global_y)
+				[out_of_range_pt];
+	const double global_z = ps.gridfn_data(nominal_gfns::gfn__global_z)
+				[out_of_range_pt];
 
 	assert(out_of_range_axis >= 0);
 	assert(out_of_range_axis < N_GRID_DIMS);
@@ -426,7 +389,7 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
 "        the point (%g,%g,%g) on the trial horizon surface\n"
 "        is outside the grid in the %c%c direction!\n"
 ,
-		   x, y, z,
+		   global_x, global_y, global_z,
 		   end, axis);					/*NOTREACHED*/
 	}
 if (status < 0)
@@ -459,9 +422,44 @@ CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
+		//
+		// compute the X_ud and X_udd derivative coefficients
+		// ... n.b. this uses the *local* (x,y,z) coordinates
+		//
+		const fp r = p.ghosted_gridfn(ghosted_gfns::gfn__h,
+					      irho,isigma);
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+		fp xx, yy, zz;
+		p.xyz_of_r_rho_sigma(r,rho,sigma, xx, yy, zz);
+
+		// 1st derivative coefficients X_ud
+		const fp X_ud_11 = p.partial_rho_wrt_x(xx, yy, zz);
+		const fp X_ud_12 = p.partial_rho_wrt_y(xx, yy, zz);
+		const fp X_ud_13 = p.partial_rho_wrt_z(xx, yy, zz);
+		const fp X_ud_21 = p.partial_sigma_wrt_x(xx, yy, zz);
+		const fp X_ud_22 = p.partial_sigma_wrt_y(xx, yy, zz);
+		const fp X_ud_23 = p.partial_sigma_wrt_z(xx, yy, zz);
+
+		// 2nd derivative coefficient gridfns X_udd
+		const fp X_udd_111 = p.partial2_rho_wrt_xx(xx, yy, zz);
+		const fp X_udd_112 = p.partial2_rho_wrt_xy(xx, yy, zz);
+		const fp X_udd_113 = p.partial2_rho_wrt_xz(xx, yy, zz);
+		const fp X_udd_122 = p.partial2_rho_wrt_yy(xx, yy, zz);
+		const fp X_udd_123 = p.partial2_rho_wrt_yz(xx, yy, zz);
+		const fp X_udd_133 = p.partial2_rho_wrt_zz(xx, yy, zz);
+		const fp X_udd_211 = p.partial2_sigma_wrt_xx(xx, yy, zz);
+		const fp X_udd_212 = p.partial2_sigma_wrt_xy(xx, yy, zz);
+		const fp X_udd_213 = p.partial2_sigma_wrt_xz(xx, yy, zz);
+		const fp X_udd_222 = p.partial2_sigma_wrt_yy(xx, yy, zz);
+		const fp X_udd_223 = p.partial2_sigma_wrt_yz(xx, yy, zz);
+		const fp X_udd_233 = p.partial2_sigma_wrt_zz(xx, yy, zz);
+
+		//
 		// "call" the Maple-generated code
-		// ... each file has a separate set of temporary variables,
+		// ... each cg/*.c file has a separate set of temp variables,
 		//     and so must be inside its own set of { } braces
+		//
 		#include "cg.hh"
 		  {
 		#include "cg/inverse_metric.c"
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 81a1c48..1b630fc 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -99,14 +99,14 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set;
     x_xyz[1] := xx;
     x_xyz[2] := yy;
     x_xyz[3] := zz;
-    x_xyz_list := [x_xyz[1], x_xyz[2], x_xyz[3]];
-    x_xyz_set := {op(x_xyz_list)};
+    x_xyz_list := [xx, yy, zz];
+    x_xyz_set := {yy, xx, zz};
     r__fnd := ssqrt(xx^2 + yy^2 + zz^2);
     y_rs := array(1 .. N_ang);
     y_rs[1] := rho;
     y_rs[2] := sigma;
-    y_rs_list := [y_rs[1], y_rs[2]];
-    y_rs_set := {op(y_rs_list)};
+    y_rs_list := [rho, sigma];
+    y_rs_set := {sigma, rho};
     xy_all_list := [op(x_xyz_list), op(y_rs_list)];
     xy_all_set := {op(xy_all_list)};
     NULL
@@ -567,9 +567,9 @@ local nn;
 end proc
 
 > 
-> read "setup_gfas.mm";
+> read "setup_gr_gfas.mm";
 setup_gr_gfas := proc()
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -599,7 +599,7 @@ end proc
 
 Diff/gridfn := proc(operand)
 local var_list, posn;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -623,7 +623,7 @@ auxiliary := proc(cg_flag::boolean)
 end proc
 
 inverse_metric := proc(cg_flag::boolean)
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -641,7 +641,7 @@ end proc
 
 extrinsic_curvature_trace_raise := proc(cg_flag::boolean)
 local i, j, m, n;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -673,7 +673,7 @@ end proc
 
 inverse_metric_gradient := proc(cg_flag::boolean)
 local k, i, j, m, n;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -699,7 +699,7 @@ end proc
 
 metric_det_gradient := proc(cg_flag::boolean)
 local k, i, j;
-global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
+global delta, N, N_ang, x, y, z, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
 rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
 fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
 g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
@@ -826,7 +826,7 @@ codegen2(g_uu) --> "cg/inverse_metric.c"
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=1000044, alloc=917336, time=0.05
+bytes used=1000180, alloc=917336, time=0.08
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
@@ -839,7 +839,7 @@ codegen2([K, K_uu]) --> "cg/extrinsic_curvature_trace_raise.c"
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=2000244, alloc=1441528, time=0.13
+bytes used=2000764, alloc=1441528, time=0.17
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -847,39 +847,39 @@ bytes used=2000244, alloc=1441528, time=0.13
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=3000564, alloc=1572576, time=0.21
+bytes used=3001128, alloc=1572576, time=0.22
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
 > metric_derivs(true);
 inverse_metric_gradient...
-bytes used=4001040, alloc=1638100, time=0.29
+bytes used=4001740, alloc=1638100, time=0.28
 codegen2(partial_d_g_uu) --> "cg/inverse_metric_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=5001252, alloc=1638100, time=0.37
+bytes used=5001904, alloc=1638100, time=0.35
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
       --> `codegen2/fix_Diff`
-bytes used=6001552, alloc=1703624, time=0.44
+bytes used=6002084, alloc=1703624, time=0.42
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=7001892, alloc=1703624, time=0.50
+bytes used=7002504, alloc=1703624, time=0.49
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
-bytes used=8002240, alloc=1703624, time=0.56
+bytes used=8002836, alloc=1703624, time=0.53
 metric_det_gradient...
 codegen2(partial_d_ln_sqrt_g) --> "cg/metric_det_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=9002448, alloc=1703624, time=0.66
+bytes used=9003068, alloc=1703624, time=0.65
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -894,40 +894,40 @@ codegen/C/expression:   Unknown function:   RATIONAL   will be left as is.
 > horizon(true);
 non_unit_normal...
 non_unit_normal_deriv...
-bytes used=10002720, alloc=1769148, time=0.73
+bytes used=10003284, alloc=1769148, time=0.74
 horizon_function...
-bytes used=11003044, alloc=1834672, time=0.80
+bytes used=11003584, alloc=1834672, time=0.81
 codegen2(H) --> "cg/horizon_function.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=12003924, alloc=1834672, time=0.87
-bytes used=13004132, alloc=2031244, time=0.97
-bytes used=14004312, alloc=2031244, time=1.04
+bytes used=12004160, alloc=1834672, time=0.88
+bytes used=13004312, alloc=2031244, time=0.97
+bytes used=14004700, alloc=2031244, time=1.05
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=15004792, alloc=2096768, time=1.08
+bytes used=15005124, alloc=2096768, time=1.13
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
-bytes used=16004952, alloc=2096768, time=1.15
+bytes used=16005504, alloc=2096768, time=1.20
       --> `codegen2/unindex`
-bytes used=17005192, alloc=2096768, time=1.21
+bytes used=17005672, alloc=2096768, time=1.27
    convert p/q --> RATIONAL(p/q)
-bytes used=18005528, alloc=2096768, time=1.27
+bytes used=18006112, alloc=2096768, time=1.33
       --> `codegen2/fix_rationals`
    writing C code
 codegen/C/expression:   Unknown function:   PARTIAL_RHO   will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_SIGMA   will be left as is.
-bytes used=19005708, alloc=2096768, time=1.32
+bytes used=19006376, alloc=2096768, time=1.40
 codegen/C/expression:   Unknown function:   PARTIAL_RHO_RHO   
 will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_RHO_SIGMA   
 will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_SIGMA_SIGMA   
 will be left as is.
-bytes used=20005860, alloc=2096768, time=1.41
+bytes used=20006612, alloc=2096768, time=1.52
 > quit
-bytes used=20241868, alloc=2096768, time=1.45
+bytes used=20238080, alloc=2096768, time=1.54