directory to hold stuff that's no longer used but that someone might want

to look at again


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

6 files changed, 229 insertions, 0 deletions

diff --git a/archive/README b/archive/README
new file mode 100644
index 0000000..6601932
--- /dev/null
+++ b/archive/README
@@ -0,0 +1,6 @@
+This directory contains stuff that's no longer used, but that I
+thought might be useful for future reference, so I didn't want to
+throw it away.
+
+Compared to the CVS attic, having a separate directory like this
+seems more convenient -- you can 'cat, awk, and grep' more easily.
diff --git a/archive/ellipsoid.c b/archive/ellipsoid.c
new file mode 100644
index 0000000..49a6457
--- /dev/null
+++ b/archive/ellipsoid.c
@@ -0,0 +1,28 @@
+fp t1, t2, t3, t5, t6, t7, t9, t10, t12, t28;
+fp t30, t33, t35, t36, t40, t42, t43, t48, t49, t52;
+fp t55;
+      t1 = a*a;
+      t2 = b*b;
+      t3 = t1*t2;
+      t5 = t3*zcos*CW;
+      t6 = c*c;
+      t7 = t1*t6;
+      t9 = t7*ycos*BV;
+      t10 = t2*t6;
+      t12 = t10*xcos*AU;
+      t28 = xcos*xcos;
+      t30 = CW*CW;
+      t33 = BV*BV;
+      t35 = t10*t28;
+      t36 = ycos*ycos;
+      t40 = AU*AU;
+      t42 = t7*t36;
+      t43 = zcos*zcos;
+      t48 = t3*t43;
+      t49 = -2.0*t1*zcos*CW*ycos*BV-2.0*t2*zcos*CW*xcos*AU-2.0*t6*ycos*BV*xcos*
+AU+t2*t28*t30+t6*t28*t33-t35+t1*t36*t30+t6*t36*t40-t42+t1*t43*t33+t2*t43*t40-
+t48;
+      t52 = sqrt(-t3*t6*t49);
+      t55 = 1/(t35+t42+t48);
+      r_plus = (t5+t9+t12+t52)*t55;
+      r_minus = (t5+t9+t12-t52)*t55;
diff --git a/archive/ellipsoid.maple b/archive/ellipsoid.maple
new file mode 100644
index 0000000..0c485fa
--- /dev/null
+++ b/archive/ellipsoid.maple
@@ -0,0 +1,35 @@
+# 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)
+# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
+# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+#
+# then the equation of the ellipsoid is
+#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*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 := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
+
+read "../maple/util.mm";
+read "../maple/codegen2.mm";
+
+[solve(eqn, r)];
+map(simplify, %);
+[r_plus = %[1], r_minus = %[2]];
+solnlist := [codegen[optimize](%)];
+
+ftruncate("ellipsoid.c");
+print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]),
+		    "fp", "ellipsoid.c");
+codegen[C](solnlist, filename="ellipsoid.c");
diff --git a/archive/ellipsoid.mm b/archive/ellipsoid.mm
new file mode 100644
index 0000000..d58b424
--- /dev/null
+++ b/archive/ellipsoid.mm
@@ -0,0 +1,35 @@
+# ellipsoid.maple -- compute equations for offset ellipsoid setup
+# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/archive/ellipsoid.mm,v 1.1 2002-07-18 17:37:11 jthorn Exp $
+
+#
+# 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)
+# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
+# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+#
+# then the equation of the ellipsoid is
+#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*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 := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
+
+read "../maple/util.mm";
+read "../maple/codegen2.mm";
+
+[solve(eqn, r)];
+map(simplify, %);
+[r_plus = %[1], r_minus = %[2]];
+solnlist := [codegen[optimize](%)];
+
+ftruncate("ellipsoid.c");
+print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]),
+		    "fp", "ellipsoid.c");
+codegen[C](solnlist, filename="ellipsoid.c");
diff --git a/archive/ellipsoid.out b/archive/ellipsoid.out
new file mode 100644
index 0000000..4773c57
--- /dev/null
+++ b/archive/ellipsoid.out
@@ -0,0 +1,38 @@
+    |\^/|     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: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/archive/ellipsoid.out,v 1.1 2002-07-18 17:37:11 jthorn Exp $
+> 
+#
+# 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)
+# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
+# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+#
+# then the equation of the ellipsoid is
+#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*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 := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
+                              2                2                2
+                 (xcos r - AU)    (ycos r - BV)    (zcos r - CW)
+          eqn := -------------- + -------------- + -------------- = 1
+                        2                2                2
+                       a                b                c
+
+> 
+> read "../maple/util.mm";
+Error, unable to read `../maple/util.mm`
+> quit
+bytes used=129844, alloc=196572, time=0.03
diff --git a/archive/setup_ellipsoid_initial_guess.cc b/archive/setup_ellipsoid_initial_guess.cc
new file mode 100644
index 0000000..f8ba391
--- /dev/null
+++ b/archive/setup_ellipsoid_initial_guess.cc
@@ -0,0 +1,87 @@
+//
+// This function sets up an ellipsoid in the gridfn h, using the
+// formulas in "ellipsoid.maple" and the Maple-generated C code in
+// "ellipsoid.c":
+//
+// 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 (xcos,ycos,zcos)
+// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+//
+// then the equation of the ellipsoid is
+//	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*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
+//
+namespace {
+void setup_ellipsoid_initial_guess
+	(patch_system& ps,
+	 fp global_center_x, fp global_center_y, fp global_center_z,
+	 fp radius_x, fp radius_y, fp radius_z)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   h = ellipsoid: global_center=(%g,%g,%g)",
+	   global_center_x, global_center_y, global_center_z);
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "                  radius=(%g,%g,%g)",
+	   radius_x, radius_y, radius_z);
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+		fp xcos, ycos, zcos;
+		p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos);
+
+		// set up variables used by Maple-generated code
+		const fp AU = global_center_x - ps.origin_x();
+		const fp BV = global_center_y - ps.origin_y();
+		const fp CW = global_center_z - ps.origin_z();
+		const fp a = radius_x;
+		const fp b = radius_y;
+		const fp c = radius_z;
+
+		// compute the solutions r_plus and r_minus
+		fp r_plus, r_minus;
+		#include "ellipsoid.c"
+
+		// 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;
+		else if ((r_plus < 0.0) && (r_minus > 0.0))
+		   then r = r_minus;
+		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "\n"
+"   expected exactly one r>0 solution, got 0 or 2!\n"
+"   %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
+"   direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n"
+"   ==> r_plus=%g r_minus=%g\n"
+				   ,
+				   p.name(), irho, isigma,
+				   double(rho), double(sigma),
+				   double(xcos), double(ycos), double(zcos),
+				   double(r_plus), double(r_minus));
+		   						/*NOTREACHED*/
+
+		// r = horizon radius at this grid point
+		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r;
+		}
+		}
+	}
+}
+	  }