Need help in plotting the 3 variable function's contour

clc
clear all
close all
% format long
step = 10;
z = linspace(0,1.2e-6,step);
x = linspace(-0.1e-6,0.1e-6,step);
y=linspace(-0.1e-6,0.1e-6,step);
%r=sqrt(x.^2+y.^2);
z1 = 1e-6;
L = 9.995498999999999e-07;
a1 = 8.999974099901931e-10;
s = linspace(-L,L,1e5);
F0 = 3.15e7;
Mat_V = zeros(length(z),length(x),length(y));
for i = 1 : length(z)
    i
    for j = 1 : length(x)
        for k= 1:length(y)
            if abs(sqrt(x(j).^2+y(k).^2)) > sqrt((a1*z1)*(1-(z(i)*z(i)/z1/z1)))
                F_f = s ./((sqrt(x(j).^2+y(k).^2) + (z(i)-s).^2).^0.5); 
                V_f = @(z,F_f) F0 * z1 / (z1 * log((L+z1)/(-L+z1)) - 2 * L) * (trapz(s,F_f)) - F0 * z;
                Mat_V(i,j,k) = V_f(z(i),F_f);
            else 
                Mat_V(i,j,k) = 0;
            end
        end
    end
end
[X,Y,Z] = ndgrid(x,y,z) ;
F =Mat_V(Z,X,Y) ;
figure
hold on
for i = 1:100
    surf(X(:,:,i),Y(:,:,i),Z(:,:,i),F(:,:,i)) ;
end
axis([-3 3 -3 3]);grid minor
view(3);
shading interp
%set(gca,'FontSize',25,'fontweight','bold')
%size = 35;
%xlabel('X ({\mu}m)','FontSize',size,'fontweight','bold')
%ylabel('Z ({\mu}m)','FontSize',size,'fontweight','bold')
%axis equal

[Walter Roberson]

Mat_V(i,j,k) = V_f(z(i),F_f);

That is an array assignment with integer indices.

F =Mat_V(Z,X,Y) ;

That is an attempt to index the array with non-integer indices. Perhaps you want to use permute()?

[Yasir Iqbal]

Can you explain where do I write this permute() in my code??

Thanks.

Answer 1

[Walter Roberson]

format long g

step = 10;
z = linspace(0,1.2e-6,step);
x = linspace(-0.1e-6,0.1e-6,step);
y=linspace(-0.1e-6,0.1e-6,step);
%r=sqrt(x.^2+y.^2);
z1 = 1e-6;
L = 9.995498999999999e-07;
a1 = 8.999974099901931e-10;
s = linspace(-L,L,1e5);
F0 = 3.15e7;
Mat_V = zeros(length(z),length(x),length(y));
for i = 1 : length(z)
    i
    for j = 1 : length(x)
        for k= 1:length(y)
            if abs(sqrt(x(j).^2+y(k).^2)) > sqrt((a1*z1)*(1-(z(i)*z(i)/z1/z1)))
                F_f = s ./((sqrt(x(j).^2+y(k).^2) + (z(i)-s).^2).^0.5); 
                V_f = @(z,F_f) F0 * z1 / (z1 * log((L+z1)/(-L+z1)) - 2 * L) * (trapz(s,F_f)) - F0 * z;
                Mat_V(i,j,k) = V_f(z(i),F_f);
            else 
                Mat_V(i,j,k) = 0;
            end
        end
    end
end

i = 
     1

i = 
     2

i = 
     3

i = 
     4

i = 
     5

i = 
     6

i = 
     7

i = 
     8

i = 
     9

i = 
    10

[X,Y,Z] = ndgrid(x,y,z) ;
F = permute(Mat_V, [2 3 1]);  %created as z x y

figure
hold on
for i = 1:size(F,3)
    surf(X(:,:,i),Y(:,:,i),Z(:,:,i),F(:,:,i)) ;
end
xlim auto; ylim auto; grid minor
view(3);
shading interp

[Yasir Iqbal]

It ran perfectly the first time but when I rerun it, its not going through. I do not know.

[Walter Roberson]

hold off

at the bottom of the code

Note that your original axis limit was -3 to 3 but the range of your variables is much smaller.

Have you considered not using surf for this purpose, and instead use several isosurface() calls with different levels, and making sure to set the Alpha (transparency) level below 1 so that the inner shells can be seen?

Do you have a sample plot that you are trying to be similar to?

[Yasir Iqbal]

I am trying to get something like this but then there were input arguments 2 dimensional so I could use contour() but here I have 3 input variables (x,y,z).

