use data from a vector to locate a row in a matrix and extract the data from the row.

So I have a 902x1 vector of unique time stamps, and I have 9 different 901x6 arrays of data (column 1 is the time stamp for each row of data)

I need to take the time stamp from the unique_time vector and find the timestamp in each 901x6 array that matches it and extract those rows of data.

i've tried for loops and if statements, but none work, I have to admit this level of complexity is above my head (programatically speaking).

The part of the code in question starts at line 70, there are a lot of functions involved in the script as a whole so i'll add them as well.

In my mind it should go something like

for i = 1:unique_time
    for j = 1:9
        if unique_time(i) == mat(j)(:, 1)
            Data_mat = [mat(j)(:, 2); mat(j)(:, 3); mat(j)(:, 4)];
        else 
            Data_mat = [0, 0, 0];
        end
    end
end

But this fails miserably. I'm sure it's some trick with the ismember function that I'm just not experienced enough to figure out.

Any help/advise/point in a general direction would be greatly appreciated.

clear all; close all; clc;

%% Setup WGS84 parameters

WGS.EQUATORIAL_RADIUS = 6378137;
WGS.FLATTENING = 1/298.257223563;
WGS.EARTH_ROTATION = 7292115 * 10^(-11); % rad/s
WGS.EARTH_GRAVITIONAL_CONSTANT = 3.9860050 * 10^(14); % m^3/s^2
WGS.ECCENTRICITY = sqrt(WGS.FLATTENING * (2-WGS.FLATTENING));
WGS.ECCENTRICITY_SQUARED = WGS.FLATTENING * (2-WGS.FLATTENING);
WGS.POLAR_RADIUS = WGS.EQUATORIAL_RADIUS * (1-WGS.FLATTENING);

%% Compute average base position

fileID = fopen('data_base\ecef_rx0.txt');
base = textscan(fileID, '%f %d %f %f %f %f %f %f %f %f %f %f %f %f %d %d %d %d');
fclose(fileID);

X = mean(base{1, 3});
Y = mean(base{1, 4});
Z = mean(base{1, 5});
xyz = [X, Y, Z];
base_lla = ecef2lla(xyz, WGS);

%% Compute Sattelite Line of Sight Vectors

% Read in icp satelite data
sat = string({'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat2.txt'; 
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat4.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat5.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat9.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat10.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat12.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat17.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat23.txt';
       'E:\AEM Grad\AEM 691\Project2\project2_resources\data_base\icp_sat25.txt'});

% Read in broadcast file
ephemeris = read_rinex_nav('brdc2930.11n');

% Assign each satelite a number
sat_num = 1:1:size(sat, 1);   

% For Loopy stuff to find Line of Sight Vectors
for i = 1:1:size(sat, 1)
    fileID = fopen(sat(i), 'r');
    current_sat = sat_num(i);
    base_obs = cell2mat(textscan(fileID, '%f %f %f %f %f %f %f %f'));
    eph_rel = ephemeris(find(ephemeris(:, 1) == current_sat), :);
    los_vec = [];
    for j = 1:1:size(base_obs, 1)
        [~, eph_rel_line] = min(abs(base_obs(j, 1) - eph_rel(:, 17)));
        ecef_pos = gpseph2ecef(base_obs(j, 1), eph_rel(eph_rel_line, :));
        sv_obs = [0, 0, 0, 0, 0, 0];
        sv_obs(1) = round(base_obs(j, 1));
        ned = ecef2ned([ecef_pos(1), ...
            ecef_pos(2), ecef_pos(3)], xyz, WGS);
        sv_obs(2) = ned(1);
        sv_obs(3) = ned(2);
        sv_obs(4) = ned(3);
        sv_obs(5) = -asin(sv_obs(4)/sqrt(sv_obs(2)^2 + sv_obs(3)^2 ...
            + sv_obs(4)^2))*180/pi;
        sv_obs(6) = base_obs(j, 8);
        los_vec = vertcat(los_vec, sv_obs);
    end
    base_data{i} = los_vec;
    fclose(fileID);
end

