From 74ef815aab12c0f9602b22b4f0542ca7c3a962be Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 22 Apr 2002 15:10:32 +0000 Subject: change names alpha --> xcos, beta --> ycos, gamma --> zcos because maple predefines gamma = Euler's constant git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@575 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- src/gr/ellipsoid.c | 20 ++--- src/gr/ellipsoid.maple | 7 +- src/gr/ellipsoid.out | 196 ++++++++++++++++++++++++------------------------- 3 files changed, 111 insertions(+), 112 deletions(-) diff --git a/src/gr/ellipsoid.c b/src/gr/ellipsoid.c index b600b76..49a6457 100644 --- a/src/gr/ellipsoid.c +++ b/src/gr/ellipsoid.c @@ -4,25 +4,25 @@ fp t55; t1 = a*a; t2 = b*b; t3 = t1*t2; - t5 = t3*0.5772156649015329*CW; + t5 = t3*zcos*CW; t6 = c*c; t7 = t1*t6; - t9 = t7*beta*BV; + t9 = t7*ycos*BV; t10 = t2*t6; - t12 = t10*alpha*AU; - t28 = alpha*alpha; + t12 = t10*xcos*AU; + t28 = xcos*xcos; t30 = CW*CW; t33 = BV*BV; t35 = t10*t28; - t36 = beta*beta; + t36 = ycos*ycos; t40 = AU*AU; t42 = t7*t36; - t43 = 0.5772156649015329*0.5772156649015329; + t43 = zcos*zcos; t48 = t3*t43; - t49 = -2.0*t1*0.5772156649015329*CW*beta*BV-2.0*t2*0.5772156649015329*CW* -alpha*AU-2.0*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; + 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/(t48+t42+t35); + t55 = 1/(t35+t42+t48); r_plus = (t5+t9+t12+t52)*t55; r_minus = (t5+t9+t12-t52)*t55; diff --git a/src/gr/ellipsoid.maple b/src/gr/ellipsoid.maple index c9b0ae4..0c485fa 100644 --- a/src/gr/ellipsoid.maple +++ b/src/gr/ellipsoid.maple @@ -6,10 +6,11 @@ # 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) +# 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+alpha*r - A)^2 (V+beta*r - B)^2 (W+gamma*r - C)^2 +# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 # ----------------- + ---------------- + ----------------- = 1 # a^2 b^2 c^2 # @@ -18,7 +19,7 @@ # 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; +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"; diff --git a/src/gr/ellipsoid.out b/src/gr/ellipsoid.out index 27b5f21..2e77483 100644 --- a/src/gr/ellipsoid.out +++ b/src/gr/ellipsoid.out @@ -4,17 +4,18 @@ <____ ____> Waterloo Maple Inc. | Type ? for help. # ellipsoid.maple -- compute equations for offset ellipsoid setup -# $Header$ +# $Header: /numrelcvs/ThornburgCVS/AHFinderDirect/src/gr/ellipsoid.maple,v 1.1 2002/04/22 11:02:27 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) -# i.e. a point has coordinates (U+alpha*r, V+beta*r, W+gamma*r) +# 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+alpha*r - A)^2 (V+beta*r - B)^2 (W+gamma*r - C)^2 +# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 # ----------------- + ---------------- + ----------------- = 1 # a^2 b^2 c^2 # @@ -23,12 +24,12 @@ # 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 +> 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"; @@ -384,136 +385,132 @@ 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 ( +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 gamma CW c beta BV + 2 a b gamma CW c alpha AU + 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 beta BV b alpha AU - b c alpha a CW - b c alpha a BV + 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 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 + + 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 beta b - a b gamma c BV - a b gamma c 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 gamma c ) ) / (a b gamma + a c beta + b c alpha ), - / + 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 - 1/2 (2 a b gamma CW + 2 a c beta BV + 2 b c alpha AU - 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 gamma CW c beta BV + 2 a b gamma CW c alpha AU + 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 beta BV b alpha AU - b c alpha a CW - b c alpha a BV + 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 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 + + 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 beta b - a b gamma c BV - a b gamma c 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 gamma c ) ) / (a b gamma + a c beta + b c alpha )] - / + 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=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 ( +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 gamma CW beta BV - 2 b gamma CW alpha AU - 2 c beta BV alpha AU + 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 - + b alpha CW + c alpha BV - b c alpha + a beta CW + 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 - + c beta AU - a c beta + a gamma BV + b gamma 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 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 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 2 2 - a b gamma CW + a c beta BV + b c alpha AU - (-a b c ( + 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 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 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 - + b alpha CW + c alpha BV - b c alpha + a beta CW + 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 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 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 gamma CW + a c beta BV + b c alpha AU + (-a b c ( + 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 gamma CW beta BV - 2 b gamma CW alpha AU - 2 c beta BV alpha AU + 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 - + b alpha CW + c alpha BV - b c alpha + a beta CW + 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 - + c beta AU - a c beta + a gamma BV + b gamma 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 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 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 2 - r_minus = (a b gamma CW + a c beta BV + b c alpha AU - (-a b c ( + 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 a gamma CW beta BV - 2 b gamma CW alpha AU - 2 c beta BV alpha AU + 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 - + b alpha CW + c alpha BV - b c alpha + a beta 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 2 2 2 - + c beta AU - a c beta + a gamma BV + b gamma AU + 2 2 2 2 2 2 2 2 2 1/2 / + + a zcos BV + b zcos AU - a b zcos )) ) / ( + / - 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 + 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 gamma CW, t6 = c , + 2 2 2 +solnlist := [t1 = a , t2 = b , t3 = t1 t2, t5 = t3 zcos CW, t6 = c , t7 = t1 t6, - 2 - t7 = t1 t6, t9 = t7 beta BV, t10 = t2 t6, t12 = t10 alpha AU, t28 = alpha , + 2 2 + t9 = t7 ycos BV, t10 = t2 t6, t12 = t10 xcos AU, t28 = xcos , t30 = CW , - 2 2 2 2 - t30 = CW , t33 = BV , t35 = t10 t28, t36 = beta , t40 = AU , t42 = t7 t36, + 2 2 2 2 + t33 = BV , t35 = t10 t28, t36 = ycos , t40 = AU , t42 = t7 t36, t43 = zcos , - 2 - t43 = gamma , t48 = t3 t43, t49 = -2 t1 gamma CW beta BV + t48 = t3 t43, t49 = -2 t1 zcos CW ycos BV - 2 t2 zcos CW xcos AU - - 2 t2 gamma CW alpha AU - 2 t6 beta BV alpha AU + t2 t28 t30 + t6 t28 t33 + - 2 t6 ycos BV xcos AU + t2 t28 t30 + t6 t28 t33 - t35 + t1 t36 t30 - - t35 + t1 t36 t30 + t6 t36 t40 - t42 + t1 t43 t33 + t2 t43 t40 - t48, + 1/2 + + t6 t36 t40 - t42 + t1 t43 t33 + t2 t43 t40 - t48, t52 = (-t3 t6 t49) , - 1/2 1 - t52 = (-t3 t6 t49) , t55 = ---------------, - t48 + t42 + t35 + 1 + t55 = ---------------, r_plus = (t5 + t9 + t12 + t52) t55, + t35 + t42 + t48 - r_plus = (t5 + t9 + t12 + t52) t55, r_minus = (t5 + t9 + t12 - t52) t55] + r_minus = (t5 + t9 + t12 - t52) t55] > > ftruncate("ellipsoid.c"); @@ -522,5 +519,6 @@ solnlist := [t1 = a , t2 = b , t3 = t1 t2, t5 = t3 gamma CW, t6 = 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 > quit -bytes used=4749864, alloc=1834672, time=0.37 +bytes used=4027884, alloc=1769148, time=0.27 -- cgit v1.2.3