Trouble using atan2 for coordinate transform, maybe recommend other tools.

I am trying to use atan2 on a 2D plane (X,Z) centered at (0,0) to define an angle phi.

I know that using atan(Z./X) or atan2(Z,X) will yield different discontinuities at either pi/2 or pi respectively. I need this angle to define a coordinate transform matrix for a vector V(X,Y,Z) into V(t,y,r).

This basically what I am trying to do, and I'm not sure

v1 = V(1,i,j,k);    % [i,j,k] represent indices within the X,Y,Z position and V(1,i,j,k) is same size.
v2 = V(2,i,j,k);    % second component
v3 = V(3,i,j,k);    % third component

The transformation matrix is defined as follows...

phi = atan2(Z(i,j,k),X(i,j,k));   % or phi = atan(Z(i,j,k)./X(i,j,k));
J = [sin(phi) 0 -cos(phi);
    0 1 0;
    cos(phi) 0 sin(phi)];
v_prime = J*[v1,v2,v3]';

... but does not take into account the discontinuity in phi. As a result, I get a similar discontinuity in v_prime, but this should not be there realistically.

Is there a way that I can avoid the discontinuity from phi showing up in v_prime? Maybe I can use a different tool that will give phi from 0 to 2*pi instead? Technically phi should be increasing infinitely but I only care about the continuous range from 0 to 2*pi. I'm not aware of any particular functions that will do what I'm looking for.

Note:

Even though I am using R2021a for these calculations, I need any possible solution to also be compatible with R2016b since that is what we use in our lab.

[David Goodmanson]

HI Nathaniel,

if you use atan2 and immediately feed that result into sine or cosine, there can't be any discontinuities because sine and cosine are immune to 2*pi discontinuities in angle. The same can't be said for atan. So the answer may lie elsewhere. For example one posssibility is , when theta = 0 one might expect the rotation to be the identity matrix and not do any rotation. but the J matrix for theta = 0 is a 90 degree rotation.

[Nathaniel H Werner]

Hi David,

Maybe you can go into more detail. If I understood you correctly, sin(atan2(Z/X)) should not exhibit any discontinuities, but in fact, it does when I multiply J*v there is a discontinuity in v_prime that I believe lines up with where the discontinuity in atan2 is. I guess either way I could try defining the angle with asin or acos instead.

The other thing you need to know, and maybe I wasn't clear enough in question, is that J is not a rotation matrix. It is a transformation of coordinates. I have a Cartesian vector in X,Y,Z, and v_prime is another vector with the same norm now X and Z are replaced with a tangential and radial component. It's basically a polar coordinate system but designed to have the properties of a Cartesian one. The explanation is more of a distraction for this question.

