# Why contourf and fcontour gives different images for the same matlab code?

2 views (last 30 days)
Athira T Das on 27 May 2023
Edited: VBBV on 29 May 2023
I have code for fcountour as;
syms x y
m =1;
sgm =1;
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2;
z=30;
Pho=0.003;
for u=1:length(z)
I=0;
Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^4));
C = (1./(4.*(lambda.^2).*(z(u).^4))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m);
fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
X = 0;
for l = 0:m
L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
sum_j = 0;
for j = 0:m
J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);
sum1 = 0.0;
for p = 0:1/2:l/2
faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
for s1 = 0:l-2*p
faktor1_s1 = nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
for q = 0:1/2:(m-l)/2
faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
sum_inner_1 = 0;
for s2 = 0:m-l-2*q
E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
end
end
end
sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
end
sum_j= sum_j + J.*(sum1);
end
X = X + (L.*sum_j);
end
X = (C.*X);
I = I + abs(vpa(X)) ;
figure;
fcontour(real(I), [-1 1 -1 1]*1E-1, 'Fill','on', 'MeshDensity',150);
colormap('hot');
axis('equal');colorbar('vert');
end
Then i have contourf plot as;
figure;
m =1;
sgm =1;
wo = 0.1;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2;
z=30
z = 30
Pho=0.003
Pho = 0.0030
[x, y] = meshgrid(-1:0.06:1);
for u=1:length(z)
I=0;
Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^4));
C = (1./(4.*(lambda.^2).*(z(u).^4))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m);
fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
X = 0;
for l = 0:m
L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
sum_j = 0;
for j = 0:m
J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);
sum1 = 0.0;
for p = 0:1/2:l/2
faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
for s1 = 0:l-2*p
faktor1_s1 = nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
for q = 0:1/2:(m-l)/2
faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
sum_inner_1 = 0;
for s2 = 0:m-l-2*q
E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
end
end
end
sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
end
sum_j= sum_j + J.*(sum1 );
end
X = X + (L.*sum_j);
end
X = (C.*X);
I = I + abs(vpa(X)) ;
contourf(x,y,real(abs(I))/max(real(abs(I(:)))))
colormap('hot');
colorbar('vert');
axis('equal')
end
According to me fcountour plot is correct, then what is wrong with countourf plot? Can anyone help me to fix the error in contourf plot?
VBBV on 28 May 2023
Edited: VBBV on 29 May 2023
What I got is plot with two extra lobes and the figure is not complete. I have done the contourf plot is to normalize my data. Why when I normalize the data, looks different?
The figures produced by both contourf & fcontour functions are same. If you look at fcontour the range limits
fcontour(real(I), [-1 1 -1 1]*1E-1, 'Fill','on', 'MeshDensity',150);
are [ -1 1 -1 1 ] multiplied by (1/10) i.e. 0.1 and fcontour doesnt normalize data when it plots. However, when you normalize the data ( as in case of contourf) it divides with maximum value from the computed vector. Also, the absoulute value was taken, this also means, the resultant value (which includes complex or imaginary values if any) is considered instead of actual value,

VBBV on 28 May 2023
If you change the meshgrid limits in the below line
[x, y] = meshgrid(-1/10:0.005:1/10);
and de-normalize the I variable in the contourf function , you will get similar figure as fcontour
contourf(x,y,real(I))
##### 2 CommentsShow 1 older commentHide 1 older comment
VBBV on 28 May 2023
m =1;
sgm =1;
% this parameter is also something different in both codes
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
Omega = .2;
z=30
z = 30
Pho=0.003
Pho = 0.0030
[x, y] = meshgrid(-1/10:0.005:1/10);
for u=1:length(z)
I=0;
Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
Gamma_star = subs(Gamma, 1i, -1i);
alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^4));
C = (1./(4.*(lambda.^2).*(z(u).^4))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m);
fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
X = 0;
for l = 0:m
L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
sum_j = 0;
for j = 0:m
J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);
sum1 = 0.0;
for p = 0:1/2:l/2
faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
for s1 = 0:l-2*p
faktor1_s1 = nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
for q = 0:1/2:(m-l)/2
faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
sum_inner_1 = 0;
for s2 = 0:m-l-2*q
E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
end
end
end
sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
end
sum_j= sum_j + J.*(sum1 );
end
X = X + (L.*sum_j);
end
X = (C.*X);
I = I + abs(vpa(X)) ;
contourf(x,y,real((I)), 'EdgeColor','none')
colormap('hot');
colorbar('vert');
axis('equal')
end
In addition to x & y range limits in meshgrid and denormailized I variable, the value of parameter wo is used differently in contourf & fcontour functions , See the below line in both codes
% this value is 0.01 in fcontour and 0.1 in contourf function codes, why ?
wo = 0.01;