Calculating closing speed between pair of aircraft

Hi, I want to calculate the closing speed of all pairs of aircraft to determine the conflict rate.

Suppose that I have these trajectory data table for 2 aircraft, how can I calculate the closing speed (in knots) of those 2 aircraft?

I saw this post (https://gamedev.stackexchange.com/q/118162/) about calculating closing speed (the first answer), but it's quite confusing to calculate like that with my data on Matlab.

I understood that calculating distance between two aircraft at each timestamp and divide by the elapsed time (in this case, 4 sec) produces closing time.

Q1) But the recorded time is not exactly same per aircraft, so this 'time discrepancy' should be considered to calculate distance.

Q2) Also, the aircraft altitude changes every timestamp. So, how to calculate closing speed between different altitudes?

If the time discrepancy matter is solved, I'll calculate the horizontal distance between pairs of aircraft using the Haversine formula, and the vertical distance from the altitudes. Then I'll use Pythagoras' theorem to calculate the diagonal, and divide the diagonal distance by the elapsed time to produce closing time. Does this make sense?

Those are my questions.. any helps will be appreciated!

+ I attached a sample data of 2 aircraft, which is a struct type data, including time, lon,lat,speed,heading,altitude.

++ The recorded time is not exactly same per aircraft, so this time discrepancy should be considered to calculate distance.

I'll use 'haversine' to calculate the distance between two aircraft.

(https://kr.mathworks.com/matlabcentral/fileexchange/27785-distance-calculation-using-haversine-formula)

Note: Since it's not allowed to upload the original data publicly, I change those values manually.

[Image Analyst]

Can you give us something to work with, like attach the two tables in a .mat file with the paperclip icon

save('answers.mat', 'table1', 'table2');

or at least MATLAB code to create the tables. I don't want to type all that in manually, and I don't even have your second table. But basically just compute the distance between the aircrafts at the same times. The change in distance between the two divided by the time difference is the closing velocity.

[Chunru]

Since the vector speeds (speed and heading) are given, you can project the speed difference along the line connecting two airplance (using their coordinates) to find out the closing speed.

[aposta95]

@Image Analyst Yeah, I attached the data with struct type.

@Chunru Thanks for your advice!

Answer 1

[Meg Noah]

Is this what you're looking for?

% Given values of positions and velocities of airplanes
t0_s = [12.0*3600.0 + 0.0*60.0 + 2.0; ...
        12.0*3600.0 + 0.0*60.0 + 6.0; ...
        12.0*3600.0 + 0.0*60.0 + 10.0]';
lon0_deg = [121.78 121.29 121.81];
lat0_deg = [38.27 38.27 38.26];
speed_knot = [513 511 510];
heading_deg = [173 173 173];
alt_m = 10e3*ones(size(lon0_deg)); % approximate with planes at the surface 
nplane = numel(lon0_deg);

% convert speed to m/s
speed_m_per_s = 0.514444*speed_knot;

% chose one time to be the snapshot of the field
T_s = t0_s(1);

% find the time difference between plane positions
delta_time_s = T_s - t0_s;

% find the position of each plane at T_s
% http://www.movable-type.co.uk/scripts/latlong.html
% given start point, initial bearing, and distance find the location
% use this to find the position of each plane at the same time
R_m = 6371.23e3 + alt_m; % Earth's radius + altitude of planes
distance_m = speed_m_per_s .* delta_time_s;
latT_deg = asind(sind(lat0_deg).*cos(distance_m/R_m) + ...
    (cosd(lat0_deg).*sin(distance_m/R_m)).*cosd(heading_deg));
lonT_deg = lon0_deg + ...
    atan2d(sind(heading_deg).*sind(distance_m/R_m).*cosd(lat0_deg), ...
    cos(distance_m/R_m)-sind(lat0_deg).*sind(latT_deg));

% find the ECEF coordinates of the planes
ecef_pos = lla2ecef([latT_deg' lonT_deg' alt_m']);

% assume the planes are traveling at a constant altitude
% velocity in ECEF coordinates is:
uNorth = speed_m_per_s.*cosd(heading_deg);
vEast = speed_m_per_s.*sind(heading_deg);
wDown = 0;
[U,V,W] = ned2ecefv(uNorth,vEast,wDown,latT_deg,lonT_deg);
ecef_vel = [U' V' W'];

% find the closing speed for each of the planes relative to plane 1
for iplane = 2:nplane
    deltaP_m = (ecef_pos(iplane,:)-ecef_pos(1,:));
    unit_deltaP = deltaP_m/norm(deltaP_m);
    deltaV_m_per_s = (ecef_vel(iplane,:)-ecef_vel(1,:));
    closing_speed_m_per_s = dot(deltaV_m_per_s,unit_deltaP);
    v_closing_m_per_s = closing_speed_m_per_s*unit_deltaP;
    fprintf(1,'Plane %d to Plane 1 Closing Speed = %f m/s\n',iplane, ...
        closing_speed_m_per_s);
end