I am answering my own question:
classdef Ball < handle
%%BAll Class
properties (GetAccess = public, SetAccess = private)
mass;
diameter;
rho;
g;
Cd;
initialPosition;
initialVelocity;
position;
time;
end
properties (GetAccess = private, SetAccess = private)
velocity;
state;
end
methods
%%constructor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function obj = Ball(ballData)
obj.initialPosition = ballData.initialPosition;
obj.initialVelocity = ballData.initialVelocity;
obj.mass = ballData.mass;
obj.diameter = ballData.diameter;
obj.rho = ballData.rho;
obj.g = ballData.g;
obj.Cd = ballData.Cd;
obj.time = 0.0;
obj.state=[obj.initialVelocity; obj.initialPosition ];
obj.position=obj.initialPosition;
end
%%update %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function update(obj, world)
%update the ball to synchronize it with the world
deltaT = world.time - obj.time;
if deltaT > 0
% Solving ODE
[~,Y] = ode45(@obj.ballDynamics,[0 deltaT],obj.state);
[i, ~] = size(Y);
obj.state = Y(i,:)';
obj.position = obj.state(4:6);
% Update ball time;
obj.time = world.time;
end
end
function x_dot = ballDynamics(obj,~,x)
A = pi*(obj.diameter^2)/4;
x_dot = [obj.g - (obj.rho*obj.Cd*A*norm(x(1:3))/(2*obj.mass)) * x(1:3);
x(1:3)];
end
end
end