function [loga, ax, ay, a_theta_x, a_theta_y, a_radial_x, a_radial_y] = lagrange_physics(Gy, Gx, tolerance)
% Do Newtonian gravity calculations.
% tolerance = log10( acceleration in Mm/s^3)
%
% data from p. A5 of
% _Fundamentals of Physics_ Halliday and Resnick 1988.
% 
% mean distance from Earth to Moon: 3.82x10^8 m = 382 Mm
% mean distance from Earth to the sun: 1.50x10^11 m
% Mass:
% Sun = 1.99x10^30 Kg
% Earth = 5.98x10^24 Kg
% Luna = 7.36x10^22 Kg
% Mean radius
% Sun = 6.96x10^8 m
% Earth = 6.37x10^6 m = 40 Mm/(2*pi) (by definition)
% Luna = 1.74x10^6 m = 1.74 Mm
% Gravitational Constant = G = 6.67x10^-11 m^3 / s^2 Kg
% = 6.67e-11 m^3 / s^2 Kg = 66.7e-12 m^3 / s^2 Kg = 66.7e-6 m^3 / s^2 Gg.
G = 66.7e-24; % G = 66.7e-24 Mm^3 /( s^2 Gg ).

% Make sure that these numbers are not exact multiples of eps
% to avoid dividing by zero.
% Perhaps also include the radius of Earth and moon.
'setting up coordinates'

Earth_mass = 5.98e18; % Gg
Luna_mass = 73.6e15 ; % Gg

Luna_to_Earth = 382e0; % Mm
Lunax = Luna_to_Earth;
Lunay = 0;

day =  24 * 60 * 60;  % s = 1 day * 24 hours/day * 60 minutes/hour * 60 seconds/minute
month = 27.3 * day;  % s = 1 month * 27.3 days/month
omega = 2*pi/month; % rad/sec

% Tweak Earth-Moon coordinates such that
% the center of mass of the entire system (what they both orbit)
% is at 0,0 of our rotating coordinate frame:
% Acceleration of earth = centroid_to_Earth_vector * omega^2 =
% = Earth_to_Luna_vector * G * Luna_mass / (Earth_to_Luna^3).

temp = G * Luna_mass /( (Luna_to_Earth)^3 * omega^2 );

Earthx = -temp * Lunax;
Earthy = -temp * Lunay;
% since we moved Earth away from Luna, move Luna the same amount back towards Earth.
Lunax = Lunax - Earthx;
Lunay = Lunay - Earthy;

% Set up x,y, and r2 coordinates relative to Earth.
[Ey, Ex] =  ndgrid( Gy - Earthy, Gx - Earthx );
Er2 = (Ex.^2) + (Ey.^2);
% Er2 = the square of the distance from Earth.

% set up x,y coordinates relative to Luna.
[Ly,Lx] =  ndgrid( Gy - Lunay, Gx - Lunax);
Lr2 = (Lx.^2) + (Ly.^2);
% Lr2 = the square of the distance from Luna
'calculating accelerations'
% acceleration due to gravity of Earth
% a = G * m / r^2 = G * m * r * (r^2)^(-3/2)
Aex = -(Ex .* (Er2 .^(-3/2) )) .* G .* Earth_mass; % Mm/s^2
Aey = -(Ey .* (Er2 .^(-3/2) )) .* G .* Earth_mass; % Mm/s^2

% acceleration due to gravity of moon
Alx = -Lx .* (Lr2 .^(-3/2)) .* G .* Luna_mass; % Mm/s^2
Aly = -Ly .* (Lr2 .^(-3/2)) .* G .* Luna_mass; % Mm/s^2

% The acceleration of objects with zero velocity in this rotating coordinate frame.
% the so-called "centrifugal force" -- only correct for
% objects with zero velocity in this rotating coordinate frame.
% FIXME: use the full dynamics formula (coriolis force and all)
% so we can investigate elliptical orbits, comet paths.
% In the rotating reference frame attached to the moon,
% a = r * omega^2.
[Cy, Cx] = ndgrid( Gy , Gx );
Acx = Cx .* omega^2; % Mm/s^2
Acy = Cy .* omega^2; % Mm/s^2


'getting total acceleration'
% total up all the accelerations on dust particles traveling at
% the specified velocity in that region
ay = (Acy + Aey + Aly);
ax = (Acx + Aex + Alx);
% Using log base 10.
loga = log( sqrt(ay.^2 + ax.^2) )./(log(10));




% Split each vector up into radial and angular components
% get unit radial vector
% (Assume that grid is *not* lined up to touch 0,0).
Uy = Cy ./( sqrt(Cx.^2 + Cy.^2) );
Ux = Cx ./( sqrt(Cx.^2 + Cy.^2) );
% dot product with unit radial vector
a_radial = (ax.*Ux + ay.*Uy);
a_radial_x = a_radial.*Ux;
a_radial_y = a_radial.*Uy;

% subtract off radial vector to leave only tangent vector
a_theta_x = ax - a_radial_x;
a_theta_y = ay - a_radial_y;

% cross product
% a_theta_x = a_x * 

if(0),
% To save RAM, throw away "uninteresting" data values
% larger than tolerance
I = find( tolerance < loga );
ay(I) = 0;
ay = sparse(ay);
ax(I) = 0;
ax = sparse(ax);
end;

% end lagrange_physics.m