% tidalaccel.m
function tidalaccel()  % (a function merely so I can use a subfunction inside)

% Simple comparison of tidal accelerations of the Earth due to various solar
% system objects and then of Europa due to Jupiter.
% Andy Ganse, APL-UW, 2005
% aganse@apl.washington.edu, http://staff.washington.edu/aganse
% Reference: 
% Section 2.3.3 of William Lowrie's Fundamentals of Geophysics, Cambridge
% University Press, 1997.
% (highly recommended undergrad or intro graduate geophysics textbook)


% First define solar system constants like radii and masses (all in MKS units).
% source: planetary "quickfacts" on NASA/JPL website

% Solar system masses in kg
Msun=1.9891e30;
Mmercury=3.3034e23;
Mvenus=4.8685e24;
Mearth=5.9736e24;
Mmoon=7.3483e22;
Mmars=1.1171e24;
Mjupiter=1.8986e27;
Msaturn=5.6846e26;
Muranus=8.6832e25;
Mneptune=1.0243e26;
Mpluto=1.2545e22;

% Mean solar system body radii in m
Rearth=6371e3;
Reuropa=1569e3;

% Mean orbital radii in m
Rorb_mercury=57.95e9;
Rorb_venus=108.11e9;
Rorb_earth=149.57e9;
Rorb_moon=382e6;
Rorb_mars=227.84e9;
Rorb_jupiter=778.14e9;
Rorb_europa=670.9e6;
Rorb_saturn=1.4270e12;
Rorb_uranus=2.8783e12;
Rorb_neptune=4.4999e12;
Rorb_pluto=5.9130e12;

name=['Sun    ';'Moon   '; 'Mercury'; 'Venus  '; 'Mars   '; 'Jupiter'; ...
      'Saturn '; 'Uranus '; 'Neptune'; 'Pluto  '];

dist=[Rorb_earth, Rorb_moon, abs(Rorb_mercury-Rorb_earth), ...
        abs(Rorb_venus-Rorb_earth), ...
        abs(Rorb_mars-Rorb_earth), abs(Rorb_jupiter-Rorb_earth), ...
        abs(Rorb_saturn-Rorb_earth), abs(Rorb_uranus-Rorb_earth), ...
        abs(Rorb_neptune-Rorb_earth), abs(Rorb_pluto-Rorb_earth)]';

mass=...
[Msun,Mmoon,Mmercury,Mvenus,Mmars,Mjupiter,Msaturn,Muranus,Mneptune,Mpluto]';

disp(' ')
disp('Tidal accels (m/s^2) at pts A & B on Earth due to solar system bodies.')
disp('----------------------------------------------------------------------')
disp('due to      a_T at A      a_T at B ')
disp('--------   -----------   -----------')
%disp('due to       a_T @ A       a_T @ B       approx @ A      approx @ B ')
%disp('--------   -----------   -----------     -----------     -----------')
[aA,aB]=accels(mass,dist,Rearth);
%[aAo,aBo]=approx(mass,dist,Rearth);
for i=1:10
  fprintf(1,'%s    %g \t %g\n',name(i,:),aA(i),aB(i));
% fprintf(1,'%s    %g \t %g \t %g \t %g\n',name(i,:),aA(i),aB(i),aAo(i),aBo(i));
end

disp(' ')
disp('Tidal accel (m/s^2) at pts A & B on Europa due to Jupiter.')
disp('----------------------------------------------------------')
disp('due to       a_T @ A       a_T @ B')
disp('--------   -----------   -----------')
[aA,aB]=accels(Mjupiter,Rorb_europa,Reuropa);
fprintf(1,'%s    %g \t %g\n','Jupiter',aA,aB);
disp(' ')




function [aA,aB]=accels(mass,separation,radius)
G=6.6726e-11;  % grav const in m^3/(s^2 kg)
% eqn2.30:
aA=G*mass.*(1./separation.^2-1./(separation+radius).^2);  % a_T at point A
% eqn2.28:
aB=G*mass.*(1./(separation-radius).^2-1./separation.^2);  % a_T at point B


function [aA,aB]=approx(mass,separation,radius)
G=6.6726e-11;  % grav const in m^3/(s^2 kg)
% eqn2.31:
aA=G*mass./separation.^2.*(2*radius./separation-3*(radius./separation).^2);
% eqn2.29:
aB=G*mass./separation.^2.*(2*radius./separation+3*(radius./separation).^2);