Plotting of a surface

Currently I am working on creating ring potential with cosine-modulated depth. Since I know, that modulation is radial invariant I split my hypothetic function on radial and angular components (like variable separation). Radial part defines mu ring with central radius r0, width w and magnitude of potential is some A:

R(r)= -A*exp(-(r-r0).^2/w.^2)

Next I need cosine modulation. I tried just multiplying on cos(atan2(y,x)) but it causes wrong plot - magnitude of a function at line x = 0 goes to some value >> A. I think that happends because multiplication of small and big power exponents (just guessing). Is there any way to define such potential? P.S. I work in cartesian coordinates, so r = sqrt(x.^2+y.^2) in R(r) function.

[Walter Roberson]

You could experiment with symbolic calculations to see if your hypothesis of exponent problems is correct.

[Andrei]

Thanks Walter on another reply. cos(atan2(y, x)) function is specific on x=0 axis. It is not smooth there

[Andrei]

Answer 1

[William Rose]

@Andrei,

My understanding is that you want to create a potential function which is the product of a radial component and an angular component.

Here is the radial part:

[X,Y]=meshgrid(-4:.04:4);
r=(X.^2+Y.^2).^0.5;
A=1; r0=2; w=2;
R=-A*exp(-(r-r0).^2/w^2);
mesh(X,Y,R); 
colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('R'); title('R(r)')

Now compute and plot the angular part, Phi(theta). Phi() is undefined at the origin, because theta is undefined, because atan2(0,0) is undefined.

theta=atan2(X,Y); 
Phi=cos(4*theta);  %4 to get 4 cycles
mesh(X,Y,Phi); 
colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('Phi'); title('\Phi(\theta)')

If you compute Z(x,y)=R(r)*Phi(theta), (where r and theta are functions of x,y, as shown above), you get an ugly surface, that is especially ugly at the origin.

Another idea: multiply Phi(theta) by a function R2(r), which is is zero at the origin and 1 far from the origin. This will suppress the ugly behavior of Phi() at the origin. w2 is a radial scaling factor that determines how rapidly R2 turns on.

w2=0.7; R2=(1-exp(-r/w2)).^2;
mesh(X,Y,R2); xlabel('X'); ylabel('Y'); zlabel('R2'); title('R2(r)')

Next, compute the modulation function A(r,theta)=(1+ a2*R2( r)*Phi(theta)). a2 (0<a2<1) is a scalar that sets the size of the angular modulation.

a2=0.3;
A=1+a2.*R2.*Phi;
mesh(X,Y,A); colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('A'); title('A(r,\theta)')

The plot of A(r,thata) above looks similar to Phi(theta), but it is not discontinuous at the origin, and the scale is different, as the colorbars indicate.

Now compute Z(r,theta)=R(r)*A(r,theta)

Z=R.*A;
mesh(X,Y,Z); colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('Z'); title('Z(r,\theta)')

This may not be quite what you are looking for, but it gives you some ideas.

Good luck.

[William Rose]

@Andrei,

If w is small compared to r0, the potential is approximately 0 at theorigin, and therefore you do not need R2( r) to suppress the discontinuity at the origin.

[X,Y]=meshgrid(-4:.04:4);
r=(X.^2+Y.^2).^0.5;
A=1; r0=3; w=0.5;
R=-A*exp(-(r-r0).^2/w^2);
mesh(X,Y,R); 
colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('R'); title('R(r)')

Now compute and plot the angular part, Phi(theta). Phi() is undefined at the origin, because theta is undefined, because atan2(0,0) is undefined.

theta=atan2(X,Y); 
Phi=cos(6*theta);  % 6 cycles this time, versus 4 in my earlier answer
mesh(X,Y,Phi); 
colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('Phi'); title('\Phi(\theta)')

Next, compute the modulation function A(r,theta)=(1+ a2*Phi(theta)). a2 (0<a2<1) is a scalar that sets the size of the angular modulation. Note that unlike my previous answer, I am not using R2(r). The plot looks just like the plot of Phi(theta), except the scale on the colorbar is different. The colorbar autoscales its range, by default.

a2=0.2;
A=1+a2.*Phi;
mesh(X,Y,A); colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('A'); title('A(r,\theta)')

Now compute Z(r,theta)=R(r)*A(r,theta)

Z=R.*A;
mesh(X,Y,Z); colorbar; view(0,90); axis equal
xlabel('X'); ylabel('Y'); zlabel('Z'); title('Z(r,\theta)')

Here is a view of the surface from underneath:

mesh(X,Y,Z); view(30,-45)
xlabel('X'); ylabel('Y'); zlabel('Z'); title('Z(r,\theta)')

Try it.

[Andrei]

Thanks a lot. You did too much, I am really grateful. I wanted just a hint, but your solution is absolutely right! Thamks!

[William Rose]

You're welcome, @Andrei.