Passing more than two inputs to a MATLAB event function

I'm attempting to pass three parameters to my MATLAB event function - as opposed to just the normal two (t,y) - in order to solve an orbital dynamics problem, but the method I currently have implemented does not work. I believe the following code should work, but when run, the event never records an iteration step, time, or state value (y) at which the event occured. I'm certain that the event does occur within the space of the simulation, because I can find where the event occurs by simply simulating over a very large time and using a for loop to find where the event occurs. (i.e. set paramater as found in the event function is reached) This is of course extremely inefficient, so utilizing an event function would be far better. Code segments are attached and notated below from my current script:

% Attempting to solve the ode and find where the event occurs. 'T_Spec_Energy_ZGP' is an important variable I need to pass through to the event function, but it's really throwing a wrench in getting this code to run.

IB_Term_Pass=@(t,y) IB_Term(t,y,T_Spec_Energy_ZGP);
options_IB=odeset('MaxStep',10,'Events',IB_Term_Pass);
[t_IB,y_IB,t_EC,ye_IB,i_IB_f]=ode45(@(t,y) rates(t,y,r_OE,r_OM,Craft_Thrust,m_craft),[0,Simulation_Time],y0,options_IB);
-----------extra unrelated code----------------
%% Functions Section:
function [value,isterminal,direction]=IB_Term(t,y,T_Spec_Energy_ZGP)
% Event function to terminate the Initial Burn Phase (1) when spacecraft
%   has reached the Tolerance specific energy required.

% Orbital Parameters:
u_M=4903*(10^3)^3; %m^3/s^2 - Lunar Grav Parameter
mM=73.48*10^21; %Mass of the Moon (kg)
aM=384.4*10^3; %Semimajor axis of Lunar orbit (km)
mE=5.974*10^24; %Mass of Earth (kg)
r_OE=(-((10^3*aM)*(mM/mE))/(1+(mM/mE))); %Distance from the origin to the Earth (m)
r_OM=((10^3*aM)+r_OE); %Distance from the origin to the Moon (m)
u_3=398600*(10^3)^3; %Earth Gravitation Constant (m^3)/(s^2)
%Spacecraft Parameters:
r_EC_IB(1)=(-r_OE+y(1));
r_EC_IB(2)=(y(2));
r_EC_IB_mag=sqrt((r_EC_IB(1))^2+(r_EC_IB(2))^2);
r_MC_IB(1)=(-r_OM+y(1));
r_MC_IB(2)=(y(2));
r_MC_IB_mag=sqrt((r_MC_IB(1))^2+(r_MC_IB(2))^2);
v_mag=sqrt((y(3)^2)+(y(4)^2));
Spec_Energy_Craft_IB=((v_mag^2)/2)-((u_M*(10^3)^3)/r_MC_IB_mag)-((u_3*(10^3)^3)/r_EC_IB_mag);
value=(T_Spec_Energy_ZGP-Spec_Energy_Craft_IB);
direction=0;
isterminal=1;
end

Essentially my main concern is retrieving the iteration step where the event occurs, so in the above code's nomenclature I really need the output 'i_IB_f'. There's a bunch of other code, I've left out what I thought was irrelevent to this specific problem, but I've run out of ideas on where the issue could be. Any help is hugely appreciated, thank you!

[Walter Roberson]

That looks okay to me at the moment. That should stop at any zero crossing where T_Spec_Energy_ZGP is close enough to Spec_Energy_Craft_IB

I have no idea if the equations are correct, though.

You might want to plot value inside IB_Term -- you could use a persistent variable to hold the handle of an axes. Just make sure not to make it the current axis -- for example,

plot(ax, x, y)

instead of

axes(ax)
plot(x, y)

Though you would probably want to use animatedline() instead of adding individual points.

[Spencer Davidson]

Plotted through the IB_Term value, in the figure below the red curve represents the T_Spec_Energy_ZGP which is where I want the event to occur, and the blue curve represents the craft's specific energy 'Spec_Energy_Craft_IB' over my simulation timespan.

Clearly these curves intersect at some time that the event function should recognize, but it just leaves my t_EC,y_EC, and i_EC arrays empty after simulating through the entire timespan. Any recommendations on where to go from here?

[Torsten]

Your code looks fine. So nothing can be said if we are not able to reproduce the failure of the event facility.

Are you sure the code to compute the blue curve above and the code to compute "Spec_Energy_Craft_IB" in the events function are really identical ?

Are you sure the value for "T_Spec_Energy_ZGP" handed to the events function equals the y-value of the red curve above ?

Answer 1

[Nipun]

Hi Spencer,

I understand that you are using the odeset function to specify an event function (IB_Term_Pass) for the ode45 solver. The event function is expected to return three outputs: value, isterminal, and direction. The event is considered to have occurred when value crosses zero (or a specified value) with the correct sign specified by direction. The integration is terminated if isterminal is set to 1.

Based on the attached code, I infer that in your IB_Term function, you have set isterminal=1, which means the integration should terminate when the event occurs. However, you are also setting direction=0, which means the event is triggered irrespective of the sign of value. This could potentially cause the event to be triggered at each integration step.

Consider appending your code as below to change the direction:

function [value, isterminal, direction] = IB_Term(t, y, T_Spec_Energy_ZGP)
    % ... (your code)
    
    value = T_Spec_Energy_ZGP - Spec_Energy_Craft_IB;
    
    % Set direction based on the condition you are checking
    if value > 0
        direction = 1;  % Event occurs when value becomes positive
    else
        direction = -1; % Event occurs when value becomes negative
    end
    
    isterminal = 1;
end

Hope this helps

Regards,

Nipun

[Walter Roberson]

This could potentially cause the event to be triggered at each integration step.

ode45() records the sign() of each value returned by the event function. The next time the event function is called, it records the sign() again. If the old sign and new sign are the same, nothing happens. If the old sign and the new sign are different, then it checks the direction: if direction is -1 and the new sign is positive, or if direction is +1 and the new sign is negative, or if direction is 0 regardless of the new sign, then an event is triggered.

With isterminal set and direction 0, then any zero crossing ends the integration. Which is fine, even if it is almost immediate.