This function predicts displacements from solid earth tides and is a wrapper function for Dennis Milbert's solid.exe program, which I believe is a wrapper for an adapted version of Dehant's DEHANTTIDEINEL.F Fortran code. For a more thoughtful explanation of the history of this program and of solid earth tides, visit Dennis Milbert's site.
[up,north,east] = earthtide(lat,lon,t) returns predicted solid earth tide displacements in the vertical, meridional, and zonal directions, in units of meters. Inputs lat and lon are decimal degrees and t is time (in normal people's civil time or UTC time, not the few-seconds-off GPS time which is required to run the solid.exe program directly). Time should can be in datenum format, or you can use a regular old string (e.g., 'March 2, 1990' ).
Calculate two weeks of vertical displacement in Austin, TX (30.25°N,97.75°W) at hourly resolution:
t = datenum('may 27, 1984'):1/24:datenum('june 10, 1984'); up = earthtide(30.25,-97.75,t);
Plot displacement in centimeters:
plot(t,100*up) box off datetick ylabel('vertical displacement (cm)') axis tight
Creating a map of displacements often requires more computation time than calculating displacements for an array of times at a single point. That's because this function runs the solid.exe program separately for each location and each unique day. In Example 1 above, 337 timesteps were calculated, but solid.exe was only called 15 times, because the time array spanned only 15 days. Here we calculate a global snapshot of solid earth tide displacements at an array of grid points, for a single snapshot in time. Although our grid is rather coarse, it contains 18x18=324 locations, so solid.exe need to run earthtide 324 times.
Note that although earthtide accepts multiple times as input in the form of an array each call of earthtide can only accept a single geographical point. Thus, a loop is required to calculate solid earth displacements for multiple points. This calculation takes about a minute and a half on my laptop computer.
% Create gridded arrays: [lon,lat] = meshgrid( -170:20:170,85:-10:-85); % Preallocate outputs for speed: up = NaN(size(lat)); north = NaN(size(lat)); east = NaN(size(lat)); for n = 1:numel(up) [up(n),north(n),east(n)] = earthtide(lat(n),lon(n),'may 16, 2007 9:23'); end
Plot the data:
figure pcolor(lon,lat,100*up) xlabel('longitude (deg)') ylabel('latitude (deg)') cb = colorbar('location','southoutside'); xlabel(cb,'vertical displacement (cm)') shading interp colormap(parula(256)) hold on quiver(lon,lat,east,north,'k');
This Matlab function was written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics (UTIG), but the real credit for this work should go to Veronique Dehant and Dennis Milbert. If any errors exist in this work (surely possible), those errors are most likely mine.