function [C,cumlen,dCdS,dS] = curvars(Xs, Ys, bc1, bcset, smoothnum) % keyboard % Inputs: % X and Y coordintes of stream centerline (or bankline) % Outputs: % C - curvature calculated via formulas of Johanesson and Parker 1985 % followed by Motta 2012 % cumlen - cumulative stream length in the streamwise coordinate direction % dCdS - dC/ds required by hydrodynamic model in calculation of a_p_5 % del_s - stream length ds between nodes along stream centerline % Xs=X;Ys=Y; % For testing only % Reshape inputs Xs = Xs(:); Ys = Ys(:); n_in = nargin; if n_in == 2 % if only X,Y coordinates are given, set upstream C=0 bcset = 2; smoothnum = 0; end %% Compute radius of circle fit through three consecutive points (other methods coded at end of function) x1 = Xs(1:end-2); y1 = Ys(1:end-2); x2 = Xs(2:end-1); y2 = Ys(2:end-1); x3 = Xs(3:end); y3 = Ys(3:end); At = 0.5 * ((x2-x1).*(y3-y1) - (y2-y1).*(x3-x1)); % signed areas of the triangle that has three consecutive points as vertices sides_prod = (((x2-x1).^2+(y2-y1).^2).*((x3-x1).^2+(y3-y1).^2).*((x3-x2).^2+(y3-y2).^2)).^.5; C = -4*At./sides_prod; % Upstream boundary condition for curvature can be set different ways if bcset == 1 CS = bc1; % 'periodic' boundary condition which imposes the curvature at the end of the river on the previous time step to the beginning of the river at the next elseif bcset == 2 CS = 0; % C(1)=0. This causes the most-upstream portion of the river to straighten and elongate over time elseif bcset == 3 CS = mean(abs(C)) + std(abs(C))*randn; % random upstream curvature (distribution has same mean and std as C) else CS = 2*C(1)-C(2); % Do a linear interpolation. This seems to have the same effect as as setting=0, but haven't tested it over large times extensively. end % Downstream boundary condition for curvature is simply set to zero: CE = 0; C = vertcat(CS, C, CE); %% Calculate stream length variables dS and cumlen dS_btwn_nodes = sqrt(diff(Xs).^2+diff(Ys).^2); dS_at_nodes = .5*(dS_btwn_nodes(1:end-1)+dS_btwn_nodes(2:end));% want dS at a node. this calculates it as the distance between a point and the next point. therefore, average either side to get dS at node. dS = [dS_at_nodes(1); dS_at_nodes; dS_at_nodes(end)]; % set dS boundaries equal to their nearest neighbors cumlen = (vertcat(0,cumsum(dS_btwn_nodes))); % this ensures that s=0 at the beginning of the curve %% dcds by central differencing dC = diff(C); dCdS = dC./dS_btwn_nodes; dCdS = [dCdS; dCdS(end)]; %% Smooth Curvature % Curvature smoothing as in Motta, eq'n 16, following Crosato 1990. for u = 1:smoothnum % for u = 1:5 Cs = C(2:end-1); % don't want boundary conditions to propogate C1 = Cs(1:end-2); C2 = Cs(2:end-1); C3 = Cs(3:end); Cs = vertcat(C(1), (2*Cs(1)+Cs(2))/3, (C1+2*C2+C3)/4, (2*Cs(end)+Cs(end-1))/3, C(end)); % reinsert boundary conditions C = Cs; end end %%----EXTRA CODE -- OTHER METHODS OF CURVATURE CALCULATION---%% %% Calculate Curvature: Central differencing method (Johanneson and Parker 1989) % Xi = Xs(2:end-1); Yi = Ys(2:end-1); % Ximinus = Xs(1:end-2); Yiminus = Ys(1:end-2); % Xiplus = Xs(3:end); Yiplus = Ys(3:end); % % del_s = sqrt((Xi-Ximinus).^2+(Yi-Yiminus).^2); % del_s1 = sqrt((Xiplus-Xi).^2+(Yiplus-Yi).^2); % dxds = (Xiplus - Ximinus) ./ (del_s1+del_s); % dyds = (Yiplus - Yiminus) ./ (del_s1+del_s); % d2xds2 = (Xiplus+Ximinus-2.*Xi) ./ (0.5*(del_s1+del_s)).^2; % d2yds2 = (Yiplus+Yiminus-2.*Yi) ./ (0.5*(del_s1+del_s)).^2; % C = -(dxds.*d2yds2 - d2xds2.*dyds); % % CS = interp1(1:2,C(1:2),0,'linear','extrap'); % 'spline' or 'linear' % CE = interp1((length(C)-1):length(C),[C(end-1),C(end)],length(C)+1,'linear','extrap'); % C = [0 C CE]; % % del_s_S = interp1(1:2,del_s(1:2),0,'linear','extrap'); % del_s_E = interp1((length(del_s)-1):length(del_s),[del_s(end-1),del_s(end)],length(del_s)+1,'linear','extrap'); % del_s = [del_s_S del_s del_s_E]; %% Calculate Curvature: least squares regression to circle fitting % num_pts_fit = 5; % must be odd! % N = floor(num_pts_fit/2); % for o = 1:length(Xs)-num_pts_fit+1 % Xls = Xs(o:o+num_pts_fit-1); % Yls = Ys(o:o+num_pts_fit-1); % [z(:,o), r(o),~] = fitcircle([Xls;Yls]); % C_fit(N+o) = 1/r(o); % end % % http://stackoverflow.com/questions/3461453/determine-which-side-of-a-line-a-point-lies % Xs_fit = Xs(N:end-N); % Ys_fit = Ys(N:end-N); % X1 = Xs_fit(1:end-1); % X2 = Xs_fit(2:end); % Y1 = Ys_fit(1:end-1); % Y2 = Ys_fit(2:end); % X3 = z(1,:); % Y3 = z(2,:); % signs = sign((X2-X1).*(Y3-Y1)-(Y2-Y1).*(X3-X1)); % % C_fit(1:floor(num_pts_fit/2))=[]; % % C_fit(end+floor(num_pts_fit/2)-2) = 0; % C_fit = C_fit.*signs; % % % % Since curvature cannot be calculated at the endpoints, fitting must be done % % Linear interpolation % % CS = interp1(1:2,C(1:2),0,'linear','extrap'); % 'spline' or 'linear' % % CE = interp1((length(C)-1):length(C),[C(end-1),C(end)],length(C)+1,'linear','extrap'); % C = -[0 0 C_fit 0 0]; % % Xi = Xs(2:end-1); Yi = Ys(2:end-1); % Ximinus = Xs(1:end-2); Yiminus = Ys(1:end-2); % Xiplus = Xs(3:end); Yiplus = Ys(3:end); % del_s = sqrt((Xi-Ximinus).^2+(Yi-Yiminus).^2); % del_s_S = interp1(1:2,del_s(1:2),0,'linear','extrap'); % del_s_E = interp1((length(del_s)-1):length(del_s),[del_s(end-1),del_s(end)],length(del_s)+1,'linear','extrap'); % del_s = [del_s_S del_s del_s_E]; %% Loop for computing area and product method % C = []; % for g = 2:numel(Xs)-1 % p1 = [Xs(g-1) Ys(g-1)]; p2 = [Xs(g) Ys(g)]; p3 = [Xs(g+1) Ys(g+1)]; % A(g) = 0.5 * ( (p2(1)-p1(1))*(p3(2)-p1(2)) - (p2(2)-p1(2))*(p3(1)-p1(1)) ); % signed area % C(g) = -2*A(g) ./ sqrt(((p2(1)-p1(1)).^2+(p2(2)-p1(2)).^2)*((p3(1)-p1(1)).^2+(p3(2)-p1(2)).^2)*((p3(1)-p2(1)).^2+(p3(2)-p2(2)).^2)); % end %% Calculate Curvature: 3-pt circle fitting method - this method is unstable and produces very jagged curvatures % P1 = [Xs(3:end) Ys(3:end)]; % P2 = [Xs(2:end-1) Ys(2:end-1)]; % P3 = [Xs(1:end-2) Ys(1:end-2)]; % % x1y1 = P2-P1; % x2y2 = P3-P2; % x3y3 = P3-P1; % % norm1 = sqrt(sum(x1y1.^2,2)); % norm3 = sqrt(sum(x3y3.^2,2)); % % determinant = x1y1(:,1).*x2y2(:,2)-x1y1(:,2).*x2y2(:,1); % term3 = sqrt(sum(norm3.^2,2)); % term2 = [x2y2(:,1).*term3 x2y2(:,2).*term3]; % normterms = sqrt(sum(term2.^2,2)); % den = norm1.*normterms; % % C = 2*determinant./den;