* expand NaN checks in gridfns to "is finite" checks

  (also catchs +/-infinity)
* move final stage of calculation of LHS function H(h)
  out of Maple-generated code (in src/gr.cg/) into hand-written code
  in src/gr/horizon_function.c (comput_H()) so we can check that HD > 0
  before taking its sqrt.
	[HD <= 0 probably means the inerpolated g_ij
	isn't positive definite]
* also systematize CCTK_VWarn() printing for various error conditions
  in horizon_function.c


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

1 files changed, 168 insertions, 170 deletions

diff --git a/src/gr.cg/horizon_function.c b/src/gr.cg/horizon_function.c
index c31778c..7e6e717 100644
--- a/src/gr.cg/horizon_function.c
+++ b/src/gr.cg/horizon_function.c
@@ -1,23 +1,22 @@
 /*
- * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, h, HA, HB, HC, HD, K}
- * outputs = {H, HA, HB, HC, HD}
- * cost = 136*assignments+7*divisions+7*functions+402*multiplications+176*additions
+ * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, h}
+ * outputs = {HA, HB, HC, HD}
+ * cost = 134*assignments+3*divisions+5*functions+401*multiplications+173*additions
  */
 fp t1, t2, t4, t5, t7, t8, t10, t11, t12, t14;
-fp t16, t18, t19, t21, t23, t25, t28, t29, t30, t32;
-fp t33, t35, t37, t43, t44, t47, t50, t51, t53, t55;
-fp t56, t57, t58, t61, t62, t63, t64, t65, t67, t69;
-fp t72, t73, t76, t78, t85, t90, t92, t94, t95, t96;
-fp t99, t100, t101, t102, t105, t106, t107, t111, t114, t115;
-fp t117, t118, t133, t134, t137, t138, t141, t142, t145, t146;
-fp t149, t150, t158, t163, t165, t169, t170, t173, t174, t179;
-fp t194, t198, t199, t204, t220, t227, t234, t237, t239, t242;
-fp t247, t250, t254, t261, t262, t265, t267, t270, t273, t276;
-fp t279, t292, t295, t296, t299, t306, t307, t312, t315, t318;
-fp t320, t323, t333, t342, t343, t346, t349, t350, t361, t362;
-fp t368, t378, t394, t413, t431, t440, t444, t447, t450, t465;
-fp t490;
-      t1 = g_uu_13;
+fp t16, t18, t19, t20, t23, t24, t25, t26, t29, t30;
+fp t31, t33, t35, t37, t38, t40, t41, t42, t43, t45;
+fp t46, t47, t48, t50, t51, t52, t53, t55, t58, t59;
+fp t61, t64, t67, t68, t70, t72, t73, t75, t77, t83;
+fp t84, t87, t90, t91, t93, t95, t112, t130, t134, t137;
+fp t139, t143, t147, t150, t152, t159, t164, t166, t168, t176;
+fp t181, t183, t199, t203, t206, t209, t213, t216, t217, t219;
+fp t228, t229, t232, t235, t236, t237, t240, t243, t247, t251;
+fp t254, t259, t260, t262, t263, t269, t272, t273, t276, t279;
+fp t282, t283, t290, t293, t296, t299, t302, t305, t312, t315;
+fp t318, t324, t325, t329, t330, t338, t358, t381, t392, t395;
+fp t396, t397, t408, t411, t431, t440, t444, t447, t450, t465;
+      t1 = g_uu_33;
       t2 = 1/r;
       t4 = X_ud_13;
       t5 = PARTIAL_RHO(h);
@@ -25,160 +24,159 @@ fp t490;
       t8 = PARTIAL_SIGMA(h);
       t10 = zz*t2-t4*t5-t7*t8;
       t11 = t1*t10;
-      t12 = partial_d_g_uu_112;
-      t14 = X_ud_12;
-      t16 = X_ud_22;
-      t18 = yy*t2-t14*t5-t16*t8;
-      t19 = t12*t18;
-      t21 = X_ud_11;
-      t23 = X_ud_21;
-      t25 = xx*t2-t21*t5-t23*t8;
-      t28 = g_uu_23;
-      t29 = t28*t28;
-      t30 = t10*t10;
-      t32 = xx*xx;
-      t33 = zz*zz;
-      t35 = r*r;
-      t37 = 1/t35/r;
-      t43 = t14*t14;
-      t44 = PARTIAL_RHO_RHO(h);
-      t47 = PARTIAL_RHO_SIGMA(h);
-      t50 = t16*t16;
-      t51 = PARTIAL_SIGMA_SIGMA(h);
-      t53 = (t32+t33)*t37-X_udd_122*t5-X_udd_222*t8-t43*t44-2.0*t16*t14*t47-t50
-*t51;
-      t55 = g_uu_33;
-      t56 = t55*t10;
-      t57 = partial_d_g_uu_312;
-      t58 = t57*t18;
-      t61 = g_uu_22;
-      t62 = t18*t18;
-      t63 = t61*t62;
-      t64 = partial_d_g_uu_212;
-      t65 = t64*t25;
-      t67 = t30*t10;
-      t69 = partial_d_g_uu_233;
-      t72 = t28*t10;
-      t73 = t64*t18;
-      t76 = t55*t55;
-      t78 = yy*yy;
-      t85 = t4*t4;
-      t90 = t7*t7;
-      t92 = (t32+t78)*t37-X_udd_133*t5-X_udd_233*t8-t85*t44-2.0*t7*t4*t47-t90*
-t51;
-      t94 = t61*t18;
-      t95 = partial_d_g_uu_213;
-      t96 = t95*t10;
-      t99 = g_uu_12;
-      t100 = t99*t18;
-      t101 = partial_d_g_uu_113;
-      t102 = t101*t10;
-      t105 = t28*t18;
-      t106 = partial_d_g_uu_313;
-      t107 = t106*t10;
-      t111 = partial_d_g_uu_133;
-      t114 = t25*t25;
-      t115 = t1*t114;
-      t117 = -t11*t19*t25-t29*t30*t53-t56*t58*t25-t63*t65+RATIONAL(-1.0,2.0)*
-t28*t67*t69-t72*t73*t25-t76*t30*t92-t94*t96*t25-t100*t102*t25-t105*t107*t25+
-RATIONAL(-1.0,2.0)*t1*t67*t111-t115*t107;
-      t118 = t55*t30;
-      t133 = -yy*zz*t37-X_udd_123*t5-X_udd_223*t8-t14*t4*t44-t16*t4*t47-t14*t7*
-t47-t16*t7*t51;
-      t134 = t28*t133;
-      t137 = partial_d_g_uu_111;
-      t138 = t137*t114;
-      t141 = partial_d_g_uu_311;
-      t142 = t141*t114;
-      t145 = partial_d_g_uu_211;
-      t146 = t145*t114;
-      t149 = g_uu_11;
-      t150 = t149*t149;
-      t158 = t21*t21;
-      t163 = t23*t23;
-      t165 = (t78+t33)*t37-X_udd_111*t5-X_udd_211*t8-t158*t44-2.0*t23*t21*t47-
-t163*t51;
-      t169 = t99*t25;
-      t170 = t69*t30;
-      t173 = t149*t25;
-      t174 = t111*t30;
-      t179 = t99*t114;
-      t194 = -xx*yy*t37-X_udd_112*t5-X_udd_212*t8-t21*t14*t44-t23*t14*t47-t21*
-t16*t47-t23*t16*t51;
-      t198 = partial_d_g_uu_333;
-      t199 = t198*t30;
-      t204 = -2.0*t118*t134+RATIONAL(-1.0,2.0)*t11*t138+RATIONAL(-1.0,2.0)*t105
-*t142+RATIONAL(-1.0,2.0)*t94*t146-t150*t114*t165+RATIONAL(-1.0,2.0)*t100*t138+
-RATIONAL(-1.0,2.0)*t169*t170+RATIONAL(-1.0,2.0)*t173*t174+RATIONAL(-1.0,2.0)*
-t100*t174-2.0*t179*t149*t194+RATIONAL(-1.0,2.0)*t105*t199+RATIONAL(-1.0,2.0)*
-t94*t170;
-      t220 = -xx*zz*t37-X_udd_113*t5-X_udd_213*t8-t21*t4*t44-t23*t4*t47-t21*t7*
-t47-t23*t7*t51;
-      t227 = t99*t194;
-      t234 = t99*t99;
-      t237 = t149*t114;
-      t239 = t28*t30;
-      t242 = t1*t1;
-      t247 = t61*t61;
-      t250 = t28*t62;
-      t254 = -2.0*t115*t149*t220-2.0*t115*t99*t133-2.0*t63*t227+RATIONAL(-1.0,
-2.0)*t72*t146+RATIONAL(-1.0,2.0)*t56*t142-t234*t114*t53-t237*t19-t239*t95*t25-
-t242*t114*t92-t234*t62*t165-t247*t62*t53-2.0*t250*t99*t220;
-      t261 = partial_d_g_uu_222;
-      t262 = t261*t62;
-      t265 = t106*t25;
-      t267 = t1*t30;
-      t270 = t169*t53;
-      t273 = t173*t220;
-      t276 = t173*t165;
-      t279 = t114*t25;
-      t292 = -2.0*t234*t18*t25*t194-t179*t96-t237*t102+RATIONAL(-1.0,2.0)*t169*
-t262-t118*t265-t267*t101*t25-2.0*t94*t270-2.0*t105*t273-2.0*t100*t276+RATIONAL(
--1.0,2.0)*t149*t279*t137+RATIONAL(-1.0,2.0)*t99*t279*t145+RATIONAL(-1.0,2.0)*
-t55*t67*t198+RATIONAL(-1.0,2.0)*t1*t279*t141;
-      t295 = partial_d_g_uu_122;
-      t296 = t295*t62;
-      t299 = t173*t194;
-      t306 = t1*t25;
-      t307 = t306*t133;
-      t312 = t169*t133;
-      t315 = t306*t92;
-      t318 = t62*t18;
-      t320 = partial_d_g_uu_322;
-      t323 = t320*t62;
-      t333 = RATIONAL(-1.0,2.0)*t173*t296-2.0*t94*t299-2.0*t242*t10*t25*t220
--2.0*t72*t307-t250*t57*t25-2.0*t56*t312-2.0*t56*t315+RATIONAL(-1.0,2.0)*t28*
-t318*t320+RATIONAL(-1.0,2.0)*t306*t323+RATIONAL(-1.0,2.0)*t99*t318*t295-2.0*t72
-*t270+RATIONAL(-1.0,2.0)*t56*t323;
-      t342 = partial_d_g_uu_223;
-      t343 = t342*t10;
-      t346 = t99*t62;
-      t349 = partial_d_g_uu_323;
-      t350 = t349*t10;
-      t361 = partial_d_g_uu_123;
-      t362 = t361*t10;
-      t368 = -t179*t73+RATIONAL(-1.0,2.0)*t61*t318*t261-2.0*t56*t273-2.0*t72*
-t299-t169*t343*t18-t346*t12*t25-t306*t350*t18-2.0*t11*t276-2.0*t94*t307-2.0*
-t105*t315+RATIONAL(-1.0,2.0)*t11*t296-t346*t362-2.0*t29*t10*t18*t133;
-      t378 = t349*t18;
-      t394 = RATIONAL(-1.0,2.0)*t72*t262-t63*t343-t250*t350-t239*t342*t18-t267*
-t361*t18-t118*t378+RATIONAL(-1.0,2.0)*t306*t199-t173*t362*t18-2.0*t105*t312-2.0
-*t100*t306*t220-2.0*t250*t61*t133-t29*t62*t92;
-      t413 = t1*t220;
-      t431 = -2.0*t56*t94*t133-2.0*t56*t105*t92-t115*t58-2.0*t11*t169*t194-t242
-*t30*t165-2.0*t11*t100*t165-2.0*t239*t1*t194-2.0*t118*t413-2.0*t11*t105*t220
--2.0*t72*t94*t53-2.0*t11*t94*t194-2.0*t56*t100*t220-2.0*t72*t100*t194;
-      HA = t117+t204+t254+t292+t333+t368+t394+t431;
-      t440 = t137*t25+t65+t265+t19+t261*t18+t378+t102+t343+t198*t10+t149*t165+
-2.0*t227+2.0*t413;
+      t12 = partial_d_g_uu_311;
+      t14 = X_ud_11;
+      t16 = X_ud_21;
+      t18 = xx*t2-t14*t5-t16*t8;
+      t19 = t18*t18;
+      t20 = t12*t19;
+      t23 = g_uu_23;
+      t24 = t23*t10;
+      t25 = partial_d_g_uu_211;
+      t26 = t25*t19;
+      t29 = g_uu_13;
+      t30 = t29*t19;
+      t31 = partial_d_g_uu_312;
+      t33 = X_ud_12;
+      t35 = X_ud_22;
+      t37 = yy*t2-t33*t5-t35*t8;
+      t38 = t31*t37;
+      t40 = g_uu_12;
+      t41 = t40*t19;
+      t42 = partial_d_g_uu_212;
+      t43 = t42*t37;
+      t45 = g_uu_11;
+      t46 = t45*t19;
+      t47 = partial_d_g_uu_112;
+      t48 = t47*t37;
+      t50 = g_uu_22;
+      t51 = t37*t37;
+      t52 = t50*t51;
+      t53 = t42*t18;
+      t55 = t40*t51;
+      t58 = partial_d_g_uu_113;
+      t59 = t58*t10;
+      t61 = t23*t51;
+      t64 = t29*t10;
+      t67 = partial_d_g_uu_213;
+      t68 = t67*t10;
+      t70 = t45*t45;
+      t72 = yy*yy;
+      t73 = zz*zz;
+      t75 = r*r;
+      t77 = 1/t75/r;
+      t83 = t14*t14;
+      t84 = PARTIAL_RHO_RHO(h);
+      t87 = PARTIAL_RHO_SIGMA(h);
+      t90 = t16*t16;
+      t91 = PARTIAL_SIGMA_SIGMA(h);
+      t93 = (t72+t73)*t77-X_udd_111*t5-X_udd_211*t8-t83*t84-2.0*t16*t14*t87-t90
+*t91;
+      t95 = RATIONAL(-1.0,2.0)*t11*t20+RATIONAL(-1.0,2.0)*t24*t26-t30*t38-t41*
+t43-t46*t48-t52*t53-t55*t47*t18-t46*t59-t61*t31*t18-t64*t48*t18-t41*t68-t70*t19
+*t93;
+      t112 = -xx*yy*t77-X_udd_112*t5-X_udd_212*t8-t14*t33*t84-t16*t33*t87-t14*
+t35*t87-t16*t35*t91;
+      t130 = -xx*zz*t77-X_udd_113*t5-X_udd_213*t8-t14*t4*t84-t16*t4*t87-t14*t7*
+t87-t16*t7*t91;
+      t134 = t40*t37;
+      t137 = t51*t37;
+      t139 = partial_d_g_uu_322;
+      t143 = partial_d_g_uu_122;
+      t147 = partial_d_g_uu_222;
+      t150 = t40*t40;
+      t152 = xx*xx;
+      t159 = t33*t33;
+      t164 = t35*t35;
+      t166 = (t152+t73)*t77-X_udd_122*t5-X_udd_222*t8-t159*t84-2.0*t35*t33*t87-
+t164*t91;
+      t168 = t29*t29;
+      t176 = t4*t4;
+      t181 = t7*t7;
+      t183 = (t152+t72)*t77-X_udd_133*t5-X_udd_233*t8-t176*t84-2.0*t7*t4*t87-
+t181*t91;
+      t199 = -yy*zz*t77-X_udd_123*t5-X_udd_223*t8-t33*t4*t84-t35*t4*t87-t33*t7*
+t87-t35*t7*t91;
+      t203 = t40*t112;
+      t206 = t50*t37;
+      t209 = -t24*t43*t18-2.0*t41*t45*t112-2.0*t30*t45*t130-t134*t59*t18+
+RATIONAL(-1.0,2.0)*t23*t137*t139+RATIONAL(-1.0,2.0)*t40*t137*t143+RATIONAL(-1.0
+,2.0)*t50*t137*t147-t150*t19*t166-t168*t19*t183-2.0*t30*t40*t199-2.0*t52*t203-
+t206*t68*t18;
+      t213 = t50*t50;
+      t216 = t10*t10;
+      t217 = t216*t10;
+      t219 = partial_d_g_uu_333;
+      t228 = partial_d_g_uu_123;
+      t229 = t228*t10;
+      t232 = partial_d_g_uu_233;
+      t235 = t23*t37;
+      t236 = partial_d_g_uu_313;
+      t237 = t236*t10;
+      t240 = t23*t23;
+      t243 = t29*t18;
+      t247 = t243*t183;
+      t251 = partial_d_g_uu_133;
+      t254 = -t150*t51*t93-t213*t51*t166+RATIONAL(-1.0,2.0)*t1*t217*t219-2.0*
+t61*t40*t130-2.0*t61*t50*t199-t55*t229+RATIONAL(-1.0,2.0)*t23*t217*t232-t235*
+t237*t18-t240*t51*t183-2.0*t134*t243*t130-2.0*t235*t247+RATIONAL(-1.0,2.0)*t29*
+t217*t251;
+      t259 = partial_d_g_uu_323;
+      t260 = t259*t10;
+      t262 = partial_d_g_uu_223;
+      t263 = t262*t10;
+      t269 = t243*t199;
+      t272 = t40*t18;
+      t273 = t272*t166;
+      t276 = t272*t199;
+      t279 = t219*t216;
+      t282 = t45*t18;
+      t283 = t251*t216;
+      t290 = t232*t216;
+      t293 = t23*t216;
+      t296 = -2.0*t150*t37*t18*t112-t61*t260-t52*t263-2.0*t168*t10*t18*t130-2.0
+*t206*t269-2.0*t206*t273-2.0*t235*t276+RATIONAL(-1.0,2.0)*t243*t279+RATIONAL(
+-1.0,2.0)*t282*t283+RATIONAL(-1.0,2.0)*t235*t279+RATIONAL(-1.0,2.0)*t134*t283+
+RATIONAL(-1.0,2.0)*t206*t290-t293*t262*t37;
+      t299 = t29*t216;
+      t302 = t147*t51;
+      t305 = t139*t51;
+      t312 = t282*t93;
+      t315 = t282*t130;
+      t318 = t282*t112;
+      t324 = t1*t216;
+      t325 = t236*t18;
+      t329 = -t299*t228*t37+RATIONAL(-1.0,2.0)*t272*t302+RATIONAL(-1.0,2.0)*t11
+*t305+RATIONAL(-1.0,2.0)*t24*t302+RATIONAL(-1.0,2.0)*t272*t290-2.0*t134*t312
+-2.0*t235*t315-2.0*t206*t318-t30*t237-2.0*t24*t269-t324*t325-2.0*t11*t247;
+      t330 = t259*t37;
+      t338 = t143*t51;
+      t358 = -t324*t330-2.0*t11*t315-t272*t263*t37-2.0*t11*t276+RATIONAL(-1.0,
+2.0)*t282*t338-t243*t260*t37-2.0*t64*t312-t282*t229*t37-t299*t58*t18-2.0*t24*
+t273-2.0*t24*t318+RATIONAL(-1.0,2.0)*t64*t338-2.0*t64*t235*t130;
+      t381 = t29*t130;
+      t392 = t23*t199;
+      t395 = -2.0*t240*t10*t37*t199-2.0*t11*t235*t183-2.0*t11*t134*t130-2.0*t64
+*t134*t93-2.0*t24*t206*t166-t11*t38*t18-2.0*t293*t29*t112-2.0*t324*t381-2.0*t11
+*t206*t199-2.0*t64*t206*t112-t168*t216*t93-2.0*t324*t392;
+      t396 = partial_d_g_uu_111;
+      t397 = t396*t19;
+      t408 = t1*t1;
+      t411 = t19*t18;
+      t431 = RATIONAL(-1.0,2.0)*t134*t397+RATIONAL(-1.0,2.0)*t206*t26+RATIONAL(
+-1.0,2.0)*t235*t20+RATIONAL(-1.0,2.0)*t64*t397-t240*t216*t166-t408*t216*t183+
+RATIONAL(-1.0,2.0)*t45*t411*t396+RATIONAL(-1.0,2.0)*t40*t411*t25+RATIONAL(-1.0,
+2.0)*t29*t411*t12-2.0*t24*t134*t112+RATIONAL(-1.0,2.0)*t243*t305-t293*t67*t18
+-2.0*t64*t272*t112;
+      HA = t95+t209+t254+t296+t329+t358+t395+t431;
+      t440 = t396*t18+t53+t325+t48+t147*t37+t330+t59+t263+t219*t10+t45*t93+2.0*
+t203+2.0*t381;
       t444 = partial_d_ln_sqrt_g_1;
       t447 = partial_d_ln_sqrt_g_2;
       t450 = partial_d_ln_sqrt_g_3;
-      t465 = t61*t53+2.0*t134+t55*t92+t149*t444*t25+t99*t447*t25+t1*t450*t25+
-t99*t444*t18+t61*t447*t18+t28*t450*t18+t1*t444*t10+t28*t447*t10+t55*t450*t10;
+      t465 = t50*t166+2.0*t392+t1*t183+t45*t444*t18+t40*t447*t18+t29*t450*t18+
+t40*t444*t37+t50*t447*t37+t23*t450*t37+t29*t444*t10+t23*t447*t10+t1*t450*t10;
       HB = t440+t465;
-      HC = K_uu_11*t114+2.0*K_uu_12*t18*t25+2.0*K_uu_13*t10*t25+K_uu_22*t62+2.0
-*K_uu_23*t10*t18+K_uu_33*t30;
-      HD = t237+2.0*t100*t25+2.0*t11*t25+t63+2.0*t72*t18+t118;
-      t490 = pow(HD,1.0*RATIONAL(1.0,2.0));
-      H = HA/t490/HD+HB/t490+HC/HD-K;
+      HC = K_uu_11*t19+2.0*K_uu_12*t37*t18+2.0*K_uu_13*t10*t18+K_uu_22*t51+2.0*
+K_uu_23*t10*t37+K_uu_33*t216;
+      HD = t46+2.0*t134*t18+2.0*t64*t18+t52+2.0*t24*t37+t324;