function [Mobj] = estimate_ts(Mobj,u,zeta) % 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 % 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. % %============================================================================== subname = 'estimate_ts'; global ftbverbose if(ftbverbose) fprintf('\n') fprintf(['begin : ' subname '\n']) end; %------------------------------------------------------------------------------ % Set constants %------------------------------------------------------------------------------ g = 9.81; %gravitational acceleration %u = 3.0; %u-velocity %zeta = 11.0; %tide amp if(~Mobj.have_bath) error('can''t estimate the time step without bathymetry') end; %------------------------------------------------------------------------------ % Compute the time step estimate %------------------------------------------------------------------------------ if Mobj.have_xy x = Mobj.x; y = Mobj.y; 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: % % http://www.mathworks.com/matlabcentral/fileexchange/27785 % 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