[Walter Roberson]

I modified the code to be slightly more efficient.

You will notice that the outputs are essentially flat. If you look at the term that I write into the temporary variable "temp" and save overall into the array Temp, you will find that that term is at most 6E-14 and so the subtraction of F0*z(i) that you were doing is contributing all of the significant value to the Mat_V result -- to withing round-off error for plotting purposes, your calculations just produce flat sheets at different z levels.

format long g

step = 10;
z = linspace(0,1.2e-6,step);
x = linspace(-0.1e-6,0.1e-6,step);
y=linspace(-0.1e-6,0.1e-6,step);
z1 = 1e-6;
L = 9.995498999999999e-07;
a1 = 8.999974099901931e-10;
s = linspace(-L,L,1e5);
F0 = 3.15e7;
Mat_V = zeros(length(z),length(x),length(y));
Temp = Mat_V;

for i = 1 : length(z)
    zi = z(i);
    for j = 1 : length(x)
        xj = x(j);
        for k= 1:length(y)
            yk = y(k);
            if abs(sqrt(xj.^2+yk.^2)) > sqrt((a1*z1)*(1-(zi*zi/z1/z1)))
                F_f = s ./((sqrt(xj.^2+yk.^2) + (zi-s).^2).^0.5);
                temp = z1 / (z1 * log((L+z1)/(-L+z1)) - 2 * L) * (trapz(s,F_f));
                V_f = F0 * (temp - zi);
                Mat_V(i,j,k) = V_f;
                Temp(i,j,k) = temp;
            else 
                Mat_V(i,j,k) = nan;
            end
        end
    end
end
[X,Y,Z] = meshgrid(x,y,z) ;
F = permute(Mat_V, [2 3 1]);  %created as z x y

fig = figure();
[Fmin, Fmax] = bounds(F, 'all');
Fvec = linspace(Fmin, Fmax, 6);

for K = 1 : length(Fvec)
    isosurface(X, Y, Z, F, Fvec(K));
end

set(findobj(fig, 'type', 'patch'), 'FaceAlpha', 0.5)

[Tmin, Tmax] = bounds(Temp, 'all')

Tmin = 
     -6.46248678947138e-26

Tmax = 
      6.33635455266913e-14

[Walter Roberson]

format long g

step = 10;
z = linspace(0,1.2e-6,step);
x = linspace(-0.1e-6,0.1e-6,step);
y=linspace(-0.1e-6,0.1e-6,step);
z1 = 1e-6;
L = 9.995498999999999e-07;
a1 = 8.999974099901931e-10;
s = linspace(-L,L,1e5);
F0 = 3.15e7;
Mat_V = zeros(length(z),length(x),length(y));
Temp = Mat_V;

for i = 1 : length(z)
    zi = z(i);
    for j = 1 : length(x)
        xj = x(j);
        for k= 1:length(y)
            yk = y(k);
            if abs(sqrt(xj.^2+yk.^2)) > sqrt((a1*z1)*(1-(zi*zi/z1/z1)))
                F_f = s ./((sqrt(xj.^2+yk.^2) + (zi-s).^2).^0.5);
                temp = z1 / (z1 * log((L+z1)/(-L+z1)) - 2 * L) * (trapz(s,F_f));
                V_f = F0 * (temp - zi);
                Mat_V(i,j,k) = V_f;
                Temp(i,j,k) = temp;
            else 
                Mat_V(i,j,k) = nan;
            end
        end
    end
end
[X,Y,Z] = meshgrid(x,y,z) ;
F = permute(Mat_V, [2 3 1]);  %created as z x y

fig = figure();
[Fmin, Fmax] = bounds(F, 'all');
Fvec = linspace(Fmin, Fmax, 6);

for K = 1 : length(Fvec)
    isosurface(X, Y, Z, F, Fvec(K));
end

set(findobj(fig, 'type', 'patch'), 'FaceAlpha', 0.5)
title('colored by MAT_V')

fig2 = figure()

fig2 = 
  Figure (2) with properties:

      Number: 2
        Name: ''
       Color: [1 1 1]
    Position: [671 661 577 433]
       Units: 'pixels'

  Show all properties

[Tmin, Tmax] = bounds(Temp, 'all');
Tvec = linspace(Tmin, Tmax, 5);
for K = 1 : length(Tvec)
    isosurface(X, Y, Z, Temp, Tvec(K));
end

set(findobj(fig2, 'type', 'patch'), 'FaceAlpha', 0.5)
title('colored by MAT_V without -F0*z term')