estimate_ts.m 2.79 KB
Newer Older
1
function [Mobj] = estimate_ts(Mobj,u,zeta)
Pierre Cazenave's avatar
Pierre Cazenave committed
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

% Estimate time step at each node
%
% [Mobj] = function estimate_ts(Mobj)
%
% DESCRIPTION:
%    Calculate barotropic time step
%
% INPUT
%    Mobj = matlab mesh object
%
% OUTPUT:
%    Mobj = matlab structure containing mesh time step
%
% EXAMPLE USAGE
%    Mobj = estimate_ts(Mobj)
%
% Author(s):
%    Geoff Cowles (University of Massachusetts Dartmouth)
%
% Revision history
23 24 25
%    2012-07-14 Add great circle approximation if only provided with
%    latitude and longitudes. Also add arguments to the function to define
%    current velocity and tidal amplitudes.
Pierre Cazenave's avatar
Pierre Cazenave committed
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
%
%==============================================================================

subname = 'estimate_ts';
global ftbverbose
if(ftbverbose)
    fprintf('\n')
    fprintf(['begin : ' subname '\n'])
end;

%------------------------------------------------------------------------------
% Set constants
%------------------------------------------------------------------------------

g    = 9.81; %gravitational acceleration
41 42
%u    = 3.0;  %u-velocity
%zeta = 11.0;  %tide amp
Pierre Cazenave's avatar
Pierre Cazenave committed
43 44 45 46 47 48 49 50

if(~Mobj.have_bath)
    error('can''t estimate the time step without bathymetry')
end;

%------------------------------------------------------------------------------
% Compute the time step estimate
%------------------------------------------------------------------------------
51 52 53
if Mobj.have_xy
    x = Mobj.x;
    y = Mobj.y;
Pierre Cazenave's avatar
Pierre Cazenave committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
else
    % Will convert to metres when calculating element edge length
    x = Mobj.lon;
    y = Mobj.lat;
end
h = max(Mobj.h,0);
tri = Mobj.tri;
nVerts = Mobj.nVerts;
nElems = Mobj.nElems;

ts = ones(nVerts,1)*1e9;
lside = zeros(nVerts,1);
for i=1:nElems
    n1 = tri(i,1);
    n2 = tri(i,2);
    n3 = tri(i,3);
    nds = [n1 n2 n3];
    % Check whether we have x and y values and use great circle
    % approximations if we don't.
    if Mobj.have_xy
        lside(i) = sqrt( (x(n1)-x(n2))^2 + (y(n1)-y(n2))^2);
    else
        lside(i) = haversine(x(n1),y(n1),x(n2),y(n2));
    end
    dpth  = max(h(nds))+zeta;
    dpth  = max(dpth,1);
    ts(nds) = min(ts(nds),lside(i)/(sqrt(g*dpth) + u));
end;
if(ftbverbose); fprintf('minimum time step: %f seconds\n',min(ts)); end;
Mobj.ts = ts;
Mobj.have_ts = true;

if(ftbverbose)
    fprintf(['end   : ' subname '\n'])
end

function [km]=haversine(lat1,lon1,lat2,lon2)
% Haversine function to calculate first order distance measurement. Assumes
% spherical Earth surface. Lifted from:
93
%
Pierre Cazenave's avatar
Pierre Cazenave committed
94
% http://www.mathworks.com/matlabcentral/fileexchange/27785
95
%
Pierre Cazenave's avatar
Pierre Cazenave committed
96 97 98 99 100 101 102 103
R = 6371000;                    % Earth's mean radius in metres
delta_lat = lat2 - lat1;        % difference in latitude
delta_lon = lon2 - lon1;        % difference in longitude
a = sin(delta_lat/2)^2 + cos(lat1) * cos(lat2) * ...
    sin(delta_lon/2)^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
km = R * c;                     % distance in metres