%
%  Matlab m-file to load and plot a topography and gravity profile.
%
   load NimmoStack.dat
   xN=NimmoStack(:,1);
   TN=NimmoStack(:,2);
%
%  now generate the topographic load
%
   xl=[-512,-20,-16,12,25,511];
   yl=[0., 0., 350., 400., 0., 0.];
%
%  plot the topography profile 
%
   plot(xN-40.,TN,'.',xl,yl);
   xlabel('distance (km)')
   ylabel('topography (m)')
%
%  now interpolate the load to a 1 km spacing
%
   x=-512:1:511;
   ax=interp1(xl,yl,x);
%
%  plot the topography profile 
%
   plot(xN-40.,TN,'.',x,ax);
   xlabel('distance (km)')
   ylabel('topography (m)')
%
%  make the flexural rigidity
%
g=1.3;
rho=900;
E=9.e9;
nu=0.3;
Te=2900;
D=E*Te.^3/(12*(1-nu*nu));
%
%  now fourier transform the load
%
   Ak=fftshift(fft(ax));
   k=-512:1:511;
   L=1000*length(x);
   k=(k./L)';
%
%  now generate the transfer function
%
   arg=(D*(2.*pi*k).^4/(rho*g) + 1.);
   phi=-1./arg;
%
%  compute topography
%
   Wk=phi'.*Ak;
   wx=real(ifft(fftshift(Wk)));
   tx=ax+wx;
%
%  plot these things
%
   plot(x,ax,x,wx,x,tx,xN-40,TN,'.')
   xlim([-100,100])
   xlabel('distance (km)')
   ylabel('topography (m)')
