How to find intersections in contour plot?
18 ビュー (過去 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!