Clear Filters
Clear Filters

How to plot xy, yz and xz plane contour with integration matrix equation

5 views (last 30 days)
I try to follw mt professor to plot contour plot of Rosenthal's equation, and here it my code
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 = 100; %W
v = 1; %m/s
u = v;
Tm = 1933; %K
d_laser = 0.01; %mm
r_laser = d_laser/2; %mm
a = r_laser;
p = D/(u*a);
%Define
x = linspace(-0.005,0.025,100);
y = linspace(-0.0005,0.0005,100);
z = linspace(0,0.005,100);
%Normalized
x_nor = x/a;
y_nor = y/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh, y_mesh] = meshgrid(x_nor, y_nor);
[x_mesh, z_mesh] = meshgrid(x_nor, z_nor);
%Calculation
r = (x_mesh.^2 + y_mesh.^2).^0.5;
Ts = A*P/(pi*rho*Cp*sqrt(D*u*a));
T2 = T0 + (A*P)./(2*pi*k*r).*exp(v.*(r+x_mesh)./(2*D));
q = (2*A*P/(pi*a^2))*exp(-2*r.^2/a^2);
syms t
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 = int(fun,t,[0 Inf]);
%Plot x'y' plane
figure
contourf(x_mesh, y_mesh, g, [303.15:100:2000])
colorbar
title('x\primey\prime plane')
xlabel('x\prime (m)')
ylabel('y\prime (m)')
xlim([-0.005 0.025])
%Plot y'z' plane
figure
contourf(x_mesh, z_mesh, g, [303.15:100:2000])
title('x\primez\prime plane')
xlabel('x\prime (m)')
ylabel('z\prime (m)')
xlim([-0.005 0.025])
  4 Comments
Torsten
Torsten on 26 Jun 2023
Edited: Torsten on 26 Jun 2023
So g should be a three-dimensional array where the first dimension refers to x, the second dimension refers to y and the third dimension refers to z ? Then, to plot in the xy plane, e.g., you want to choose a value "z(iz)" and plot g(:,:,iz) against the x- and y-array ? That's not what you do in your code.
Pannawat Pongsriassawin
Pannawat Pongsriassawin on 26 Jun 2023
yeah, actualy output (g) should be 100x100x100, but it is g(:,:,1), g(:,:,2), ..., g(:,:,100) then it can plot it into 2D plane contour. I don't know how to solve it. I think g(:,:,1) mean g(1:100,1:100,1) then it is g in z=1 plane? If it is, that what I desired.Actual it can't plot or Should I use aother command to plot it?

Sign in to comment.

Accepted Answer

Torsten
Torsten on 26 Jun 2023
Edited: Torsten on 26 Jun 2023
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 = 100; %W
v = 1; %m/s
u = v;
Tm = 1933; %K
d_laser = 0.01; %mm
r_laser = d_laser/2; %mm
a = r_laser;
p = D/(u*a);
%Define
x = linspace(-0.005,0.025,100);
y = linspace(-0.0005,0.0005,100);
z = linspace(0,0.005,100);
%Normalized
x_nor = x/a;
y_nor = y/a;
z_nor = z/(D*a/u).^0.5;
[x_mesh,y_mesh,z_mesh] = ndgrid(x_nor,y_nor,z_nor);
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);
figure(1)
iz = 1;
contourf(y_nor,x_nor,squeeze(g(:,:,iz)))
colorbar
figure(2)
iy = 1;
contourf(z_nor,x_nor,squeeze(g(:,iy,:)))
colorbar
figure(3)
ix = 1;
contourf(z_nor,y_nor,squeeze(g(ix,:,:)))
colorbar

More Answers (1)

Sarthak
Sarthak on 27 Jun 2023
Hello Pannawat,
You can plot the contour by using the meshgrid and surf functions which are available in MATLAB.
[X,Y] = meshgrid(x,y); % Create a grid from x and y vectors
z = zeros(size(x, 1)); % Generate z data
surf(x, y, z); % Plot the contour
You can also refer to the MATLAB Answer provided below for further information.
Hope this helps!!

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!