How to find intersections in contour plot?
17 次查看(过去 30 天)
显示 更早的评论
I've been trying to find the points where this contour intersects the line y = -3, but it's not working. I tried just counting the intersection points and got 0, when it should clearly be 2. How do I fix this?
This is my entire code btw, in case you want to try running it. The part where I'm trying to find the intersections starts at line 51, where I've divided the code sections.
clear all
close all
clc
%Constant
rho = 4420; %kg/m^3
Cp = 550; %J/kg?K
T0 = 303.15; %K
A = 0.5; %[Absorbtivity]
k = 7.2; %W/m/K
alpha = 2.96*10^-6; %m^2/s
D = alpha;
P = 500; %W
u = 0.5; %m/s
Tm = 1933; %K
d_laser = 0.0001; %m
r_laser = d_laser/2; %m
a = r_laser;
p = D/(u*a);
%Define
n = 100;
% x = linspace(-0.00015,0.00025,n);
% z = linspace(-0.000005,0,n);
x = linspace(-0.00025,0.0012,n);
z = linspace(-0.0001,0,n);
%Normalized
x_nor = x/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh,z_mesh] = meshgrid(x_nor,z_nor);
y_mesh = 0;
fun = @(t) exp((-z_mesh.^2./(4*t))-((y_mesh.^2+(x_mesh-t).^2)./(4*p.*t+1)))./((4.*p.*t+1).*sqrt(t));
g = integral(fun,0,Inf,'ArrayValued',true);
% Ts = (A*P)/(pi*rho*Cp*sqrt(D*u*(a^3)))
Ts = Tm/0.16713
g2 = Tm/Ts
figure
hold on
[C,h] = contourf(x_mesh,z_mesh,g, 'ShowText', 'On')
[C1,h1] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList=0.16713,Color='red')
[C2,h2] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList = g2,Color='white')
hold off
clabel(C,h,'Color','white')
clabel(C1,h1,'Color','white')
clabel(C2,h2,'Color','white')
% xticks(-5:5:25)
% yticks(-10:2:0)
% axis([-5 25 -10 0])
title('x\primez\prime plane')
xlabel('x\prime')
ylabel('z\prime')
colorbar
%Extract data from contour plot
contourdata = contourc(x_nor,z_nor,g,[g2 g2])
% Initialize arrays to store the boundary coordinates
x_prime = [contourdata(1,2:size(contourdata,2))]
z_prime = [contourdata(2,2:size(contourdata,2))]
outlineplot = plot(x_prime,z_prime)
x_outline = outlineplot.XData
z_outline = outlineplot.YData
outline = [x,z]
z_line = yline(-3)
% Find the indices where the curve crosses the line y = 3
crossing_indices = find(outlineplot == -3)
% intersections = intersect(outline, x_line)
% Count the number of intersections
num_intersections = numel(crossing_indices)
% Display the result
disp(['The curve intersects the line z = -3 at ', num2str(num_intersections), ' point(s).']);
2 个评论
Dyuman Joshi
2023-7-10
编辑:Dyuman Joshi
2023-7-10
You are assuming that the plot data will have a point with y-coordinate equals to -3, which is not the case.
Also, you are comparing a graphics object, a plot (i.e. outlineplot variable), with a numerical value via this line of code, which does not make sense.
crossing_indices = find(outlineplot == -3)
And as a plot is not equal to a number, the output of the comparison is 0 thus find() returns an empty array.
If you want to compare, use the YData from the plot to compare.
There are some FEX files you can look into - Fast and Robust Curve Intersections and Curve Intersections
采纳的回答
Chunru
2023-7-10
rho = 4420; %kg/m^3
Cp = 550; %J/kg?K
T0 = 303.15; %K
A = 0.5; %[Absorbtivity]
k = 7.2; %W/m/K
alpha = 2.96*10^-6; %m^2/s
D = alpha;
P = 500; %W
u = 0.5; %m/s
Tm = 1933; %K
d_laser = 0.0001; %m
r_laser = d_laser/2; %m
a = r_laser;
p = D/(u*a);
%Define
n = 100;
% x = linspace(-0.00015,0.00025,n);
% z = linspace(-0.000005,0,n);
x = linspace(-0.00025,0.0012,n);
z = linspace(-0.0001,0,n);
%Normalized
x_nor = x/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh,z_mesh] = meshgrid(x_nor,z_nor);
y_mesh = 0;
fun = @(t) exp((-z_mesh.^2./(4*t))-((y_mesh.^2+(x_mesh-t).^2)./(4*p.*t+1)))./((4.*p.*t+1).*sqrt(t));
g = integral(fun,0,Inf,'ArrayValued',true);
% Ts = (A*P)/(pi*rho*Cp*sqrt(D*u*(a^3)))
Ts = Tm/0.16713
g2 = Tm/Ts
figure
hold on
[C,h] = contourf(x_mesh,z_mesh,g, 'ShowText', 'On')
[C1,h1] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList=0.16713,Color='red')
[C2,h2] = contour(x_mesh,z_mesh,g, 'ShowText', 'On',LevelList = g2,Color='white')
hold off
clabel(C,h,'Color','white')
clabel(C1,h1,'Color','white')
clabel(C2,h2,'Color','white')
% xticks(-5:5:25)
% yticks(-10:2:0)
% axis([-5 25 -10 0])
title('x\primez\prime plane')
xlabel('x\prime')
ylabel('z\prime')
colorbar
%Extract data from contour plot
contourdata = contourc(x_nor,z_nor,g,[g2 g2])
% Initialize arrays to store the boundary coordinates
x_prime = [contourdata(1,2:size(contourdata,2))]
z_prime = [contourdata(2,2:size(contourdata,2))]
outlineplot = plot(x_prime,z_prime)
x_outline = outlineplot.XData
z_outline = outlineplot.YData
outline = [x,z]
z_line = yline(-3)
% Find the indices where the curve crosses the line y = 3
crossing_indices = find(outlineplot == -3)
% intersections = intersect(outline, x_line)
% Count the number of intersections
num_intersections = numel(crossing_indices)
% Display the result
disp(['The curve intersects the line z = -3 at ', num2str(num_intersections), ' point(s).']);
% Intersection of polygon and line
figure
plot(x_prime, z_prime);
hold on
yline(-3)
p = polyshape(x_prime, z_prime);
l = [-5 -3; 15 -3];
[in out] = intersect(p, l);
in(:, 1) % the intersect point
plot(in(:, 1), [-3 -3], 'ro')
更多回答(1 个)
Sandeep Mishra
2023-7-10
Hello Sarinya,
I understand that you want to find the intersection points of you contour plot for y==-3, But you are using the find operator with plot object, so you are getting an empty array as result.
Also in your code the outlinePlot.YData conatins double data type & it doesn't contain exact -3.0000 value so you will not get your desried result.
You can try the below code to get your desired result.
% Indices of value > -3.01
ind1 = find(outlineplot.YData>-3.01)
% Indices of value < -2.99
ind2 = find(outlineplot.YData<-2.99)
% Indices of value > -3.01 and < -2.99
ind = intersect(ind1, ind2)
% X coordinates
xCoord = outlineplot.XData(ind)
You can refer to the below documentation to learn more about "find" in MATLAB.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!