How do I combine my function's output matrices into a single matrix?

function [] = integrateQuaternions(BR,EA)

% Parameters

A = 0; % Initial Conditions (angles halved and in rad)
B = 0;
C = -45 * (pi/180);
p = 0 * (pi/180); % Body rates in radians (0)
q = 5.5 * (pi/180); % (5.N8 deg/s)
r = 3.5 *(pi/180); % ((1+0.5*N9) deg/s)


BR = [p q r];
EA = [0 0 (-45*(pi/180))];

for i = 1:0.1:30 % for 30 second period with 0.1s timestep

% Initial Attitude

q0 = i.*(cosd(C)*cosd(B)*cosd(A) + sind(C)*sind(B)*sind(A));
q1 = i.*(cosd(C)*cosd(B)*sind(A) - sind(C)*sind(B)*cosd(A));
q2 = i.*(cosd(C)*sind(B)*cosd(A) + sind(C)*cosd(B)*sind(A));
q3 = i.*(-cosd(C)*sind(B)*sind(A) + sind(C)*cosd(B)*cosd(A));

Q = [q0 q1 q2 q3];

%Euler Interation

q0dot = i.*(-0.5*((q1*p)+(q2*q)+(q3*r)));
q1dot = i.*(0.5*((q0*p)-(q3*q)+(q2*r)));
q2dot = i.*(0.5*((q3*p)+(q0*q)-(q1*r)));
q3dot = i.*(-0.5*((q2*p)-(q1*q)-(q0*r)));

% Normalizing

Qdot = [q0dot q1dot q2dot q3dot];

mu = sqrt(q0^2 + q1^2 + q2^2 + q3^2);

qnplus1 = Q + Qdot;

Qx = (qnplus1./mu);


Qfinal = [Qx]

end
end

My function is working correctly as it gives the matrices I want, although I'd like to combine them into a single output matrix but am not sure how to. Any advice would be greatly appreciated!

Answer 1

[Dave B]

Your function loops over some values and for each one computes a row vector.

To store the row vector in a matrix, specify an integer row. I've also adjusted your function to return the value:

q = integrateQuaternions;
size(q)

ans = 1×2

291 4

function Qfinal = integrateQuaternions(BR,EA)
    % Parameters
    A = 0; % Initial Conditions (angles halved and in rad)
    B = 0;
    C = -45 * (pi/180);
    p = 0 * (pi/180); % Body rates in radians (0)
    q = 5.5 * (pi/180); % (5.N8 deg/s)
    r = 3.5 *(pi/180); % ((1+0.5*N9) deg/s)
    BR = [p q r];
    EA = [0 0 (-45*(pi/180))];
    
    QFinal = nan(length(1:.1:30),4);
    row = 1;
    for i = 1:0.1:30 % for 30 second period with 0.1s timestep
        % Initial Attitude
        q0 = i.*(cosd(C)*cosd(B)*cosd(A) + sind(C)*sind(B)*sind(A));
        q1 = i.*(cosd(C)*cosd(B)*sind(A) - sind(C)*sind(B)*cosd(A));
        q2 = i.*(cosd(C)*sind(B)*cosd(A) + sind(C)*cosd(B)*sind(A));
        q3 = i.*(-cosd(C)*sind(B)*sind(A) + sind(C)*cosd(B)*cosd(A));
        Q = [q0 q1 q2 q3];
        %Euler Interation
        q0dot = i.*(-0.5*((q1*p)+(q2*q)+(q3*r)));
        q1dot = i.*(0.5*((q0*p)-(q3*q)+(q2*r)));
        q2dot = i.*(0.5*((q3*p)+(q0*q)-(q1*r)));
        q3dot = i.*(-0.5*((q2*p)-(q1*q)-(q0*r)));
        % Normalizing
        Qdot = [q0dot q1dot q2dot q3dot];
        mu = sqrt(q0^2 + q1^2 + q2^2 + q3^2);
        qnplus1 = Q + Qdot;
        Qx = (qnplus1./mu);
        Qfinal(row,:) = Qx;
        row = row + 1;
    end
end

Note that a more common MATLAB approach looks more like

% timesteps = 1:.1:30;
% for i = 1:length(timesteps)
    % t = timesteps(i);
    % q0 = t .* ...
    % ...
    % Qfinal(i,:) = ...
% end

[jrz]

Thank you so much for your reply!

Answer 2

[Walter Roberson]

function Qfinal = integrateQuaternions(BR,EA)
    % Parameters
    
    A = 0; % Initial Conditions (angles halved and in rad)
    B = 0;
    C = -45 * (pi/180);
    p = 0 * (pi/180); % Body rates in radians (0)
    q = 5.5 * (pi/180); % (5.N8 deg/s)
    r = 3.5 *(pi/180); % ((1+0.5*N9) deg/s)
    
    
    BR = [p q r];
    EA = [0 0 (-45*(pi/180))];
    
    ivals = 1:0.1:30; % for 30 second period with 0.1s timestep
    num_i = length(ivals);
    Qfinal = zeros(1, num_i);
    
    for iidx = 1 : num_i % for 30 second period with 0.1s timestep
        i = ivals(iidx);
        % Initial Attitude
        
        q0 = i.*(cosd(C)*cosd(B)*cosd(A) + sind(C)*sind(B)*sind(A));
        q1 = i.*(cosd(C)*cosd(B)*sind(A) - sind(C)*sind(B)*cosd(A));
        q2 = i.*(cosd(C)*sind(B)*cosd(A) + sind(C)*cosd(B)*sind(A));
        q3 = i.*(-cosd(C)*sind(B)*sind(A) + sind(C)*cosd(B)*cosd(A));
        
        Q = [q0 q1 q2 q3];
        
        %Euler Interation
        
        q0dot = i.*(-0.5*((q1*p)+(q2*q)+(q3*r)));
        q1dot = i.*(0.5*((q0*p)-(q3*q)+(q2*r)));
        q2dot = i.*(0.5*((q3*p)+(q0*q)-(q1*r)));
        q3dot = i.*(-0.5*((q2*p)-(q1*q)-(q0*r)));
        
        % Normalizing
        
        Qdot = [q0dot q1dot q2dot q3dot];
        
        mu = sqrt(q0^2 + q1^2 + q2^2 + q3^2);
        
        qnplus1 = Q + Qdot;
        
        Qx = (qnplus1./mu);
        
        
        Qfinal(iidx) = [Qx];
    
    end
end