Row versus column benchmark

rowVersusColumnBenchmark.m

function [rowTime, colTime] = rowVersusColumnBenchmark(n, m)
massRatio = 0.012;

rng(10)
states = rand(n, 6);

tic()
for k = 1:m
    pp = pseudopotentialRowStates(massRatio, states);
end
rowTime = toc();

states = transpose(states);

tic()
for k = 1:m
    pp = pseudopotentialColumnStates(massRatio, states);
end
colTime = toc();

end


function omega = pseudopotentialRowStates(massRatio, q)
%CRTBPPSEUDOPOTENTIAL compute the CRTBP pseudopotential
%
%   Author: Rolfe Power <rpower@purdue.edu>
arguments
    massRatio (:, 1) double
    q (:, 6) double
end

% Concision
mu = massRatio;
omm = 1 - mu;
q1 = q(:, 1); % X
q2 = q(:, 2); % Y
q3 = q(:, 3); % Z

yz2 = sum(q(:, 2:3).^2, 2);
d1 = sqrt((0 - massRatio - q(:, 1)).^2 + yz2);
d2 = sqrt((1 - massRatio - q(:, 1)).^2 + yz2);

omega = omm ./ d1 + mu ./ d2 + 0.5 * (q1.^2 + q2.^2);
end

function omega = pseudopotentialColumnStates(massRatio, q)
%CRTBPPSEUDOPOTENTIAL compute the CRTBP pseudopotential
%
%   Author: Rolfe Power <rpower@purdue.edu>
arguments
    massRatio (:, 1) double
    q (6, :) double
end

% Concision
mu = massRatio;
omm = 1 - mu;
q1 = q(1, :);
q2 = q(2, :);
q3 = q(3, :);

yz2 = sum(q(2:3, :).^2, 2);
d1 = sqrt((0 - massRatio - q(1, :)).^2 + yz2);
d2 = sqrt((1 - massRatio - q(1, :)).^2 + yz2);

omega = omm ./ d1 + mu ./ d2 + 0.5 * (q1.^2 + q2.^2);
end