Converting from Cartesian rates to spherical rates.

I have the following data of a track from a radar sensor:

• Time (UTC)
• Position North from the sensor in meters
• Position East from the sensor in meters
• Position Down from the sensor in meters
• Velocity North in meters/second
• Velocity East in meters/second
• Velocity Down in meters/second

I can easily convert the postion into an Azimuth Elevation and Range using cart2sph and back again with sph2cart.

I would like to convert the velocities into Az Rate, El Rate, and Range Rate. I could do this by converting the positions and differentiating the anwser but that would only give me the 'inferred' velocity between detections rather than the actual velocity meausred by the sensor.

%% SET UP AXES
figure('color', 'w')
hold on
xline(0)
yline(0)
plot3([0,0],[0,0],[-5,5],'k')

%% DETECTION AT EAST NORTH UP POSITION WITH ENU VELOCITIES
% PSON
E = 2;
N = 3;
U = -2;

% Vel
E_vel = 1;
N_vel = -2;
U_vel = 1;

%% PLOT DETECTION AND VELOCITY
quiver3(E,N,U,E_vel, N_vel, U_vel ,'AutoScaleFactor',1,'color','r', 'Marker','x')

%% ADD REFERENCE LINES
plot3([0,E],[0,N],[0,U],'color',[0.8,0.8,0.8])
plot3([0,E+E_vel],[0,N+N_vel],[0,U+U_vel],'color',[0.8,0.8,0.8])

%% TIDY UP
xlabel("East")
ylabel("North")
zlabel("Up")

grid minor
box on
axis equal

xlim([-5,5])
ylim([-5,5])
zlim([-5,5])
view(125,25)

view(0,90)

The Az rate is the angle between the two lines when viewing from the top down

Answer 1

[Adam Danz]

To convert the three Cartesian velocity components to spherical velocity values, you can follow these steps.

I assume East is the x-axis, North is the y-axis, and Down the z-axis.

0. Create demo data

NorthVelocity = 2; 
EastVelocity = -4;
DownVelocity = -3; 

1. Compute the magnitude of the velocity vector. This is the radial velocity.

radialVelocity = sqrt(EastVelocity^ + NorthVelocity^2 + DownVelocity^2);

2. Compute the azimuthal angle. This is the azimuthal velocity over the same unit of time.

azimuthVelocity = atan2(NorthVelocity, EastVelocity);

3. Compute the elevation angle. This is the elevational velocity over the same unit of time.

elevationVelocity = atan2(DownVelocity, sqrt(EastVelocity^2 + NorthVelocity^2));

4. Convert radians to degrees if needed: rad2deg

fprintf('Azimuth: %.1f\nElevation: %.1f\nRadius: %.1f', ...
    azimuthVelocity*180/pi, elevationVelocity*180/pi, radialVelocity);

Azimuth: 153.4
Elevation: -33.9
Radius: 16.3

Plot results

1. Show the original vector and its component vectors

figure()
% Velocity Vector
quiver3(0,0,0,EastVelocity,NorthVelocity,DownVelocity,'off','-ob','LineWidth',2)
hold on
axis equal padded
% Velocity components
quiver3(0,NorthVelocity,DownVelocity,EastVelocity,0,0,'off','--om')
quiver3(EastVelocity,0,DownVelocity,0,NorthVelocity,0,'off','--om')
quiver3(EastVelocity,NorthVelocity,0,0,0,DownVelocity,'off','--om')
xlabel('x (East)')
ylabel('y (North)')
zlabel('z (Down)')
box on
set(gca,'boxstyle','full')
constantplane([1 0 0],0)

2. Ideally we could plot the spherical results in a spherical graphics system but MATLAB doesn't have that. To visually verify the converted components, we could convert them back into Cartesian values and plot another vector but some people may rightfully be suspicious of that approach. Here, I created a wedge to show the azimuthal component and then ran out of time 😀.

% distance from (0,0) to end point
endDist = hypot(NorthVelocity, EastVelocity);

% Azimuth 
thAz = linspace(0,azimuthVelocity);
xAz = endDist * [cos(thAz),0,1];
yAz = endDist * [sin(thAz),0,0]; 
zAz = DownVelocity * ones(size(xAz)); 
patch(xAz, yAz, zAz, [0 .5 0],'FaceAlpha',0.3,'EdgeColor','none')

●

[George Hackney]

Hi Adam, Not sure I worded the question perfectly. you have provided the spehical angles of the track from with the track at the origin, I require the sensor to be at the origin. I have added an example to the question which might help.

[Adam Danz]

It seemed like you intended to focus solely on the velocity data, making the starting points arbitrary. I now believe the primary goal is to compute the offset of the coordinate in spherical coordinates (azimuth, elevation) by using the velocity component data. The starting position is also required.

Start off with your demo

%% SET UP AXES
figure('color', 'w')
hold on
xline(0)
yline(0)
plot3([0,0],[0,0],[-5,5],'k')

%% DETECTION AT EAST NORTH UP POSITION WITH ENU VELOCITIES
% PSON
E = 2;
N = 3;
U = -2;

% Vel
E_vel = 1;
N_vel = -2;
U_vel = 1;

%% PLOT DETECTION AND VELOCITY
quiver3(E,N,U,E_vel, N_vel, U_vel, 'off','color','r', 'Marker','x')

%% ADD REFERENCE LINES
plot3([0,E],[0,N],[0,U],'color',[0.8,0.8,0.8])
plot3([0,E+E_vel],[0,N+N_vel],[0,U+U_vel],'color',[0.8,0.8,0.8])

%% TIDY UP
xlabel("East")
ylabel("North")
zlabel("Up")

grid minor
box on
axis equal

xlim([-5,5])
ylim([-5,5])
zlim([-5,5])
view(125,25)

Compute end point rather than using the stored coordinates.

t = 1;   % time duration
E1 = E + E_vel*t; 
N1 = N + N_vel*t;
U1 = U + U_vel*t;

Plot the end point for verification

plot3(E1,N1,U1,'b*')

Compute the azimuth and elevation angles for starting point and the computed end point. Negative values are moving toward the origin.

az0 = atan2(N,E);
el0 = atan2(U,sqrt(E^2 + N^2));

az1 = atan2(N1,E1);
el1 = atan2(U1,sqrt(E1^2 + N1^2));

Compute angular velocity between the pair of angles and convert to degrees/time

azVel = (az1-az0)/t * 180/pi

azVel = -37.8750

elVel = (el1-el0)/t * 180/pi

elVel = 11.4687

The depth component in spherical coordinates is the distance from (0,0,0). Compute the distance for the starting point and converted end point.

d0 = sqrt(E^2 + N^2 + U^2);
d1 = sqrt(E1^2 + N1^2 + U1^2);
dVel = (d1-d0)/t

dVel = -0.8065