The output size of the lambert_equation function that follows should be 1x1x3 but it is 1x3x3. Does anyone identify the problem? Thanks, in advance!

% The following syntax is used inside a double-for loop in the main script:

[v1_2(i,j,:), v2_prime_2(i,j,:)] = lambert_equation(r1(i,:), r2_prime(i,j,:), delta_t(i,j), string);

% lambert_euation function:

function [V1, V2] = lambert_equation(R1, R2, t, str)
%   This function solves Lambert's problem.
%   INPUTS: 
%       R1 = position vector at position 1 (km)
%       R2 = position vector at position 2 (km)
%       t  = time of flight  (s) 
%   OUTPUTS: 
%       V1 = velocity vector at position 1 (km/s)
%       V2 = velocity vector at position 2 (km/s)
%
%   VARIABLES DESCRIPTION:
%       r1, r2     - magnitudes of R1 and R2 
%       c12        - cross-product of R1 and R2 
%       theta      - change in true anomaly between position 1 and 2 
%       z          - alpha*x^2, where alpha is the reciprocal of the
%                    semimajor axis and x is the universal anomaly
%       y(z)       - a function of z
%       F(z,t)     - a function of the variable z and constant t
%       dFdz(z)    - the derivative of F(z,t)
%       ratio      - F/dFdz
%       eps        - tolerence on precision of convergence
%       C(z), S(z) - Stumpff functions

%% Gravitational parameter:    
global mu

%% Magnitudes of R1 and R2:

r1 = norm(R1);
r2 = norm(R2,'fro');

%% Theta:

c12 = cross(R1, squeeze(R2));
theta = acos(dot(R1,squeeze(R2))/r1/r2);
%Assume prograde or retrograde orbit:
if strcmp(str,'pro')
    if c12(3) <= 0
        theta = 2*pi - theta;
    end
elseif strcmp(str,'retro')
    if c12(3) >= 0
        theta = 2*pi - theta;
    end
end

%% Solve equation for z:

A = sin(theta)*sqrt(r1*r2/(1 - cos(theta)));
%Check where F(z,t) changes sign, and use that 
%value of z as a starting point:
z = -100;
while F(z,t) < 0
    z = z + 0.1;
end
%Tolerence and number of iterations:
eps= 1.e-8;
nmax = 5000;
%Iterate until z is found under the imposed tolerance:
ratio = 1;
n = 0;
while (abs(ratio) > eps) && (n <= nmax)
    n = n + 1;
    ratio = F(z,t)/dFdz(z);
    z = z - ratio;
end
%Maximum number of iterations exceeded
if n >= nmax
    fprintf('\n\n **Number of iterations exceeds %g \n\n ',nmax)
end

%% Lagrange coefficents:
f = 1 - y(z)/r1;
g = A*sqrt(y(z)/mu);
g_dot = 1 - y(z)/r2;

%% V1 and V2:

V1 =1/g*(R2 - f*R1);
V2 =1/g*(g_dot*R2 - R1);

return

%% Subfunctions used: 
function dum = y(z)
    dum = r1 + r2 + A*(z*S(z) - 1)/sqrt(C(z));
end

function dum = F(z,t)
    dum = (y(z)/C(z))^1.5*S(z) + A*sqrt(y(z)) - sqrt(mu)*t;
end

function dum = dFdz(z)
    if z == 0
        dum = sqrt(2)/40*y(0)^1.5 + A/8*(sqrt(y(0)) + A*sqrt(1/2/y(0)));
    else
        dum = (y(z)/C(z))^1.5*(1/2/z*(C(z) - 3*S(z)/2/C(z)) ...
        + 3*S(z)^2/4/C(z)) + A/8*(3*S(z)/C(z)*sqrt(y(z)) ...
        + A*sqrt(C(z)/y(z)));
    end
end

%Stumpff functions:
function c = C(z)
    if z > 0
        c = (1 - cos(sqrt(z)))/z;
    elseif z < 0
        c = (cosh(sqrt(-z)) - 1)/(-z);
    else
        c = 1/2;
    end
end

function s = S(z)
    if z > 0
        s = (sqrt(z) - sin(sqrt(z)))/(sqrt(z))^3;
    elseif z < 0
        s = (sinh(sqrt(-z)) - sqrt(-z))/(sqrt(-z))^3;
    else
        s = 1/6;
    end
end
end

[Sam Chak]

I haven't delved into the equations in the code yet as the description of the technical issue might be insufficient. However, I believe that Lambert's Problem ultimately relates to determining the orbit based on two position vectors and the time taken to travel between them.

From a mathematical standpoint, this can be seen as an advanced curve-fitting problem where you aim to fit an orbit (defined by the six orbital elements) to the given two position vectors (two data points), although the output of your code ends with returning only two velocity vectors.

Answer 1

[SAI SRUJAN]

Hi Ioannis,

I undertand that you are facing an issue with the output dimensions returned by the function 'lambert_equation'.

In the provided function, 'R2' is being squeezed using 'squeeze(R2)' in two places. This squeezing operation changes the dimensions of 'R2' from potentially 1x3 to 3, affecting subsequent operations that expect 'R2' to maintain its original dimensions.

The real issue likely comes from how the outputs 'V1' and 'V2' are calculated:

V1 = 1/g*(R2 - f*R1);
V2 = 1/g*(g_dot*R2 - R1);

Ensure that the inputs 'R1' and 'R2' (particularly 'R2' as it's being squeezed) are of correct dimensions when they are passed into the 'lambert_equation' function. Make sure that the dimensions of 'R2' is as expected after it is being squeezed.

You can go through the following work around to manipulate the dimensions of 'V1' and 'V2' to get the desired dimension output.

% Preallocate V1 and V2
V1 = zeros(1, 1, 3);
V2 = zeros(1, 1, 3);

[V1(i,j,:), V2(i,j,:)] = lambert_equation(r1(i,:), r2_prime(i,j,:), delta_t(i,j), string);

For a comprehensive understanding of the 'squeeze' function in MATLAB, please go through the following documentation.

https://www.mathworks.com/help/matlab/ref/squeeze.html

I hope this helps!