[John D'Errico]

@Nathaniel H Werner

Um sorry, but J IS a rotation matrix. It is also a transformation of coordinates, but the transformation it performs is to effecively rotate the coodinate system, by an angle (pi/2-phi). The rotation is around one of the axes. The location of the 1 indicates the rotation is performed around the y axis, leaving Y unchanged, while X and Z spin.

Why do I say the rotation is by the angle pi/2-phi? Because a simple rotation matrix that will leave everything unchanged is the identity matrix. And you will get that for J when phi = pi/2 radians.

Call it what you want, but one of the things you COULD call it is a rotation matrix, and that is arguably the best description. It is the best description, because a general transformation of coordinates could elongate or compact one of the dimensions. It could make the original axes no longer perpendicular to each other. And this matrix does not do any of that. So a rotation matrix is most descriptive.

Disregarding all of that, I'm also not sure what you mean by a discontinuity in phi. Let me see if I can write an answer that may possibly clear things up.

[David Goodmanson]

Hi Nathaniel,

I wanted to make the point that if Z, X and all of the v_i are smooth functions (no discontinuities) of whatever variable is being varied, then although theta = atan2(Z,X) can have discontinuities, cos(theta) and sin(theta) cannot, meaning that J cannot and neither will J*[v1,v2,v3]. If any of Z, X or v_i have discontinuities then of course it's a different story.

There are exceptions. For example if either Z or X goes through zero and has a bit of noise, then the angle can vary wildly near the zero. Another happens if Z and X are smooth but the path goes right through the origin in the X_Z plane. For example

u = -1:.001:1;
X = u;
Z = u.^2;
theta = atan2(Z,X);

has discontinuities in cos(theta) and sin(theta).

[Nathaniel H Werner]

Hi David,

Thanks I think that makes sense.

X and Z that do cross over zero. X ranges from -30 to 40, and Z ranges from -60 to 60. So I guess that could be resulting in the discontinuities in v_prime. I wonder if it could also be due to calculating the angle ahead of time. Maybe I could try to do it in the same for loop as the multiplication of J*[v1,v2,v3].

[David Goodmanson]

Hi Nathanial,

For the type of discontinuity I was talking about both Z and X have to cross the origin 'at the same time'. As long as the appropriate index is kept track of correctly, I doubt that the issue comes from calculating beforehand, but it's worth a try.

Answer 1

[John D'Errico]

The transformation you are using is a simple tranformation into polar coordinates. When you have y in there also, this becomes what is commonly known as cylindrical coordinates. So we will have a transformation of your coordinate system from (X,Y,Z) into (THETA,Y,R). Totally common in mathematics, and often very useful.

A useful tool in this is to just use cart2pol. It will convert the pair (X,Z) into (theta,R). And since you know that Y is unchanged, that will make things almost trivial, simpler even than using atan2 and then computing R.

A problem is, you seem to think there is a problem of a discontinuity. In reality, yes, we will expect to see a discontinuity, in the sense the angle theta will always lie in the (half open) interval [-pi,pi). Or [-180,+180), if you prefer to think in degrees.

That is, you can imagine the rotation being done smoothly in a circle. At some point, the rotation angle starts to approach 180 degrees. Just at the point where it tips over, past 180 degrees say to 181 degrees, the angle reported will revert to -179 degrees. The rotation is the same. Everything is happy. But because the angles as reported only live in the interval [-180,180), you must see what you are describing as a discontinuity, because at that point, you are visualizing this as an angle in a polar coordinate system. You can still have angles in that coordinate system that are larger than 180, or less than -180. But when you compute that rotation angle from X and Z, it will be contained in that limited interval. And that is what you are thinking of as a discontinuity.

(Another sometimes useful tool here is the unwrap function. You really don't need it in this context though.)

But, really, is this a discontinuity in the real world? NO! The rotation matrix you have does not truly convert things into cylindrical coordinates at all. You can think of it as an implicit transformation. All it does is to rotate the X and Z axes around the Y axis, so leaving the Y axis unchanged. If you think of the result after rotation, there is no discontinuity. If you spin a top around, there is no huge bump as it spins. You do not hear a thump, thump, thump,... as it spins. Just a faint whir as it disturbs the air. The discontinuity lies only in your mind as you try to visualize the cylindrical coordinate system.

Hmm. Maybe I need to expain this as an animation. That may take some time to draw pictures. So first, I'll try to draw a mental picture:

Consider a point that rotates around a circle. So imagine this as the point (1,Y,0). Think of that point rotating around the y axis. At a rotation of 90 degrees, it will become (0,Y,1), etc. Eventiually, it gets to the point at 180 degrees, where the coordinates after rotation are now (-1,Y,0). If we rotate just a bit more, just 1 degree more, things are still moving smoothly. There is no huge bump in the night. As I have written it here, as a function of the rotation angle phi, that original point will have rotated to: (cos(phi),Y,sin(phi)). Would you agree with that?

In fact, for ANY smoothly incrementing angle phi, that rotation will be smooth and continuous. Again, no bumps will ever be seen, none felt. As you smoothly change phi, the transformation is smooth and continuous. The only discontinuity you will ever see is when you visualize the transformation in terms of that polar coordinte system. And this is because the rotation angle will always be reported in a finite range like [-pi,pi). So the discontinuity lies only in the mind, not in reality.

I'll hope for now this is sufficient, because it is Sunday morning here and I am feeling far, far too lazy to draw anything more than mental pictures and thought experiments. :)

[Nathaniel H Werner]

It is not cylindrical because each of the components retains the same units of length. J has similarities to rotation matrices, but since it's not acting to rotate anything, this semantic argument is rather fruitless. So no, I don't agree with your illustration at all because each particle at any specific grid point is not undergoing any rotation. It remains stationary and the components effectively become aligned with a cylindrical coordinate system. That is the primary similarity between the two. You can do the derivation yourself of the Jacobian and see that this matrix is not identical to cylindrical, but I really don't see the need to go into that much detail here.

I have no problem with expecting the discontinuity mathematically, however my data is a continuous spatially varying vector valued function. Any discontinuity is an indication of non-physicality and is therefore not permitted.

As I have mentioned in a previous comment, the discontinuities are not in my mind. MatLab is actually plotting the transformed data in 3d space with a discontinuity. Hence, my calculation of the angle needed to perform the transformation is causing the problem. So I need a new method.

Anyway, since I am primarily interested in the radial coordinate (sorry if I didn't mention that in the question, I was try to be brief), maybe I can use the built in tools for a cylindrical transformation. It had not occurred to me, but that's probably a good idea. I'll try that.

When I get back to my computer, I'll update my question with some figures showing the discontinuity in my attempt at calculating the radial component.

[Nathaniel H Werner]

Apparently I left the drive with my data in the lab yesterday, so I can't update the question to show the figures. Sorry about that.