%  Matlab m-file to take the derivative of a geoid height 
%
%  The first column of the file is time (not used here).
%  The second column of the file is latitude in degrees.
%  The third column of the file is longitude in degrees.
%  The fourth column is geoid height in meters.
%  The fifth column is gravity anomaly in milligals.
%  The sixth column is uncertainty in gravity anomaly.
%  (note this is a descending track so latitude decreases
%  with increasing time.)
%
%  There are 512 evenly spaced points.  The total length of the
%  profile is 1705000 m.
%
%  Get the hw1part.m file and the data file using fetch or ftp
%  topex.ucsd.edu
%  cd pub/class/hw1
%  mget *
%
   load geosatd.dat
   lat=geosatd(:,2);
   lon=geosatd(:,3);
   geoid=geosatd(:,4);
%
%  remove the mean and window the profile to minimize edge effects
%
   nx=length(lat);
   window=hanning(nx);
   geoid=(geoid-mean(geoid)).*window;
%
%  plot geoid height
%
   subplot(3,1,1), plot(lat,geoid);
   ylabel('geoid height (m)')
%
%  get the profile length and data spacing
%
   L=1705000.;
   dx=L/nx;
%  
%  compute the derivative using the first difference formula.
%  duplicate the last point to extend the length from 511 to 512.
%  note there is a shift
%
   vd = diff(geoid)/dx;
   nd = length(vd);
   vd = [vd;vd(nd:nd)];
%
%  now compute the derivative using the derivative property.
%  first compute FFT and do an fftshift to make it easy to gerenate
%  wavenumbers.
%
   cg = fftshift(fft(geoid));
%
%  generate the wavenumber.  if these wavenumbers are wrong the
%  answer will be wrong and the imaginary part of the result will
%  be large.
%
   k = -nx/2:(nx/2-1);
   k = k'./L;
%  
% compute derivative using fft.  don't forget to shift the
% memory locations back after the multiply
%
   cdg = i*2*pi*k.*cg;
   rdg = real(ifft(fftshift(cdg)));
%
%  plot the two derivatives as well as their difference
%
   subplot(3,1,2), plot(lat,vd,lat,rdg)
   ylabel(' geoid slope ')
   subplot(3,1,3), plot(lat,vd-rdg)
   ylabel(' slope difference ')
   xlabel('latitude (deg)')
   pause
%
% now see if we can do a better job by shifting 1/2 data point
%
   cdg = i*2*pi*k.*cg.*exp(i*2*pi*k*dx/2);
   rdg = real(ifft(fftshift(cdg)));
%
%  plot the curves again
%
   subplot(3,1,1), plot(lat,geoid);
   ylabel('geoid height (m)')
   subplot(3,1,2), plot(lat,vd,lat,rdg)
   ylabel(' geoid slope ')
   subplot(3,1,3), plot(lat,vd-rdg)
   ylabel(' slope difference ')
   xlabel('latitude (deg)')