%% Combine 9 different size arrays into 1 single array 
%% with zeroes filling the extra values.
mat1 = base_data{1, 1};
mat1(901, 6) = 0;
mat2 = base_data{1, 2};
mat2(901, 6) = 0;
mat3 = base_data{1, 3};
mat3(901, 6) = 0;
mat4 = base_data{1, 4};
mat4(901, 6) = 0;
mat5 = base_data{1, 5};
mat5(901, 6) = 0;
mat6 = base_data{1, 6};
mat6(901, 6) = 0;
mat7 = base_data{1, 7};
mat7(901, 6) = 0;
mat8 = base_data{1, 8};
mat8(901, 6) = 0;
mat9 = base_data{1, 9};
mat9(901, 6) = 0;

%% Extract X, Y, Z coordinate columns into individual arrays
X_mat = [mat1(:, 2), mat2(:, 2), mat3(:, 2), mat4(:, 2), mat5(:, 2)...
    mat6(:, 2), mat7(:, 2), mat8(:, 2), mat9(:, 2)];
Y_mat = [mat1(:, 3), mat2(:, 3), mat3(:, 3), mat4(:, 3), mat5(:, 3)...
    mat6(:, 3), mat7(:, 3), mat8(:, 3), mat9(:, 3)];
Z_mat = [mat1(:, 4), mat2(:, 4), mat3(:, 4), mat4(:, 4), mat5(:, 4)...
    mat6(:, 4), mat7(:, 4), mat8(:, 4), mat9(:, 4)];

%% Find each unique time step
time_mat = [mat1(:, 1), mat2(:, 1), mat3(:, 1), mat4(:, 1), mat5(:, 1)...
    mat6(:, 1), mat7(:, 1), mat8(:, 1), mat9(:, 1)];
time_linear = time_mat(:);
time_sort = sort(time_linear);
time_unique = unique(time_sort);

%% Find highest value for X, Y, and Z matrices at each unique time step

Functions used in script:

function lla = ecef2lla(xyz, WGS)
% Converts the WGS Earth-Centered, Earth-Fixed (ECEF) Cartesian 
% coordinates to geodetic latitude, longitude, altitude (LLA).
    MU = 3.986005 * 10^14;  % m^3/s^2
    SPEED_OF_LIGHT = 2.99792458 * 10^8;  % m/s
    EARTH_ROT_RATE = 7.2921151467 * 10^-5;  % rad/s
    PI = 3.1415926535898;
    FLATTENING = 1 / 298.257223563;
    ECCENTRICITY = sqrt(FLATTENING * (2 - FLATTENING));
    EQ_RAD = 6378137;                      % m
    POL_RAD = EQ_RAD * (1 - FLATTENING);
    WGS.EQUATORIAL_RADIUS = 6378137;
    WGS.FLATTENING = 1/298.257223563;
    WGS.EARTH_ROTATION = 7292115 * 10^(-11); % rad/s
    WGS.EARTH_GRAVITIONAL_CONSTANT = 3.9860050 * 10^(14); % m^3/s^2

    WGS.ECCENTRICITY = sqrt(WGS.FLATTENING * (2-WGS.FLATTENING));
    WGS.ECCENTRICITY_SQUARED = WGS.FLATTENING * (2-WGS.FLATTENING);
    WGS.POLAR_RADIUS = WGS.EQUATORIAL_RADIUS * (1-WGS.FLATTENING);

	RAD2DEG = 180/pi;
    c = 2.99792458 * 10^8;  % m/s
    x = xyz(1);
    y = xyz(2);
    z = xyz(3);
    
    % Calculate longitude.
	if ((xyz(1) == 0.0) && (xyz(2) == 0.0))

		long = 0.0;

    else

		long = atan2(xyz(2), xyz(1))*RAD2DEG;

    end

    % Compute altitude and latitude.
	if ((xyz(1) == 0.0) && (xyz(2) == 0.0) && ...
            (xyz(3) == 0.0))
        error('XYZ at center of earth');

    else 
        
        % Compute altitude.
		p = norm([xyz(1) xyz(2)]);
        E = sqrt(WGS.EQUATORIAL_RADIUS^2 - WGS.POLAR_RADIUS^2);
        F = 54 * (WGS.POLAR_RADIUS*xyz(3))^2;
        G = p^2 + (1-WGS.ECCENTRICITY_SQUARED)*xyz(3)^2 ...
            - (WGS.ECCENTRICITY*E)^2;
        s = (1 + WGS.ECCENTRICITY_SQUARED^2*F*p^2/G^3 ...
            + sqrt(c^2 + 2*c))^(1/3);
        P = (F/(3*G^2)) / ((s + (1/s) + 1)^2);
        Q = sqrt(1 + 2*WGS.ECCENTRICITY_SQUARED^2*P);

        k_1 = -P*WGS.ECCENTRICITY_SQUARED*p / (1 + Q);
        k_2 = 0.5*WGS.EQUATORIAL_RADIUS^2*(1 + 1/Q);
        k_3 = -(1 - WGS.ECCENTRICITY_SQUARED)*P*z^2 / (Q*(1 + Q));
        k_4 = -0.5*P*p^2;
        r_0 = k_1 + sqrt(k_2 + k_3 + k_4);
        k_5 = (p - WGS.ECCENTRICITY_SQUARED*r_0);

        U = sqrt(k_5^2 + xyz(1)^2);
        V = sqrt(k_5^2 + (1 - WGS.ECCENTRICITY_SQUARED)*xyz(3)^2);
        
        alt = U*(1 - (WGS.POLAR_RADIUS^2 ...
            /(WGS.EQUATORIAL_RADIUS*V)));
        
        % Compute latitude.
        
        z_0 = (WGS.POLAR_RADIUS^2*xyz(3)) ...
            / (WGS.EQUATORIAL_RADIUS*V);
        e_p = (WGS.EQUATORIAL_RADIUS/WGS.POLAR_RADIUS) ...
            * WGS.ECCENTRICITY;
        
        lat = atan((z + z_0*(e_p)^2)/p) * RAD2DEG;

	end

    lla = [lat; long; alt];

