
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the gravitational potential      %
% relative to the rotating earth.            % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   Large scale view.
grav0 = 9.8066;
a=2*10^7/pi;
Omega=2*pi/(24*60*60);
x = linspace(0,8,101);
y = linspace(-4,4,100);
xx = a*x; yy = a*y;
[XX,YY] = meshdom(xx,yy);
r = sqrt(XX.^2+YY.^2);
Phi_e = -grav0*a^2./r;
Phi_c = -0.5*Omega^2*XX.^2;
Phi_00 = -grav0*a-0.5*Omega^2*a^2;
Phi = Phi_e+Phi_c;
Phi = max(Phi,Phi_00);
LogPhi = log(Phi);
XXP = XX/a; YYP = YY/a;
subplot(2,2,1)
contour(XXP,YYP,LogPhi,50)
axis([0,8,-4,4]);
axis('square')
title('GRAVITY POTENTIAL')
xlabel('Earth Radii')
%%%%%%%%%%%%%%%%%   Close-in view.
x = linspace(eps,1.1,101);
y = linspace(eps,1.1,101);
xx = a*x; yy = a*y;
[XX,YY] = meshdom(xx,yy);
r = sqrt(XX.^2+YY.^2);
Phi_e = -grav0*a^2./r;
Phi_c = -0.5*Omega^2*XX.^2;
Phi = Phi_e+Phi_c;
Phi = max(Phi,Phi_00);
Phi = Phi - min(min(Phi)) + 1; %%% to force contour at 0.
XXP = XX/a; YYP = YY/a;
subplot(2,2,2)
contour(XXP,YYP,Phi,20)
axis('square')
title('GRAVITY POTENTIAL')
xlabel('Earth Radii')
%%%%%%%%%%%%%%%%%   Horizontal variation.
 
z0 = 0; z1=10^4; Delz = z1-z0;
nump = 101;
phi=linspace(0,pi/2,nump);
phideg=(180/pi)*phi;
Delphi = (pi/2)/(nump-1);
phimid=linspace(Delphi/2,pi/2-Delphi/2,nump-1);
r1= a+z1;
RR = r1+0*phi;
xx = r1.*cos(phi);
Phi_e = -grav0*a^2./RR;
Phi_c = -0.5*Omega^2*xx.^2;
Phi = Phi_e+Phi_c;
Phi = max(Phi,Phi_00);
%%% Potential at surface.
r0= a+z0;
RR = r0+0*phi;
xx = r0.*cos(phi);
Phi_e = -grav0*a^2./RR;
Phi_c = -0.5*Omega^2*xx.^2;
Phi_surf = Phi_e+Phi_c;
Phi_surf = max(Phi_surf,Phi_00);
Phi_diff = Phi-Phi_surf;
subplot(2,2,3)
plot(phideg,Phi_diff)
axis([0,90,9.75*10^4, 9.80*10^4]);
title('Delta Phi(10km)')
xlabel('Degrees North')
%%%%%%%%%%%%%%%%%%%%  Horizontal gravity component.
phimid = phimid*(180/pi);
subplot(2,2,4)
  grav_H = -(1/a)*diff(Phi_diff)/Delphi;
  plot(phimid,grav_H)
end
axis([0,90,-6*10^(-5),0]);
title('Horizontal g (10km)')
xlabel('Degrees North')
print -deps 'gravxx.eps'
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%