Negative Gravity Gradient

[CalLightmanCode]

Dear Jason M. Pearl:
I respect this project you've developed. And I have learned about your paper "A fast quadrature-based gravity model for the homogeneous polyhedron". It is a wonderful work.
However, an unexpected Negative Gravity Gradient may appear in the function "gravityGradient" in the catalog src/GravityModels. To avoid confusion, I mean negative gravitational potential, for example, U = -Mu / R, where U is the gravity potential, Mu is the gravity parameter, and R = [X, Y, Z].
I suspect it may originate from an unusual definition of r = x - p which is defined in the aforementioned paper.
I offer my suggestions for you to look over and refer to. And I hope I don't bother you.
Best Wishes.

[jmpearl]

I may be misunderstanding the issue (please let me know if this is the case) but the gravitational potential is negative by common convention since it takes work to remove a mass from a gravitation field. It is certainly possible to define r as p-x, I went with r = x-p because that is what Werner and Scheeres did in there 1996 work "Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia".

By gravitational gradient did you mean there's an issue with the gravitational gradient tensor?

[CalLightmanCode]

Yes, I mean that there's an issue with the gravitational gradient tensor.
There is an opposite sign in the gravitational gradient tensor.
This is a validation code for your reference::

clear all
clc
close all
addpath(genpath("GaMA-main"))
% This is a validation code for the gravitational gradient tensor.
%% A benchmarking model: mass point
syms mu % the gravity parameter of the asteroid
syms x y z % the components of the position vector
r = sqrt(x * x + y * y + z * z);
V = -mu / r; % the negative gravitational potential of a mass point model
% In this part, I have derived the partial derivative of the potential.
q = 2; % differential order
pV = sym(zeros(q + 1, q + 1, q + 1));
for i = 0 : q
for j = 0 : (q - i)
for k = 0 : (q - i - j)
% the 3-dimension array contains the partial derivative of the potential.
% for example, pV(2, 2, 1) refers to $\frac{\partial^{2}V}{\partial x\partial y}$
pV(i + 1, j + 1, k + 1) = diff(diff(diff(V, x, i), y, j), z, k);
end
end
end
func_GravityTensor = matlabFunction(pV);
%% The gravity models in GaMA-main
G = 6.6710^-11; % gravitational constant (N m2/kg2)
load('Eros.mat'); % load stored Eros properties
Mu = bodyProperties.massG; % set stand grav parameter
meshfile = 'Eros_7624.obj'; % mesh file to load
mesh = SurfaceMesh(meshfile); % create the surface mesh object
mesh.coarsen(3000); % coarsen it to 3k faces
numMascons = 1000; % number of mascons

quadratureModel = ApproximatePolyhedralModel(mesh, Mu, 'B2');
polyhedralModel = AnalyticPolyhedralModel(mesh, Mu);
masconModel = MasconModel(mesh, Mu, numMascons);
%% the mass point model
R = [1 2 -3] * 1e6; % if the field point is far enough, the above gravity models should approximate to a mass point model
GravityTensor_MassPoint = func_GravityTensor(Mu, R(1), R(2), R(3));

disp('================================================================')
disp('Gravity potential')
potQ = quadratureModel.potential(R);
potP = polyhedralModel.potential(R);
potM = masconModel.potential(R);
potMP = GravityTensor_MassPoint(1, 1, 1);

disp(['quadratureModel::Potential = ', num2str(potQ,'%.8e')])
disp(['polyhedralModel::Potential = ', num2str(potP,'%.8e')])
disp(['masconModel::Potential = ', num2str(potM,'%.8e')])
disp(['MassPoint::Potential = ', num2str(potMP,'%.8e')])

disp('================================================================')
disp('Gravity Acceleration')
accQ = quadratureModel.acceleration(R);
accP = polyhedralModel.acceleration(R);
accM = masconModel.acceleration(R);
accMP = -[GravityTensor_MassPoint(2, 1, 1) GravityTensor_MassPoint(1, 2, 1) GravityTensor_MassPoint(1, 1, 2)];
% It should be noted that the acceleration equals to the nagative potential
% gradient. And the ``minus" sign is present.
disp(['quadratureModel::Acceleration = ',num2str(accQ,'%.8e,')])
disp(['polyhedralModel::Acceleration = ',num2str(accP,'%.8e,')])
disp(['masconModel::Acceleration = ',num2str(accM,'%.8e,')])
disp(['MassPoint::Potential = ',num2str(accMP,'%.8e,')])

disp('================================================================')
disp('GravGradient')
gravGradQ = quadratureModel.gravityGradient(R);
gravGradP = polyhedralModel.gravityGradient(R);
gravGradM = masconModel.gravityGradient(R);

gravGradMP = [GravityTensor_MassPoint(3, 1, 1) GravityTensor_MassPoint(2, 2, 1) GravityTensor_MassPoint(2, 1, 2);...
GravityTensor_MassPoint(2, 2, 1) GravityTensor_MassPoint(1, 3, 1) GravityTensor_MassPoint(1, 2, 2);...
GravityTensor_MassPoint(2, 1, 2) GravityTensor_MassPoint(1, 2, 2) GravityTensor_MassPoint(1, 1, 3)];

disp("quadratureModel::")
disp(['GravGradient = |',num2str(gravGradQ(1,[1,2,3]),'%.8e, '),'|'])
disp([' |',num2str(gravGradQ(1,[2,4,5]),'%.8e, '),'|'])
disp([' |',num2str(gravGradQ(1,[3,5,6]),'%.8e, '),'|'])

disp("polyhedralModel::")
disp(['GravGradient = |',num2str(gravGradP(1,[1,2,3]),'%.8e, '),'|'])
disp([' |',num2str(gravGradP(1,[2,4,5]),'%.8e, '),'|'])
disp([' |',num2str(gravGradP(1,[3,5,6]),'%.8e, '),'|'])

disp("masconModel::")
disp(['GravGradient = |',num2str(gravGradM(1,[1,2,3]),'%.8e, '),'|'])
disp([' |',num2str(gravGradM(1,[2,4,5]),'%.8e, '),'|'])
disp([' |',num2str(gravGradM(1,[3,5,6]),'%.8e, '),'|'])

disp("MassPoint::")
disp(['GravGradient = |',num2str(gravGradMP([1,2,3]),'%.8e, '),'|'])
disp([' |',num2str(gravGradMP([4,5,6]),'%.8e, '),'|'])
disp([' |',num2str(gravGradMP([7,8,9]),'%.8e, '),'|'])

disp("Conclusion: There is an opposite sign in GravGradient between the mass point model and the GaMA models.")
disp("I provide this code for your refence.")
disp("I really respect this project you've developed. It is an excellent work!")
disp("Best wishes.")