Meshing a complex function with limits

Hi People

Completely new to matlab coding so will try my best to explain my problems with this code so far:

I am trying to create a mesh in matlab that will reflect the complex function of a optical lens. The code i have so far is this:

%Matlab coding - OTF of ideal lens

%define x and y boundaries, z will be the function

x=(-1:0.03:1);

y=(-1:0.03:1);

%create meshgrid

[X,Y] = meshgrid(x,y);

%define the frequency

rho = (X.^2+Y.^2);

rho0 = 8;

v = rho./2*rho0;

syms z

z=(acos(v) - (v).*sqrt(1-(v).^2));

%surf plot with OTF as Z axis - only interested in real values of z

surf(x,y,real(z));

This code is to reflect the equation in the photo H(rho), the big equation:

My problem is the following:

The wrong shape is simulated, as it has to be the following shape in this final photo:

I know i need to include the limits, but I've no idea how to code this in, whenever i try nd use the "limit" function matlab throws an error saying it can't use the limit function for variable type double.

[Haider Anjum]

if you change rho to 0.25 it will give roughly the shape but again it's not the spread out conical shape I require.

Super confused....If you check the photo from the textbook the 2-D line, the x axis is rho/2rho0 and the y axis is the function H...

Answer 1

[David Goodmanson]

Hello Haider, for someone who is new to Matlab you are doing well. You just need a couple of changes.

First, you forgot that rho = sqrt(x.^2+y.^2) and that makes all the difference.

The statement v = rho./2*rho0 actually multiplies by rho0 and needs to be rho/(2*rho0) instead. [the dot in ./ doesn't hurt anything but it does not have to be there because you are dividing by a scalar quantity].

You are calculating z numerically and there is no need to make z a sym, that just slows things down.

rho <= 2* rho0 is equivalent to v <= 1. The meshgrid limits are large enough that so far v does have v > 1 values, in which case z becomes complex. The 'include the limits' requirement you were looking for is the command v(v>1) = 1. Just like it looks, this finds all the matrix indices where v>1 and sets v = 1 in those cases. Quick, easy, clear. No wonder people like Matlab.

I inserted the (2/pi) and changed the colormap to something different, since the default 'parula' colormap didn't look very good in this case.

%define x and y boundaries, z will be the function
xmax = 22;
ymax = 22;
x = (-xmax:1:ymax);
y = (-ymax:1:ymax);
%create meshgrid
[X,Y] = meshgrid(x,y); 
%define the frequency
rho = sqrt((X.^2+Y.^2));
rho0 = 8;
v = rho./(2*rho0);  
v(v>1) = 1;
z = (2/pi)*(acos(v) - (v).*sqrt(1-(v).^2));
%surf plot with OTF as Z axis
colormap('cool')
surf(x,y,real(z));
colorbar

[Haider Anjum]

Has anyone ever told you, you were their hero?

Because you're mine! Thanks!

The biggest problem I had was setting the limit for v, who knew such a simple line of code was the answer. again thanks!