end

function ephemeris = read_rinex_nav( filename )
% Read the RINEX navigation file.  The header is skipped because
% information in it (a0, a1, iono alpha and beta parameters) is not 
% currently needed for orbit computation.  This can be easily amended to
% include navigation header information by adding lines in the 'while' loop
% where the header is currently skipped.
% 
% Input:        - filename - enter the filename to be read.  If filename
%                            exists, the orbit will be calculated.
% 
% Output:       - ephemeris - Output is a matrix with rows for each PRN and
%                             columns as follows:
% 
%                  col  1:    PRN    ....... satellite PRN          
%                  col  2:    M0     ....... mean anomaly at reference time
%                  col  3:    delta_n  ..... mean motion difference
%                  col  4:    e      ....... eccentricity
%                  col  5:    sqrt(A)  ..... where A is semimajor axis
%                  col  6:    OMEGA  ....... LoAN at weekly epoch
%                  col  7:    i0     ....... inclination at reference time
%                  col  8:    omega  ....... argument of perigee
%                  col  9:    OMEGA_dot  ... rate of right ascension 
%                  col 10:    i_dot  ....... rate of inclination angle
%                  col 11:    Cuc    ....... cosine term, arg. of latitude
%                  col 12:    Cus    ....... sine term, arg. of latitude
%                  col 13:    Crc    ....... cosine term, radius
%                  col 14:    Crs    ....... sine term, radius
%                  col 15:    Cic    ....... cosine term, inclination
%                  col 16:    Cis    ....... sine term, inclination
%                  col 17:    toe    ....... time of ephemeris
%                  col 18:    IODE   ....... Issue of Data Ephemeris
%                  col 19:    GPS_wk ....... GPS week

    fid = fopen(filename);

    if fid == -1
        errordlg(['The file ''' filename ''' does not exist.']);
        return;
    end


    % skip through header
    end_of_header = 0;
    while end_of_header == 0
        current_line = fgetl(fid);
        if strfind(current_line,'END OF HEADER')
            end_of_header=1;
        end
    end

    j = 0;
    while feof(fid) ~= 1
        j = j+1;
        current_line = fgetl(fid);
        
        % parse epoch line (ignores SV clock bias, drift, and drift rate)
        line_data = cell2mat(textscan(current_line, '%2.0f %3.0f %3.0f %3.0f %3.0f %3.0f %5.1f %19.12f %19.12f %19.12f'));
        Y = line_data(2) + 2000;
        M = line_data(3);
        D = line_data(4);
        H = line_data(5);
        min = line_data(6);
        sec = line_data(7);
        [ gps_week, tow ] = gpsweek(Y, M, D, H, min, sec);
        PRN = line_data(1);
        af0 = line_data(8);
        af1 = line_data(9);
        af2 = line_data(10);
        % Broadcast orbit line 1
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        IODE = line_data(1);
        Crs = line_data(2);
        delta_n = line_data(3);
        M0 = line_data(4);

        % Broadcast orbit line 2
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        Cuc = line_data(1);
        e = line_data(2);
        Cus = line_data(3);
        sqrtA = line_data(4);

        % Broadcast orbit line 3
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        toe = line_data(1);
        Cic = line_data(2);
        OMEGA = line_data (3);
        Cis = line_data(4);

        % Broadcast orbit line 4
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        i0 = line_data(1);
        Crc = line_data(2);
        omega = line_data (3);
        OMEGA_dot = line_data(4);

        % Broadcast orbit line 5
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        i_dot = line_data(1);
        L2_codes = line_data(2);
        GPS_wk = line_data (3);
        L2_dataflag = line_data(4);
        
        % Broadcast orbit line 6
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        SV_acc = line_data(1);
        SV_health = line_data(2);
        TGD = line_data (3);
        IODC = line_data(4);

        % Broadcast orbit line 7
        current_line = fgetl(fid);
        line_data = cell2mat(textscan(current_line, '%22.12f %19.12f %19.12f %19.12f'));
        msg_trans_t = line_data(1);
        fit_int = line_data(2);

        ephemeris(j,:) = [PRN, M0, delta_n, e, sqrtA, OMEGA, i0, omega, OMEGA_dot, i_dot, Cuc, Cus, Crc, Crs, Cic, Cis, toe, IODE, GPS_wk];
    end

    fclose(fid);

end

function ecef = gpseph2ecef(time, ephemeris)
% Converts a GPS ephemeris to the WGS84 Earth-Centered, Earth-Fixed
% (ECEF) Cartesian coordinates at the specified time.
%
    M0 = ephemeris(2);
    delta_n = ephemeris(3);
    e = ephemeris(4);
    sqrtA = ephemeris(5);
    OMEGA = ephemeris(6);
    i0 = ephemeris(7);
    omega = ephemeris(8);
    OMEGA_dot = ephemeris(9);
    i_dot = ephemeris(10);
    Cuc = ephemeris(11);
    Cus = ephemeris(12);
    Crc = ephemeris(13);
    Crs = ephemeris(14);
    Cic = ephemeris(15);
    Cis = ephemeris(16);
    toe = ephemeris(17);
    
    % WGS84 Parameters
    MU = 3.986005 * 10^14;  % m^3/s^2
    EARTH_ROT_RATE = 7.2921151467 * 10^-5;  % rad/s


    % Compute the semimajor axis.
    A = sqrtA^2;
    
    % Compute the corrected mean motion (rad/s)
    n = sqrt(MU/A^3) + delta_n;
    
    % Compute the time that has passed since the time of ephemeris.
    delta_time = time - toe;

    % Compue the mean anomaly
    M = M0 + n*delta_time;

    % Compute the eccentric anomaly using Kepler's equation
    E = M;
    ratio = 1;
    
    while abs(ratio) > 10^-8
        E_error = E - e*sin(E) - M;
        E_deriv = 1 - e*cos(E);
        ratio = E_error/E_deriv;
        E = E - ratio;
    end
    E = E + ratio;

    % Compute the true anomaly.
    nu = atan2(sqrt(1 - e^2)*sin(E) / (1 - e*cos(E)),...
        (cos(E) - e) / (1 - e*cos(E)) );

    % Compute the nominal argument of latitude.
    phi = nu + omega;

    % Compute the corrected argument of latitude.
    arg_lat = phi + Cus*sin(2*phi) + Cuc*cos(2*phi);

    % Compute the corrected radius.
    radius = A*(1 - e*cos(E)) + Crs*sin(2*phi) ...
        + Crc*cos(2*phi);

    % Compute the corrected inclination.
    incline = i0 + i_dot*delta_time ...
        + Cis*sin(2*phi) + Cic*cos(2*phi);

    % Compute the corrected longitude of the ascending node.
    loan = OMEGA + (OMEGA_dot - EARTH_ROT_RATE)...
        *delta_time - EARTH_ROT_RATE*toe;

    % Compute the orbital in-plane x-y position (ECI)
    x_eci = radius * cos(arg_lat);
    y_eci = radius * sin(arg_lat);

    % Compute the ECEF WGS84 coordinates
    ecef(1, 1) = x_eci * cos(loan) - y_eci * cos(incline) * sin(loan);
    ecef(1, 2) = x_eci * sin(loan) + y_eci * cos(incline) * cos(loan);
    ecef(1, 3) = y_eci * sin(incline);

end

function ned = ecef2ned(ecef, ref_ecef, WGS)
% Converts the WGS84 Earth-Centered, Earth-Fixed (ECEF) Cartesian 
% coordinates to the North-East-Down (NED) Cartesian coordinates using
% the reference position and latitude and longitude for the NED frame.
    
   ref_lla = ecef2lla(ref_ecef, WGS);

   delta_ecef = ecef - ref_ecef;

   north = -sind(ref_lla(1)).*cosd(ref_lla(2)).*delta_ecef(1) ...
       - sind(ref_lla(1)).*sind(ref_lla(2)).*delta_ecef(2) ...
       + cosd(ref_lla(1)).*delta_ecef(3);
   
   east = -sind(ref_lla(2)).*delta_ecef(1) ...
       + cosd(ref_lla(2)).*delta_ecef(2);

   down = -cosd(ref_lla(1)).*cosd(ref_lla(2)).*delta_ecef(1) ...
       - cosd(ref_lla(1)).*sind(ref_lla(2)).*delta_ecef(2) ...
       - sind(ref_lla(1)).*delta_ecef(3);

   ned = [north, east, down];

end

[Michael Van de Graaff]

That's waaaay to much to look though (points for not giving too little info though). how about you give like 3 rows of one of your nine 901x6 arrays

so i'm expecting something like

sample_data = [timestamp1,a1,b1,c1,d1,e1; timestamp2,a2,b2,c2,d2,e2; timestamp3,a3,b3,c3,d3,e3];

and somehting similar (without a,b,c,d,e) for the unique_time vector. (make sure there's at least one match) Then we can at least zero in on what's actually going on.

Also, you could extract the Second, Minute, Hour,Day, Month, Year, from the timestamps and just do a nested loop. Have you looked at the datetime dovumentation?

Answer 1

[Gyan Vaibhav]

Hi Terry,

I understand that you are trying to find the occurrences of your unique timestamps in the mentioned array. I assume that you are trying to extract the 2nd to 4th columns for the matched timestamps, into Data_mat.

The given code can be modified as follows to achieve the expected results:

%% Find each unique time step 
time_mat = [mat1(:, 1), mat2(:, 1), mat3(:, 1), mat4(:, 1), mat5(:, 1),mat6(:, 1), mat7(:, 1), mat8(:, 1), mat9(:, 1)]; 
time_linear = time_mat(:); 
time_sort = sort(time_linear); 
time_unique = unique(time_sort); 

% Initialize Data_mat as an empty matrix 
Data_mat = []; 

% Loop through each unique time stamp 
for i = 1:length(time_unique) 
    % Initialize a temporary matrix to store rows matching the current time stamp 
    temp_mat = []; 
    
    % Loop through each data matrix 
    for j = 1:9 
        % Find rows with matching time stamp in the j-th matrix 
        matching_rows = mat{j}(mat{j}(:, 1) == time_unique(i), :); 
        
        % Append matching rows to the temporary matrix 
        temp_mat = [temp_mat; matching_rows(:, 2:end)]; 
    end 
    
    % If there are matching rows, append them to Data_mat 
    % Otherwise, append [0, 0, 0] to Data_mat 
    if ~isempty(temp_mat) 
        Data_mat = [Data_mat; temp_mat]; 
    else 
        Data_mat = [Data_mat; 0, 0, 0]; 
    end 
end 

This approach gives us Data_mat that now contains the extracted rows of data based on the matching timestamps as expected.

Hope this helps.

Thanks

Gyan