MYTE

clear all; close all; clc;

%% Define a random velocity field
[X,Y] = meshgrid(-3*pi:pi/4:3*pi,-3*pi:pi/4:3*pi);
U = sin(Y);
V = cos(X);

%% Plot the velocity field
figure(1),
quiver(X,Y,U,V,0,'autoscale','off'); 
title('Vector Field with Particle Trajectory starting from point (x,y)', 'fontsize',20)
axis equal; grid on; hold on

%%  Initial
% 
%   for x = 1:10
%       disp(x)
%   end
% 
%   for x = 1:10
%       disp(x)
%   end
% 
% 
a = [0 0; 5 5;];
%% For loop ALL
tF_test = [15, 5];
for i = 1:2
    % Solving with ODE
    p0(1) = a(i,1);
    p0(2) = a(i,2);
    f = @(t,p) odefun_p(t,p,U,V,X,Y);
    tspan = [0,tF_test(i)];

    % assign options for solver
    options = odeset('RelTol',1e-5,'AbsTol',1e-5);
    [t,p] = ode45(f,tspan,p0, options);

    % overlay the particle track
    plot(p(:,1),p(:,2),'b-','LineWidth',4,'MarkerSize',10)
    i;
end

%% Function: PARTICLE 1
function dpdt = odefun_p(~,p,U,V,X,Y,extrap)

xq = p(1);
yq = p(2);

% interpolate to get velocity components at query points
u = interp2(X,Y,U,xq,yq,"cubic");
v = interp2(X,Y,V,xq,yq,"cubic");

% assign derivative vector
dpdt = [u;v];
end