/* # PDS_residual: maxima script for evaluating # the residual gravity acceleration effect # induced by spatial topology in the Poincare # dodecahedral space # # Copyright (C) 2009 Boud Roukema # # This program is free software; you can # redistribute it and/or modify it under the # terms of the GNU General Public License as # published by the Free Software Foundation; # either version 2, or (at your option) any # later version. # # This program is distributed in the hope that # it will be useful, but WITHOUT ANY WARRANTY; # without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR # PURPOSE. See the GNU General Public License # for more details. # # You should have received a copy of the GNU # General Public License along with this # program; if not, see: # http://www.gnu.org/licenses/gpl.html */ /* Version: 1.0 */ /* Feb 2009 Many improvements thanks to Piotr T. Rozanski */ /* REFERENCES: Roukema & Rozanski 2009, Astronomy & Astrophysics, xxx, xxx http://arXiv.org/abs/0902.3402 */ /* avoid imaginary sqrt ambiguities */ radexpand: false; normalise_v4(x) := [ x[1]/sqrt(x.x), x[2]/sqrt(x.x), x[3]/sqrt(x.x), x[4]/sqrt(x.x) ]; /* normalised tangent at 4-vector a pointing towards 4-vec b */ tangent_S3_at_a(a, b) := normalise_v4(b - ((a . b) * a)); /* spherical-newtonian acceleration at a from image at b */ accel_at_a(a,b) := tangent_S3_at_a(a, b) / (1- (a . b)^2); /* zeroth image of massive particle is at [0,0,0,1] */ /* number related to the Golden Ratio */ phi_m14 : (sqrt(5)-1)/4 ; phi_p14 : (sqrt(5)+1)/4 ; /* = cos(%pi/5) */ /* 12 adjacent images of massive particle; modulus must be 1 since the physical space is S^3, not R^4 */ image1 : [ 0 , 1/2, -phi_m14, phi_p14 ]; image2 : [ 0 , -1/2, -phi_m14, phi_p14 ]; image3 : [ -1/2, phi_m14, 0, phi_p14 ]; image4 : [ 1/2, phi_m14, 0, phi_p14 ]; image5 : [ -phi_m14, 0, -1/2, phi_p14 ]; image6 : [ -phi_m14, 0, 1/2, phi_p14 ]; image7 : [ 0 , -1/2, phi_m14, phi_p14 ]; image8 : [ 0 , 1/2, phi_m14, phi_p14 ]; image9 : [ 1/2, -phi_m14, 0, phi_p14 ]; image10 : [ -1/2, -phi_m14, 0, phi_p14 ]; image11 : [ phi_m14, 0, 1/2, phi_p14 ]; image12 : [ phi_m14, 0, -1/2, phi_p14 ]; /* sum of accelerations */ accel_full(vec_m) := accel_at_a(vec_m,image1) + accel_at_a(vec_m,image2) + accel_at_a(vec_m,image3) + accel_at_a(vec_m,image4) + accel_at_a(vec_m,image5) + accel_at_a(vec_m,image6) + accel_at_a(vec_m,image7) + accel_at_a(vec_m,image8) + accel_at_a(vec_m,image9) + accel_at_a(vec_m,image10) + accel_at_a(vec_m,image11) + accel_at_a(vec_m,image12); /* test particle near to zeroth image of massive particle, in arbitrary direction (x,y,z); modulus must be 1 since the physical space is S^3, not R^4 */ vec_m : [ sin(rr)*x, sin(rr)*y , sin(rr)*z, cos(rr) ]; /* Taylor expand to fifth order */ tay : taylor(accel_full(vec_m),rr,0,5); /* simplify x, y components */ tay_x : gfactor(rectform(expand(tay[1]))), z : sqrt(1-x^2-y^2); tay_y : gfactor(rectform(expand(tay[2]))), x : sqrt(1-y^2-z^2); /* simplify z component; to retain symmetry, reverse alphabetical order of (x,z) by temporarily substituting (b,a); */ subst(a,z,subst(b,x,tay[3])); gfactor(rectform(expand(%))), y : sqrt(1-a^2-b^2); tay_z : subst(z,a,subst(x,b,%)); /* simplify w component */ tay_w : expand(tay[4]), x : sqrt(1-y^2-z^2); /* still fully algebraic */ residual_fully_alg : [tay_x, tay_y, tay_z, tay_w ]; /* 212-bit precision for big floats 1b0 etc. */ fpprec : 64; ratepsilon : 1b-64; /* numerical coefficients in front of x,y,z */ residual_xyz_alg : residual_fully_alg, bfloat; /* put some weakly pseudo-random numbers into x,y,z, while requiring the R^3 modulus to be 1 so that vec_m is on S^3; */ x: 2b0*random(1.0)-1b0; theta: 2b0*%pi *random(1.0); y: sin(acos(x)) * cos(theta), bfloat; z: sin(acos(x)) * sin(theta), bfloat; residual : residual_xyz_alg, bfloat; /* components */ radial_n: -tangent_S3_at_a(vec_m, [0, 0, 0, 1]), rr : 1b-3; residual_radial: expand(residual . radial_n); ortho: expand(residual - residual_radial * radial_n); residual_ortho: expand(sqrt(ortho.ortho)); /* result summary */ print("vector residual (coeff, x,y,z algeb.)"); residual_fully_alg ; print("vector residual (x,y,z algebraic)"); residual_xyz_alg; print("vector residual (x,y,z numeric)"); residual; print("radial residual (x,y,z numeric)"); residual_radial; print("orthogonal residual (x,y,z numeric)"); residual_ortho; /* ------------------------------------------- */