Quivers/streamlines not perpendicular to contours

I've created a mesh representing a superelevation transition along a highway - basically a rotating plane around X=0 as Y increases.

I'm trying to find the longest raindrop flow line, but the problem I'm having is that currently the contours and streamlines/quivers do not align. THat is, the streamlines are not perpendicular/normal to the contours. I'm not sure if the contours are incorrect or if the streamlines are incorrect, or if I'm interpreting the use of a function incorrectly. Some help would be much appreciated!

The code is below:

% Hardcoded Input Variables

left_width = 13.5; % Left width of the road in meters

right_width = 3; % Right width of the road in meters

grade = 0.02; % Longitudinal grade (e.g., 2%)

normal_crossfall = -0.03; % Normal crossfall (e.g., -3%)

superelevated_crossfall = 0.06; % Superelevated crossfall (e.g., 6%)

transition_length = 70; % Length of transition in meters

contour_interval = 0.20; % Contour interval (e.g., 200 mm)

% Mesh Parameters

num_segments = 50; % Number of segments along the transition length

num_width_points = 20; % Number of points across the road width (ensure this is even)

% Generate Mesh Points

x = linspace(-left_width, right_width, num_width_points); % Evenly spaced points across the total width

y = linspace(0, transition_length, num_segments + 1); % Only transition length

% Initialize Elevation Matrix

elevation = zeros(length(x), length(y)); % Adjust size based on unique x

% Compute Elevations

for i = 1:length(y)

segment_length = y(i);

for j = 1:length(x)

% Longitudinal grade component

longitudinal_elevation = grade * segment_length; % Extend grade throughout

% During transition

start_elevation = longitudinal_elevation + (normal_crossfall * (-x(j)));

end_elevation = longitudinal_elevation + (superelevated_crossfall * (-x(j)));

elevation(j, i) = interp1([0, transition_length], [start_elevation, end_elevation], segment_length);

end

end

% Create 2D Contour Plot

[X, Y] = meshgrid(x, y);

figure;

contourf(X, Y, elevation', 'LineColor', 'none'); % Filled contour plot

hold on;

% Add contour lines using the specified contour interval

contour_levels = min(elevation(:)):contour_interval:max(elevation(:));

contour(X, Y, elevation', contour_levels, 'LineColor', 'k'); % Contours in black

colorbar;

title('2D Contour Plot of Superelevation Transition');

xlabel('Width of Road (m)');

ylabel('Length of Transition (m)');

axis equal; % Set equal scaling for both axes

% Calculate gradients for flow direction

[grad_x, grad_y] = gradient(elevation'); % Calculate gradients of elevation

% Normalize gradients for better visualization

magnitude = sqrt(grad_x.^2 + grad_y.^2);

grad_x = -grad_x ./ magnitude; % Reverse direction for upstream to downstream

grad_y = -grad_y ./ magnitude; % Reverse direction for upstream to downstream

% Create Streamline starting from the highest values of Y along the Y-axis

lowest_x = min(x); % Find the lowest value of X

startY = transition_length:-5:0; % Start from downstream to upstream

streamlines = streamline(X, Y, grad_x, grad_y, lowest_x * ones(size(startY)), startY);

% Add quiver plot for direction visualization

quiver(X, Y, grad_x, grad_y, 'Color', 'r', 'AutoScale', 'off'); % Add quivers

hold off;

[Trav]

This is a snapshot of the figure I get - the thick blue lines are drawn on by me manually and are what I would expect the streamlines and quivers to follow, according to the contours (in black).

Answer 1

[Bruno Luong]

Add the command

axis equal

Different aspect ratio in x and y can skew the perpendicularity visually.

[Trav]

axis equal is already being used. I think I figured it out - had some incorrect calcs which were skewing the results.

Thanks