Maple code to compute r(angular coords) for points on offset ellipsoid

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

1 files changed, 526 insertions, 0 deletions

diff --git a/src/gr/ellipsoid.out b/src/gr/ellipsoid.out
new file mode 100644
index 0000000..27b5f21
--- /dev/null
+++ b/src/gr/ellipsoid.out
@@ -0,0 +1,526 @@
+    |\^/|     Maple 7 (IBM INTEL LINUX)
+._|\|   |/|_. Copyright (c) 2001 by Waterloo Maple Inc.
+ \  MAPLE  /  All rights reserved. Maple is a registered trademark of
+ <____ ____>  Waterloo Maple Inc.
+      |       Type ? for help.
+# ellipsoid.maple -- compute equations for offset ellipsoid setup
+# $Header$
+> 
+#
+# ellipsoid has center (A,B,C), radius (a,b,c)
+# angular coordinate system has center (U,V,W)
+#
+# direction cosines wrt angular coordinate center are (alpha,beta,gamma)
+# i.e. a point has coordinates (U+alpha*r, V+beta*r, W+gamma*r)
+#
+# then the equation of the ellipsoid is
+#	(U+alpha*r - A)^2     (V+beta*r - B)^2     (W+gamma*r - C)^2
+#	-----------------  +  ----------------  +  -----------------  =  1
+#	        a^2                  b^2                   c^2
+#
+# to solve this, we introduce intermediate variables
+#	AU = A - U
+#	BV = B - V
+#	CW = C - W
+#
+> eqn := (alpha*r - AU)^2/a^2 + (beta*r - BV)^2/b^2 + (gamma*r - CW)^2/c^2 = 1;
+                              2                2                 2
+                (alpha r - AU)    (beta r - BV)    (gamma r - CW)
+         eqn := --------------- + -------------- + --------------- = 1
+                       2                 2                2
+                      a                 b                c
+
+> 
+> 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=1000108, alloc=917336, time=0.09
+bytes used=2000516, alloc=1376004, time=0.15
+         2  2               2  2              2  2
+[1/2 (2 a  b  gamma CW + 2 a  c  beta BV + 2 b  c  alpha AU + 2 (
+
+       4  2           2              2  4           2
+    2 a  b  gamma CW c  beta BV + 2 a  b  gamma CW c  alpha AU
+
+          2  4          2             4  2      2  2   2    2  4      2  2   2
+     + 2 a  c  beta BV b  alpha AU - b  c  alpha  a  CW  - b  c  alpha  a  BV
+
+        4  4      2  2    4  2     2  2   2    2  4     2  2   2
+     + b  c  alpha  a  - a  c  beta  b  CW  - a  c  beta  b  AU
+
+        4  4     2  2    4  2      2  2   2    2  4      2  2   2
+     + a  c  beta  b  - a  b  gamma  c  BV  - a  b  gamma  c  AU
+
+        4  4      2  2 1/2    /   2  2      2    2  2     2    2  2      2
+     + a  b  gamma  c )   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha ),
+                            /
+
+            2  2               2  2              2  2
+    1/2 (2 a  b  gamma CW + 2 a  c  beta BV + 2 b  c  alpha AU - 2 (
+
+       4  2           2              2  4           2
+    2 a  b  gamma CW c  beta BV + 2 a  b  gamma CW c  alpha AU
+
+          2  4          2             4  2      2  2   2    2  4      2  2   2
+     + 2 a  c  beta BV b  alpha AU - b  c  alpha  a  CW  - b  c  alpha  a  BV
+
+        4  4      2  2    4  2     2  2   2    2  4     2  2   2
+     + b  c  alpha  a  - a  c  beta  b  CW  - a  c  beta  b  AU
+
+        4  4     2  2    4  2      2  2   2    2  4      2  2   2
+     + a  c  beta  b  - a  b  gamma  c  BV  - a  b  gamma  c  AU
+
+        4  4      2  2 1/2    /   2  2      2    2  2     2    2  2      2
+     + a  b  gamma  c )   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha )]
+                            /
+
+> map(simplify, %);
+bytes used=3007136, alloc=1769148, time=0.22
+bytes used=4007324, alloc=1834672, time=0.29
+   2  2             2  2            2  2               2  2  2
+[(a  b  gamma CW + a  c  beta BV + b  c  alpha AU + (-a  b  c  (
+
+        2                       2                        2
+    -2 a  gamma CW beta BV - 2 b  gamma CW alpha AU - 2 c  beta BV alpha AU
+
+        2      2   2    2      2   2    2  2      2    2     2   2
+     + b  alpha  CW  + c  alpha  BV  - b  c  alpha  + a  beta  CW
+
+        2     2   2    2  2     2    2      2   2    2      2   2
+     + c  beta  AU  - a  c  beta  + a  gamma  BV  + b  gamma  AU
+
+        2  2      2  1/2    /   2  2      2    2  2     2    2  2      2
+     - a  b  gamma ))   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha ), (
+                          /
+
+     2  2             2  2            2  2               2  2  2
+    a  b  gamma CW + a  c  beta BV + b  c  alpha AU - (-a  b  c  (
+
+        2                       2                        2
+    -2 a  gamma CW beta BV - 2 b  gamma CW alpha AU - 2 c  beta BV alpha AU
+
+        2      2   2    2      2   2    2  2      2    2     2   2
+     + b  alpha  CW  + c  alpha  BV  - b  c  alpha  + a  beta  CW
+
+        2     2   2    2  2     2    2      2   2    2      2   2
+     + c  beta  AU  - a  c  beta  + a  gamma  BV  + b  gamma  AU
+
+        2  2      2  1/2    /   2  2      2    2  2     2    2  2      2
+     - a  b  gamma ))   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha )]
+                          /
+
+> [r_plus = %[1], r_minus = %[2]];
+            2  2             2  2            2  2               2  2  2
+[r_plus = (a  b  gamma CW + a  c  beta BV + b  c  alpha AU + (-a  b  c  (
+
+        2                       2                        2
+    -2 a  gamma CW beta BV - 2 b  gamma CW alpha AU - 2 c  beta BV alpha AU
+
+        2      2   2    2      2   2    2  2      2    2     2   2
+     + b  alpha  CW  + c  alpha  BV  - b  c  alpha  + a  beta  CW
+
+        2     2   2    2  2     2    2      2   2    2      2   2
+     + c  beta  AU  - a  c  beta  + a  gamma  BV  + b  gamma  AU
+
+        2  2      2  1/2    /   2  2      2    2  2     2    2  2      2
+     - a  b  gamma ))   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha ),
+                          /
+
+                2  2             2  2            2  2               2  2  2
+    r_minus = (a  b  gamma CW + a  c  beta BV + b  c  alpha AU - (-a  b  c  (
+
+        2                       2                        2
+    -2 a  gamma CW beta BV - 2 b  gamma CW alpha AU - 2 c  beta BV alpha AU
+
+        2      2   2    2      2   2    2  2      2    2     2   2
+     + b  alpha  CW  + c  alpha  BV  - b  c  alpha  + a  beta  CW
+
+        2     2   2    2  2     2    2      2   2    2      2   2
+     + c  beta  AU  - a  c  beta  + a  gamma  BV  + b  gamma  AU
+
+        2  2      2  1/2    /   2  2      2    2  2     2    2  2      2
+     - a  b  gamma ))   )  /  (a  b  gamma  + a  c  beta  + b  c  alpha )]
+                          /
+
+> solnlist := [codegen[optimize](%)];
+                   2        2                                      2
+solnlist := [t1 = a , t2 = b , t3 = t1 t2, t5 = t3 gamma CW, t6 = c ,
+
+                                                                             2
+    t7 = t1 t6, t9 = t7 beta BV, t10 = t2 t6, t12 = t10 alpha AU, t28 = alpha ,
+
+            2          2                           2          2
+    t30 = CW , t33 = BV , t35 = t10 t28, t36 = beta , t40 = AU , t42 = t7 t36,
+
+               2
+    t43 = gamma , t48 = t3 t43, t49 = -2 t1 gamma CW beta BV
+
+     - 2 t2 gamma CW alpha AU - 2 t6 beta BV alpha AU + t2 t28 t30 + t6 t28 t33
+
+     - t35 + t1 t36 t30 + t6 t36 t40 - t42 + t1 t43 t33 + t2 t43 t40 - t48,
+
+                      1/2               1
+    t52 = (-t3 t6 t49)   , t55 = ---------------,
+                                 t48 + t42 + t35
+
+    r_plus = (t5 + t9 + t12 + t52) t55, 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");
+> quit
+bytes used=4749864, alloc=1834672, time=0.37