/*
# 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;

/* ------------------------------------------- */