mirror of
https://github.com/KevinMidboe/linguist.git
synced 2025-10-29 17:50:22 +00:00
Added matlab code samples.
All of these code samples currently are mis-identified in my repositories. I'm donating them to the cause.
This commit is contained in:
Michele Mastropietro
12
samples/Matlab/Check_plot.m
Normal file
12
samples/Matlab/Check_plot.m
Normal file
@@ -0,0 +1,12 @@
|
|||||||
|
x_0=linspace(0,100,101);
|
||||||
|
vx_0=linspace(0,100,101);
|
||||||
|
z=zeros(101,101);
|
||||||
|
for i=1:101
|
||||||
|
for j=1:101
|
||||||
|
z(i,j)=x_0(i)*vx_0(j);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
figure
|
||||||
|
pcolor(x_0,vx_0,z)
|
||||||
|
shading flat
|
||||||
149
samples/Matlab/FTLEH.m
Normal file
149
samples/Matlab/FTLEH.m
Normal file
@@ -0,0 +1,149 @@
|
|||||||
|
tic
|
||||||
|
clear all
|
||||||
|
%% Choice of the mass parameter
|
||||||
|
mu=0.1;
|
||||||
|
|
||||||
|
%% Computation of Lagrangian Points
|
||||||
|
[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu);
|
||||||
|
|
||||||
|
%% Computation of initial total energy
|
||||||
|
E_L1=-Omega(xl1,yl1,mu);
|
||||||
|
E=E_L1+0.03715; % Offset as in figure 2.2 "LCS in the ER3BP"
|
||||||
|
|
||||||
|
%% Initial conditions range
|
||||||
|
x_0_min=-0.8;
|
||||||
|
x_0_max=-0.2;
|
||||||
|
|
||||||
|
vx_0_min=-2;
|
||||||
|
vx_0_max=2;
|
||||||
|
|
||||||
|
y_0=0;
|
||||||
|
|
||||||
|
% Elements for grid definition
|
||||||
|
n=200;
|
||||||
|
|
||||||
|
% Dimensionless integrating time
|
||||||
|
T=2;
|
||||||
|
|
||||||
|
% Grid initializing
|
||||||
|
[x_0,vx_0]=ndgrid(linspace(x_0_min,x_0_max,n),linspace(vx_0_min,vx_0_max,n));
|
||||||
|
vy_0=sqrt(2*E+2*Omega(x_0,y_0,mu)-vx_0.^2);
|
||||||
|
|
||||||
|
% Kinetic energy computation
|
||||||
|
E_cin=E+Omega(x_0,y_0,mu);
|
||||||
|
|
||||||
|
%% Transforming into Hamiltonian variables
|
||||||
|
px_0=vx_0-y_0;
|
||||||
|
py_0=vy_0+x_0;
|
||||||
|
|
||||||
|
% Inizializing
|
||||||
|
x_T=zeros(n,n);
|
||||||
|
y_T=zeros(n,n);
|
||||||
|
px_T=zeros(n,n);
|
||||||
|
py_T=zeros(n,n);
|
||||||
|
filtro=ones(n,n);
|
||||||
|
E_T=zeros(n,n);
|
||||||
|
a=zeros(n,n); % matrix of numbers of integration steps for each integration
|
||||||
|
np=0; % number of integrated points
|
||||||
|
|
||||||
|
fprintf(' con n = %i\n',n)
|
||||||
|
|
||||||
|
%% Energy tolerance setting
|
||||||
|
energy_tol=inf;
|
||||||
|
|
||||||
|
%% Computation of the Jacobian of the system
|
||||||
|
options=odeset('Jacobian',@cr3bp_jac);
|
||||||
|
|
||||||
|
%% Parallel integration of equations of motion
|
||||||
|
parfor i=1:n
|
||||||
|
for j=1:n
|
||||||
|
if E_cin(i,j)>0 && isreal(vy_0(i,j)) % Check for real velocity and positive Kinetic energy
|
||||||
|
[t,Y]=ode45(@fH,[0 T],[x_0(i,j); y_0; px_0(i,j); py_0(i,j)],options);
|
||||||
|
% Try to obtain the name of the solver for a following use
|
||||||
|
% sol=ode45(@f,[0 T],[x_0(i,j); y_0; vx_0(i,j); vy_0(i,j)],options);
|
||||||
|
% Y=sol.y';
|
||||||
|
% solver=sol.solver;
|
||||||
|
a(i,j)=length(Y);
|
||||||
|
%Saving solutions
|
||||||
|
x_T(i,j)=Y(a(i,j),1);
|
||||||
|
px_T(i,j)=Y(a(i,j),3);
|
||||||
|
y_T(i,j)=Y(a(i,j),2);
|
||||||
|
py_T(i,j)=Y(a(i,j),4);
|
||||||
|
%Computation of final total energy and difference with
|
||||||
|
%initial one
|
||||||
|
E_T(i,j)=EnergyH(x_T(i,j),y_T(i,j),px_T(i,j),py_T(i,j),mu);
|
||||||
|
delta_E=abs(E_T(i,j)-E);
|
||||||
|
if delta_E > energy_tol; %Check of total energy conservation
|
||||||
|
fprintf(' Ouch! Wrong Integration: i,j=(%i,%i)\n E_T=%.2f \n delta_E=%.2f\n\n',i,j,E_T(i,j),delta_E);
|
||||||
|
filtro(i,j)=2; %Saving position of the point
|
||||||
|
end
|
||||||
|
np=np+1;
|
||||||
|
else
|
||||||
|
filtro(i,j)=0; % 1=interesting point; 0=non-sense point; 2= bad integration point
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
t_integrazione=toc;
|
||||||
|
fprintf(' n = %i\n',n)
|
||||||
|
fprintf(' energy_tol = %.2f\n',energy_tol)
|
||||||
|
fprintf('total \t%i\n',n^2)
|
||||||
|
fprintf('nunber \t%i\n',np)
|
||||||
|
fprintf('time to integrate \t%.2f s\n',t_integr)
|
||||||
|
|
||||||
|
%% Back to Lagrangian variables
|
||||||
|
vx_T=px_T+y_T;
|
||||||
|
vy_T=py_T-x_T;
|
||||||
|
%% FTLE Computation
|
||||||
|
fprintf('adesso calcolo ftle\n')
|
||||||
|
tic
|
||||||
|
dphi=zeros(2,2);
|
||||||
|
ftle=zeros(n-2,n-2);
|
||||||
|
|
||||||
|
for i=2:n-1
|
||||||
|
for j=2:n-1
|
||||||
|
if filtro(i,j) && ... % Check for interesting point
|
||||||
|
filtro(i,j-1) && ...
|
||||||
|
filtro(i,j+1) && ...
|
||||||
|
filtro(i-1,j) && ...
|
||||||
|
filtro(i+1,j)
|
||||||
|
|
||||||
|
dphi(1,1)=(x_T(i-1,j)-x_T(i+1,j))/(x_0(i-1,j)-x_0(i+1,j));
|
||||||
|
|
||||||
|
dphi(1,2)=(x_T(i,j-1)-x_T(i,j+1))/(vx_0(i,j-1)-vx_0(i,j+1));
|
||||||
|
|
||||||
|
dphi(2,1)=(vx_T(i-1,j)-vx_T(i+1,j))/(x_0(i-1,j)-x_0(i+1,j));
|
||||||
|
|
||||||
|
dphi(2,2)=(vx_T(i,j-1)-vx_T(i,j+1))/(vx_0(i,j-1)-vx_0(i,j+1));
|
||||||
|
|
||||||
|
if filtro(i,j)==2 % Manual setting to visualize bad integrated points
|
||||||
|
ftle(i-1,j-1)=-Inf;
|
||||||
|
else
|
||||||
|
ftle(i-1,j-1)=1/(2*T)*log(max(abs(eig(dphi'*dphi))));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Plotting results
|
||||||
|
% figure
|
||||||
|
% plot(t,Y)
|
||||||
|
% figure
|
||||||
|
% plot(Y(:,1),Y(:,2))
|
||||||
|
% figure
|
||||||
|
|
||||||
|
xx=linspace(x_0_min,x_0_max,n);
|
||||||
|
vvx=linspace(vx_0_min,vx_0_max,n);
|
||||||
|
[x,vx]=ndgrid(xx(2:n-1),vvx(2:n-1));
|
||||||
|
figure
|
||||||
|
pcolor(x,vx,ftle)
|
||||||
|
shading flat
|
||||||
|
|
||||||
|
t_ftle=toc;
|
||||||
|
fprintf('tempo per integrare \t%.2f s\n',t_integrazione)
|
||||||
|
fprintf('tempo per calcolare ftle \t%.2f s\n',t_ftle)
|
||||||
|
|
||||||
|
% save(['var_' num2str(n) '_' num2str(clock(4)])
|
||||||
|
|
||||||
|
nome=['var_xvx_', 'ode00', '_n',num2str(n),'_e',num2str(energy_tol),'_H'];
|
||||||
|
save(nome)
|
||||||
178
samples/Matlab/FTLE_reg.m
Normal file
178
samples/Matlab/FTLE_reg.m
Normal file
@@ -0,0 +1,178 @@
|
|||||||
|
tic
|
||||||
|
clear all
|
||||||
|
%% Elements for grid definition
|
||||||
|
n=100;
|
||||||
|
|
||||||
|
%% Dimensionless integrating time
|
||||||
|
T=2;
|
||||||
|
|
||||||
|
%% Choice of the mass parameter
|
||||||
|
mu=0.1;
|
||||||
|
|
||||||
|
%% Computation of Lagrangian Points
|
||||||
|
[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu);
|
||||||
|
|
||||||
|
%% Computation of initial total energy
|
||||||
|
E_L1=-Omega(xl1,yl1,mu);
|
||||||
|
C_L1=-2*E_L1; % C_L1 = 3.6869532299 from Szebehely
|
||||||
|
E=E_L1+0.03715; % Offset as in figure 2.2 "LCS in the ER3BP"
|
||||||
|
|
||||||
|
%% Initial conditions range
|
||||||
|
x_0_min=-0.8;
|
||||||
|
x_0_max=-0.2;
|
||||||
|
|
||||||
|
vx_0_min=-2;
|
||||||
|
vx_0_max=2;
|
||||||
|
|
||||||
|
y_0=0;
|
||||||
|
|
||||||
|
% Grid initializing
|
||||||
|
[x_0,vx_0]=ndgrid(linspace(x_0_min,x_0_max,n),linspace(vx_0_min,vx_0_max,n));
|
||||||
|
vy_0=sqrt(2*E+2.*Omega(x_0,y_0,mu)-vx_0.^2);
|
||||||
|
% Kinetic energy computation
|
||||||
|
E_cin=E+Omega(x_0,y_0,mu);
|
||||||
|
|
||||||
|
% Inizializing
|
||||||
|
x_T=zeros(n,n);
|
||||||
|
y_T=zeros(n,n);
|
||||||
|
vx_T=zeros(n,n);
|
||||||
|
vy_T=zeros(n,n);
|
||||||
|
filtro=ones(n,n);
|
||||||
|
E_T=zeros(n,n);
|
||||||
|
delta_E=zeros(n,n);
|
||||||
|
a=zeros(n,n); % matrix of numbers of integration steps for each integration
|
||||||
|
np=0; % number of integrated points
|
||||||
|
|
||||||
|
fprintf('integro con n = %i\n',n)
|
||||||
|
|
||||||
|
%% Energy tolerance setting
|
||||||
|
energy_tol=0.1;
|
||||||
|
|
||||||
|
%% Setting the options for the integrator
|
||||||
|
RelTol=1e-12;AbsTol=1e-12; % From Short
|
||||||
|
% RelTol=1e-13;AbsTol=1e-22; % From JD James Mireles
|
||||||
|
% RelTol=3e-14;AbsTol=1e-16; % HIGH accuracy from Ross
|
||||||
|
options=odeset('AbsTol',AbsTol,'RelTol',RelTol);
|
||||||
|
%% Parallel integration of equations of motion
|
||||||
|
h=waitbar(0,'','Name','Integration in progress, please wait!');
|
||||||
|
S=zeros(n,n);
|
||||||
|
r1=zeros(n,n);
|
||||||
|
r2=zeros(n,n);
|
||||||
|
g=zeros(n,n);
|
||||||
|
for i=1:n
|
||||||
|
waitbar(i/n,h,sprintf('Computing i=%i',i));
|
||||||
|
parfor j=1:n
|
||||||
|
r1(i,j)=sqrt((x_0(i,j)+mu).^2+y_0.^2);
|
||||||
|
r2(i,j)=sqrt((x_0(i,j)-1+mu).^2+y_0.^2);
|
||||||
|
g(i,j)=((1-mu)./(r1(i,j).^3)+mu./(r2(i,j).^3));
|
||||||
|
if E_cin(i,j)>0 && isreal(vy_0(i,j)) % Check for real velocity and positive Kinetic energy
|
||||||
|
S(i,j)=g(i,j)*T;
|
||||||
|
[s,Y]=ode45(@f_reg,[0 S(i,j)],[x_0(i,j); y_0; vx_0(i,j); vy_0(i,j)],options,mu);
|
||||||
|
a(i,j)=length(Y);
|
||||||
|
% if s(a(i,j)) < 2
|
||||||
|
% filtro(i,j)=3;
|
||||||
|
% end
|
||||||
|
% Saving solutions
|
||||||
|
x_T(i,j)=Y(a(i,j),1);
|
||||||
|
vx_T(i,j)=Y(a(i,j),3);
|
||||||
|
y_T(i,j)=Y(a(i,j),2);
|
||||||
|
vy_T(i,j)=Y(a(i,j),4);
|
||||||
|
|
||||||
|
% Computation of final total energy and difference with
|
||||||
|
% initial one
|
||||||
|
E_T(i,j)=Energy(x_T(i,j),y_T(i,j),vx_T(i,j),vy_T(i,j),mu);
|
||||||
|
delta_E(i,j)=abs(E_T(i,j)-E);
|
||||||
|
if delta_E(i,j) > energy_tol; % Check of total energy conservation
|
||||||
|
fprintf(' Ouch! Wrong Integration: i,j=(%i,%i)\n E_T=%.2f \n delta_E=%f\n\n',i,j,E_T(i,j),delta_E(i,j));
|
||||||
|
filtro(i,j)=2; % Saving position of the point
|
||||||
|
end
|
||||||
|
np=np+1;
|
||||||
|
else
|
||||||
|
filtro(i,j)=0; % 1 = interesting point; 0 = non-sense point; 2 = bad integration point
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
close(h);
|
||||||
|
t_integrazione=toc;
|
||||||
|
%%
|
||||||
|
filtro_1=filtro;
|
||||||
|
for i=2:n-1
|
||||||
|
for j=2:n-1
|
||||||
|
if filtro(i,j)==2 || filtro (i,j)==3
|
||||||
|
filtro_1(i,j)=2;
|
||||||
|
filtro_1(i+1,j)=2;
|
||||||
|
filtro_1(i-1,j)=2;
|
||||||
|
filtro_1(i,j+1)=2;
|
||||||
|
filtro_1(i,j-1)=2;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
fprintf('integato con n = %i\n',n)
|
||||||
|
fprintf('integato con energy_tol = %f\n',energy_tol)
|
||||||
|
fprintf('numero punti totali \t%i\n',n^2)
|
||||||
|
fprintf('numero punti integrati \t%i\n',np)
|
||||||
|
fprintf('tempo per integrare \t%.2f s\n',t_integrazione)
|
||||||
|
|
||||||
|
%% FTLE Computation
|
||||||
|
fprintf('adesso calcolo ftle\n')
|
||||||
|
tic
|
||||||
|
dphi=zeros(2,2);
|
||||||
|
ftle=zeros(n-2,n-2);
|
||||||
|
ftle_norm=zeros(n-2,n-2);
|
||||||
|
|
||||||
|
ds_x=(x_0_max-x_0_min)/(n-1);
|
||||||
|
ds_vx=(vx_0_max-vx_0_min)/(n-1);
|
||||||
|
|
||||||
|
for i=2:n-1
|
||||||
|
for j=2:n-1
|
||||||
|
if filtro_1(i,j) && ... % Check for interesting point
|
||||||
|
filtro_1(i,j-1) && ...
|
||||||
|
filtro_1(i,j+1) && ...
|
||||||
|
filtro_1(i-1,j) && ...
|
||||||
|
filtro_1(i+1,j)
|
||||||
|
% La direzione dello spostamento la decide il denominatore
|
||||||
|
|
||||||
|
% TODO spiegarsi teoricamente come mai la matrice pu<EFBFBD>
|
||||||
|
% essere ridotta a 2x2
|
||||||
|
dphi(1,1)=(x_T(i+1,j)-x_T(i-1,j))/(2*ds_x); %(x_0(i-1,j)-x_0(i+1,j));
|
||||||
|
|
||||||
|
dphi(1,2)=(x_T(i,j+1)-x_T(i,j-1))/(2*ds_vx); %(vx_0(i,j-1)-vx_0(i,j+1));
|
||||||
|
|
||||||
|
dphi(2,1)=(vx_T(i+1,j)-vx_T(i-1,j))/(2*ds_x); %(x_0(i-1,j)-x_0(i+1,j));
|
||||||
|
|
||||||
|
dphi(2,2)=(vx_T(i,j+1)-vx_T(i,j-1))/(2*ds_vx); %(vx_0(i,j-1)-vx_0(i,j+1));
|
||||||
|
|
||||||
|
if filtro_1(i,j)==2 % Manual setting to visualize bad integrated points
|
||||||
|
ftle(i-1,j-1)=0;
|
||||||
|
else
|
||||||
|
ftle(i-1,j-1)=(1/abs(T))*log(max(sqrt(abs(eig(dphi*dphi')))));
|
||||||
|
ftle_norm(i-1,j-1)=(1/abs(T))*log(norm(dphi));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Plotting results
|
||||||
|
% figure
|
||||||
|
% plot(t,Y)
|
||||||
|
% figure
|
||||||
|
% plot(Y(:,1),Y(:,2))
|
||||||
|
% figure
|
||||||
|
|
||||||
|
xx=linspace(x_0_min,x_0_max,n);
|
||||||
|
vvx=linspace(vx_0_min,vx_0_max,n);
|
||||||
|
[x,vx]=ndgrid(xx(2:n-1),vvx(2:n-1));
|
||||||
|
figure
|
||||||
|
pcolor(x,vx,ftle)
|
||||||
|
shading flat
|
||||||
|
|
||||||
|
t_ftle=toc;
|
||||||
|
fprintf('tempo per integrare \t%.2f s\n',t_integrazione)
|
||||||
|
fprintf('tempo per calcolare ftle \t%.2f s\n',t_ftle)
|
||||||
|
|
||||||
|
% ora=fstringf %TODO
|
||||||
|
% save(['var_' num2str(n) '_' num2str(clock(4)])
|
||||||
|
|
||||||
|
nome=['var_xvx_', 'ode00', '_n',num2str(n)];
|
||||||
|
save(nome)
|
||||||
40
samples/Matlab/Integrate1.m
Normal file
40
samples/Matlab/Integrate1.m
Normal file
@@ -0,0 +1,40 @@
|
|||||||
|
function [ x_T, y_T, vx_T, e_T, filter ] = Integrate_FILE( x_0, y_0, vx_0, e_0, T, N, mu, options)
|
||||||
|
%Integrate
|
||||||
|
% This function performs Runge-Kutta-Fehlberg integration for given
|
||||||
|
% initial conditions to compute FILE
|
||||||
|
nx=length(x_0);
|
||||||
|
ny=length(y_0);
|
||||||
|
nvx=length(vx_0);
|
||||||
|
ne=length(e_0);
|
||||||
|
vy_0=zeros(nx,ny,nvx,ne);
|
||||||
|
x_T=zeros(nx,ny,nvx,ne);
|
||||||
|
y_T=zeros(nx,ny,nvx,ne);
|
||||||
|
vx_T=zeros(nx,ny,nvx,ne);
|
||||||
|
vy_T=zeros(nx,ny,nvx,ne);
|
||||||
|
e_T=zeros(nx,ny,nvx,ne);
|
||||||
|
%% Look for phisically meaningful points
|
||||||
|
filter=zeros(nx,ny,nvx,ne); %0=meaningless point 1=meaningful point
|
||||||
|
|
||||||
|
%% Integrate only meaningful points
|
||||||
|
h=waitbar(0,'','Name','Integration in progress, please wait!');
|
||||||
|
for i=1:nx
|
||||||
|
waitbar(i/nx,h,sprintf('Computing i=%i',i));
|
||||||
|
for j=1:ny
|
||||||
|
parfor k=1:nvx
|
||||||
|
for l=1:ne
|
||||||
|
vy_0(i,j,k,l)=sqrt(2*Potential(x_0(i),y_0(j),mu)+2*e_0(l)-vx_0(k)^2);
|
||||||
|
if isreal(vy_0(i,j,k,l))
|
||||||
|
filter(i,j,k,l)=1;
|
||||||
|
ci=[x_0(i), y_0(j), vx_0(k), vy_0(i,j,k,l)];
|
||||||
|
[t,Y,te,ye,ie]=ode45(@f,[0 T], ci, options, mu);
|
||||||
|
x_T(i,j,k,l)=ye(N+1,1);
|
||||||
|
y_T(i,j,k,l)=ye(N+1,2);
|
||||||
|
vx_T(i,j,k,l)=ye(N+1,3);
|
||||||
|
vy_T(i,j,k,l)=ye(N+1,4);
|
||||||
|
e_T(i,j,k,l)=0.5*(vx_T(i,j,k,l)^2+vy_T(i,j,k,l)^2)-Potential(x_T(i,j,k,l),y_T(i,j,k,l),mu);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
close(h);
|
||||||
60
samples/Matlab/Integrate2.m
Normal file
60
samples/Matlab/Integrate2.m
Normal file
@@ -0,0 +1,60 @@
|
|||||||
|
function [ x_T, y_T, vx_T, e_T, filter, delta_e ] = Integrate_FTLE_Gawlick_ell( x_0, y_0, vx_0, e_0, T, mu, ecc, nu, options)
|
||||||
|
%Integrate
|
||||||
|
% This function performs Runge-Kutta-Fehlberg integration for given
|
||||||
|
% initial conditions to compute FTLE to obtain the image in the Gawlick's
|
||||||
|
% article "Lagrangian Coherent Structures in the Elliptic Restricted
|
||||||
|
% Three-Body Problem".
|
||||||
|
nx=length(x_0);
|
||||||
|
ny=length(y_0);
|
||||||
|
nvx=length(vx_0);
|
||||||
|
ne=length(e_0);
|
||||||
|
vy_0=zeros(nx,ny,nvx,ne);
|
||||||
|
x_T=zeros(nx,ny,nvx,ne);
|
||||||
|
y_T=zeros(nx,ny,nvx,ne);
|
||||||
|
vx_T=zeros(nx,ny,nvx,ne);
|
||||||
|
vy_T=zeros(nx,ny,nvx,ne);
|
||||||
|
e_T=zeros(nx,ny,nvx,ne);
|
||||||
|
delta_e=zeros(nx,ny,nvx,ne);
|
||||||
|
%% Look for phisically meaningful points
|
||||||
|
filter=zeros(nx,ny,nvx,ne); %0=meaningless point 1=meaningful point
|
||||||
|
useful=ones(nx,ny,nvx,ne);
|
||||||
|
%% Integrate only useful points
|
||||||
|
useful(:,1,:,1)=0;
|
||||||
|
useful(:,1,:,3)=0;
|
||||||
|
useful(:,3,:,1)=0;
|
||||||
|
useful(:,3,:,3)=0;
|
||||||
|
|
||||||
|
%% Integrate only meaningful points
|
||||||
|
h=waitbar(0,'','Name','Integration in progress, please wait!');
|
||||||
|
for i=1:nx
|
||||||
|
waitbar(i/nx,h,sprintf('Computing i=%i',i));
|
||||||
|
for j=1:ny
|
||||||
|
parfor k=1:nvx
|
||||||
|
for l=1:ne
|
||||||
|
if useful(i,j,k,l)
|
||||||
|
vy_0(i,j,k,l)=-sqrt(2*(Omega(x_0(i),y_0(j),mu)/(1+ecc*cos(nu)))+2*e_0(l)-vx_0(k)^2);
|
||||||
|
if isreal(vy_0(i,j,k,l))
|
||||||
|
filter(i,j,k,l)=1;
|
||||||
|
|
||||||
|
ci=[x_0(i), y_0(j), vx_0(k), vy_0(i,j,k,l)];
|
||||||
|
[t,Y]=ode45(@f_ell,[0 T], ci, options, mu, ecc);
|
||||||
|
|
||||||
|
if abs(t(end)) < abs(T) % Consider also negative time
|
||||||
|
filter(i,j,k,l)=3
|
||||||
|
end
|
||||||
|
|
||||||
|
x_T(i,j,k,l)=Y(end,1);
|
||||||
|
y_T(i,j,k,l)=Y(end,2);
|
||||||
|
vx_T(i,j,k,l)=Y(end,3);
|
||||||
|
vy_T(i,j,k,l)=Y(end,4);
|
||||||
|
e_T(i,j,k,l)=0.5*(vx_T(i,j,k,l)^2+vy_T(i,j,k,l)^2)-Omega(x_T(i,j,k,l),y_T(i,j,k,l),mu);
|
||||||
|
|
||||||
|
% Compute the goodness of the integration
|
||||||
|
delta_e(i,j,k,l)=abs(e_T(i,j,k,l)-e_0(l));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
close(h);
|
||||||
28
samples/Matlab/Lagr.m
Normal file
28
samples/Matlab/Lagr.m
Normal file
@@ -0,0 +1,28 @@
|
|||||||
|
function [xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu)
|
||||||
|
% [xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu)
|
||||||
|
% Lagr This function computes the coordinates of the Lagrangian points,
|
||||||
|
% given the mass parameter
|
||||||
|
yl1=0;
|
||||||
|
yl2=0;
|
||||||
|
yl3=0;
|
||||||
|
yl4=sqrt(3)/2;
|
||||||
|
yl5=-sqrt(3)/2;
|
||||||
|
c1=roots([1 mu-3 3-2*mu -mu 2*mu -mu]);
|
||||||
|
c2=roots([1 3-mu 3-2*mu -mu -2*mu -mu]);
|
||||||
|
c3=roots([1 2+mu 1+2*mu mu-1 2*mu-2 mu-1]);
|
||||||
|
xl1=0;
|
||||||
|
xl2=0;
|
||||||
|
for i=1:5
|
||||||
|
if isreal(c1(i))
|
||||||
|
xl1=1-mu-c1(i);
|
||||||
|
end
|
||||||
|
if isreal(c2(i))
|
||||||
|
xl2=1-mu+c2(i);
|
||||||
|
end
|
||||||
|
if isreal(c3(i))
|
||||||
|
xl3=-mu-c3(i);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
xl4=0.5-mu;
|
||||||
|
xl5=xl4;
|
||||||
|
end
|
||||||
16
samples/Matlab/Lagrangian_points.m
Normal file
16
samples/Matlab/Lagrangian_points.m
Normal file
@@ -0,0 +1,16 @@
|
|||||||
|
% Plot dei Lagrangian points
|
||||||
|
n=5;
|
||||||
|
mu=linspace(0,0.5,n);
|
||||||
|
for i=1:n
|
||||||
|
[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu(i));
|
||||||
|
figure (1)
|
||||||
|
hold all
|
||||||
|
plot(xl1, yl1, 's')
|
||||||
|
plot(xl2, yl2, 's')
|
||||||
|
plot(xl3, yl3, 's')
|
||||||
|
plot(xl4, yl4, 's')
|
||||||
|
plot(xl5, yl5, 's')
|
||||||
|
plot(-mu,0,'o')
|
||||||
|
plot(1-mu,0, 'o')
|
||||||
|
plot([-mu(i) xl4],[0 yl4])
|
||||||
|
end
|
||||||
18
samples/Matlab/Poincare.m
Normal file
18
samples/Matlab/Poincare.m
Normal file
@@ -0,0 +1,18 @@
|
|||||||
|
clear
|
||||||
|
%% Initial Conditions
|
||||||
|
mu=0.012277471;
|
||||||
|
T=10;
|
||||||
|
N=5;
|
||||||
|
C=3.17;
|
||||||
|
x_0=0.30910452642073;
|
||||||
|
y_0=0.07738174525518;
|
||||||
|
vx_0=-0.72560796964234;
|
||||||
|
vy_0=sqrt(-C-vx_0^2+2*Potential(x_0,y_0,mu));
|
||||||
|
k=0;
|
||||||
|
%% Integration
|
||||||
|
options=odeset('AbsTol',1e-22,'RelTol',1e-13,'Events',@cross_y);
|
||||||
|
[t,y,te,ye,ie]=ode113(@f,[0 T],[x_0; y_0; vx_0; vy_0],options,mu);
|
||||||
|
|
||||||
|
figure
|
||||||
|
%plot(ye(:,1),ye(:,3),'rs')
|
||||||
|
plot(ye(:,1),0,'rs')
|
||||||
24
samples/Matlab/RK4.m
Normal file
24
samples/Matlab/RK4.m
Normal file
@@ -0,0 +1,24 @@
|
|||||||
|
function x = RK4( fun, tspan, ci, mu )
|
||||||
|
%RK4 4th-order Runge Kutta integrator
|
||||||
|
% Detailed explanation goes here
|
||||||
|
h=1e-5;
|
||||||
|
t=tspan(1);
|
||||||
|
T=tspan(length(tspan));
|
||||||
|
dim=length(ci);
|
||||||
|
%x=zeros(l,dim);
|
||||||
|
x(:,1)=ci;
|
||||||
|
i=1;
|
||||||
|
while t<T
|
||||||
|
k1=fun(t,x(:,i),mu);
|
||||||
|
k2=fun(t+h/2,x(:,i)+k1*h/2,mu);
|
||||||
|
k3=fun(t+h/2,x(:,i)+k2*h/2,mu);
|
||||||
|
k4=fun(t+h,x(:,i)+h*k3,mu);
|
||||||
|
x(:,i+1)=x(:,i)+(h/6*(k1+2*k2+2*k3+k4));
|
||||||
|
t=t+h;
|
||||||
|
i=i+1;
|
||||||
|
end
|
||||||
|
x=x';
|
||||||
|
% function events(x)
|
||||||
|
% dist=
|
||||||
|
% return
|
||||||
|
end
|
||||||
13
samples/Matlab/distance.m
Normal file
13
samples/Matlab/distance.m
Normal file
@@ -0,0 +1,13 @@
|
|||||||
|
function [ value,isterminal,direction ] = distance( t,y,mu )
|
||||||
|
% DISTANCE compute the distance from the attactors
|
||||||
|
% [ value,terminal,direction ] = distance( t,y )
|
||||||
|
|
||||||
|
d=1e-2; % FIXME
|
||||||
|
|
||||||
|
% TODO mettere if se tolleranza D-d<tol -> value=0
|
||||||
|
D=sqrt((y(1)+mu).^2+y(2).^2); % distance from the largest primary
|
||||||
|
|
||||||
|
value=d-D;
|
||||||
|
isterminal=1;
|
||||||
|
direction=0;
|
||||||
|
end
|
||||||
49
samples/Matlab/double_gyre.m
Normal file
49
samples/Matlab/double_gyre.m
Normal file
@@ -0,0 +1,49 @@
|
|||||||
|
clear all
|
||||||
|
tic
|
||||||
|
% initialize integration time T, f(x,t), discretization size n ----------------
|
||||||
|
T = 8;
|
||||||
|
x_min=0;
|
||||||
|
x_max=2;
|
||||||
|
y_min=0;
|
||||||
|
y_max=1;
|
||||||
|
n=50; % how many points per one measure unit (both in x and in y)
|
||||||
|
ds=1/(n-1);
|
||||||
|
x_res=(x_max-x_min)*n;
|
||||||
|
y_res=(y_max-y_min)*n;
|
||||||
|
grid_x=linspace(x_min,x_max,x_res);
|
||||||
|
grid_y=linspace(y_min,y_max,y_res);
|
||||||
|
|
||||||
|
advected_x=zeros(x_res,y_res);
|
||||||
|
advected_y=zeros(x_res,y_res);
|
||||||
|
% integrate all initial points for t in [0,T] --------------------------------
|
||||||
|
parfor i = 1:x_res
|
||||||
|
for j = 1:y_res
|
||||||
|
[t,X] = ode45(@dg,[0,T],[grid_x(i),grid_y(j)]);
|
||||||
|
% store advected positions as they would appear in (x,y) coords ------
|
||||||
|
advected_x(i,j) = X(length(X(:,1)),1);
|
||||||
|
advected_y(i,j) = X(length(X(:,2)),2);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
%% Compute FTLE
|
||||||
|
sigma=zeros(x_res,y_res);
|
||||||
|
% at each point in interior of grid, store FTLE ------------------------------
|
||||||
|
for i = 2:x_res-1
|
||||||
|
for j = 2:y_res-1
|
||||||
|
% compute Jacobian phi -----------------------------------------------
|
||||||
|
phi(1,1) = (advected_x(i+1,j)-advected_x(i-1,j))/(2*ds);
|
||||||
|
phi(1,2) = (advected_x(i,j-1)-advected_x(i,j+1))/(2*ds);
|
||||||
|
phi(2,1) = (advected_y(i+1,j)-advected_y(i-1,j))/(2*ds);
|
||||||
|
phi(2,2) = (advected_y(i,j-1)-advected_y(i,j+1))/(2*ds);
|
||||||
|
% find max eigenvalue of phi'*phi ------------------------------------
|
||||||
|
lambda_max = max(abs(eig(phi'*phi)));
|
||||||
|
% store FTLE ---------------------------------------------------------
|
||||||
|
sigma(i,j) = log(lambda_max)/abs(2*T);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
toc
|
||||||
|
%% plot FTLE field ------------------------------------------------------------
|
||||||
|
figure
|
||||||
|
contourf(grid_x,grid_y,sigma');
|
||||||
|
colorbar('location','EastOutside');
|
||||||
|
axis equal
|
||||||
|
shading flat
|
||||||
41
samples/Matlab/example.m
Normal file
41
samples/Matlab/example.m
Normal file
@@ -0,0 +1,41 @@
|
|||||||
|
clear all
|
||||||
|
tic
|
||||||
|
% initialize integration time T, f(x,t), discretization size n ----------------
|
||||||
|
T = 20;
|
||||||
|
f_x_t = inline('[v(2);-sin(v(1))]','t','v');
|
||||||
|
grid_min = -3.4;
|
||||||
|
grid_max = 3.4;
|
||||||
|
grid_width = grid_max-grid_min;
|
||||||
|
n = 35;
|
||||||
|
grid_spacing = grid_min:(grid_width/(n-1)):grid_max;
|
||||||
|
advected_x=zeros(n,n);
|
||||||
|
advected_y=zeros(n,n);
|
||||||
|
% integrate all initial points for t in [0,T] --------------------------------
|
||||||
|
for i = 1:n
|
||||||
|
for j = 1:n
|
||||||
|
[t,x] = ode45(f_x_t,[0,T],[grid_spacing(i),grid_spacing(j)]);
|
||||||
|
% store advected positions as they would appear in (x,y) coords ------
|
||||||
|
advected_x(n-j+1,i) = x(length(x(:,1)),1);
|
||||||
|
advected_y(n-j+1,i) = x(length(x(:,2)),2);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
sigma=zeros(n,n);
|
||||||
|
% at each point in interior of grid, store FTLE ------------------------------
|
||||||
|
for i = 2:n-1
|
||||||
|
for j = 2:n-1
|
||||||
|
% compute Jacobian phi -----------------------------------------------
|
||||||
|
phi(1,1) = (advected_x(i,j+1)-advected_x(i,j-1))/(2*grid_width/(n-1));
|
||||||
|
phi(1,2) = (advected_x(i-1,j)-advected_x(i+1,j))/(2*grid_width/(n-1));
|
||||||
|
phi(2,1) = (advected_y(i,j+1)-advected_y(i,j-1))/(2*grid_width/(n-1));
|
||||||
|
phi(2,2) = (advected_y(i-1,j)-advected_y(i+1,j))/(2*grid_width/(n-1));
|
||||||
|
% find max eigenvalue of phi'*phi ------------------------------------
|
||||||
|
lambda_max = max(abs(eig(phi'*phi)));
|
||||||
|
% store FTLE ---------------------------------------------------------
|
||||||
|
sigma(i,j) = log(lambda_max)/abs(T);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
toc
|
||||||
|
%% plot FTLE field ------------------------------------------------------------
|
||||||
|
figure
|
||||||
|
contourf(grid_spacing,grid_spacing,sigma);
|
||||||
|
colorbar('location','EastOutside');
|
||||||
156
samples/Matlab/gpu_RKF45_FILE.m
Normal file
156
samples/Matlab/gpu_RKF45_FILE.m
Normal file
@@ -0,0 +1,156 @@
|
|||||||
|
tic
|
||||||
|
clear
|
||||||
|
%% Range definition
|
||||||
|
n=200;
|
||||||
|
|
||||||
|
mu=0.1;
|
||||||
|
[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5]=Lagr(mu);
|
||||||
|
C_L1=2*Omega(xl1,yl1,mu);
|
||||||
|
E_0=-C_L1/2+0.03715;
|
||||||
|
Y_0=0;
|
||||||
|
|
||||||
|
nx=n;
|
||||||
|
x_0_min=-0.8;
|
||||||
|
x_0_max=-0.15;
|
||||||
|
x_0=linspace(x_0_min, x_0_max, nx);
|
||||||
|
dx=(x_0_max-x_0_min)/(nx-1);
|
||||||
|
|
||||||
|
nvx=n;
|
||||||
|
vx_0_min=-2;
|
||||||
|
vx_0_max=2;
|
||||||
|
vx_0=linspace(vx_0_min, vx_0_max, nvx);
|
||||||
|
dvx=(vx_0_max-vx_0_min)/(nvx-1);
|
||||||
|
|
||||||
|
ny=3;
|
||||||
|
dy=(dx+dvx)/2;
|
||||||
|
y_0=[Y_0-dy Y_0 Y_0+dy];
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
ne=3;
|
||||||
|
de=dy;
|
||||||
|
e_0=[E_0-de E_0 E_0+de];
|
||||||
|
|
||||||
|
%% Definition of arrays of initial conditions
|
||||||
|
|
||||||
|
%In this approach, only useful pints are stored and integrated
|
||||||
|
|
||||||
|
m=1;
|
||||||
|
% x=zeros(1,nx*ny*nvx*ne);
|
||||||
|
% y=zeros(1,nx*ny*nvx*ne);
|
||||||
|
% vx=zeros(1,nx*ny*nvx*ne);
|
||||||
|
% e=zeros(1,nx*ny*nvx*ne);
|
||||||
|
% vy=zeros(1,nx*ny*nvx*ne);
|
||||||
|
filter=zeros(nx,3,nvx,3);
|
||||||
|
|
||||||
|
for i=1:nx
|
||||||
|
for j=1:ny
|
||||||
|
for k=1:nvx
|
||||||
|
for l=1:ne
|
||||||
|
v_y=-sqrt(2*Omega(x_0(i),y_0(j),mu)+2*e_0(l)-vx_0(k)^2);
|
||||||
|
if ~((j~=2) && (l~=2)) && isreal(v_y)
|
||||||
|
x(m)=x_0(i);
|
||||||
|
y(m)=y_0(j);
|
||||||
|
vx(m)=vx_0(k);
|
||||||
|
e(m)=e_0(l);
|
||||||
|
vy(m)=v_y;
|
||||||
|
filter(i,j,k,l)=1;
|
||||||
|
m=m+1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Selection of useful points
|
||||||
|
|
||||||
|
%% Data transfer to GPU
|
||||||
|
x_gpu=gpuArray(x);
|
||||||
|
y_gpu=gpuArray(y);
|
||||||
|
vx_gpu=gpuArray(vx);
|
||||||
|
vy_gpu=gpuArray(vy);
|
||||||
|
|
||||||
|
%% Integration on GPU
|
||||||
|
N=1;
|
||||||
|
t0=0;
|
||||||
|
|
||||||
|
[x_f,y_f,vx_f,vy_f]=arrayfun(@RKF45_FILE_gpu,t0,N,x_gpu,y_gpu,vx_gpu,vy_gpu,mu);
|
||||||
|
|
||||||
|
%% Data back to CPU and GPU memory cleaning
|
||||||
|
clear x_gpu y_gpu vx_gpu vy_gpu
|
||||||
|
x_T=gather(x_f);
|
||||||
|
clear x_f
|
||||||
|
y_T=gather(y_f);
|
||||||
|
clear y_f
|
||||||
|
vx_T=gather(vx_f);
|
||||||
|
clear vx_f
|
||||||
|
vy_T=gather(vy_f);
|
||||||
|
clear vy_f
|
||||||
|
|
||||||
|
%% Construction of matrix for FTLE computation
|
||||||
|
|
||||||
|
X_T=zeros(nx,ny,nvx,ne);
|
||||||
|
Y_T=zeros(nx,ny,nvx,ne);
|
||||||
|
VX_T=zeros(nx,ny,nvx,ne);
|
||||||
|
VY_T=zeros(nx,ny,nvx,ne);
|
||||||
|
E_T=zeros(nx,ny,nvx,ne);
|
||||||
|
m=1;
|
||||||
|
for i=1:nx
|
||||||
|
for j=1:ny
|
||||||
|
for k=1:nvx
|
||||||
|
for l=1:ne
|
||||||
|
if filter(i,j,k,l)==1
|
||||||
|
X_T(i,j,k,l)=x_T(m);
|
||||||
|
Y_T(i,j,k,l)=y_T(m);
|
||||||
|
VX_T(i,j,k,l)=vx_T(m);
|
||||||
|
VY_T(i,j,k,l)=vy_T(m);
|
||||||
|
E_T(i,j,k,l)=0.5*(VX_T(i,j,k,l)^2+VY_T(i,j,k,l)^2)-Omega(X_T(i,j,k,l),Y_T(i,j,k,l),mu);
|
||||||
|
m=m+1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Compute filter for FTLE
|
||||||
|
filter_ftle=filter;
|
||||||
|
for i=2:(nx-1)
|
||||||
|
for j=2:(ny-1)
|
||||||
|
for k=2:(nvx-1)
|
||||||
|
for l=2:(ne-1)
|
||||||
|
if filter(i,j,k,l)==0 || filter (i,j,k,l)==3
|
||||||
|
filter_ftle(i,j,k,l)=0;
|
||||||
|
|
||||||
|
filter_ftle(i+1,j,k,l)=0;
|
||||||
|
filter_ftle(i-1,j,k,l)=0;
|
||||||
|
|
||||||
|
filter_ftle(i,j+1,k,l)=0;
|
||||||
|
filter_ftle(i,j-1,k,l)=0;
|
||||||
|
|
||||||
|
filter_ftle(i,j,k+1,l)=0;
|
||||||
|
filter_ftle(i,j,k-1,l)=0;
|
||||||
|
|
||||||
|
filter_ftle(i,j,k,l+1)=0;
|
||||||
|
filter_ftle(i,j,k,l-1)=0;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
%% FTLE computation
|
||||||
|
|
||||||
|
[ftle, dphi]=Compute_FILE_gpu( X_T, Y_T, VX_T, E_T, dx, dy, dvx, de, N, filter_ftle);
|
||||||
|
|
||||||
|
%% Plot results
|
||||||
|
figure
|
||||||
|
FTLE=squeeze(ftle(:,2,:,2));
|
||||||
|
FTLE(1,:)=[];
|
||||||
|
% FTLE(2,:)=[];
|
||||||
|
FTLE(:,1)=[];
|
||||||
|
% FTLE(:,2)=[];
|
||||||
|
x_0(1)=[];
|
||||||
|
vx_0(1)=[];
|
||||||
|
pcolor(x_0, vx_0, FTLE')
|
||||||
|
shading flat
|
||||||
|
toc
|
||||||
19
samples/Matlab/test_rk_par.m
Normal file
19
samples/Matlab/test_rk_par.m
Normal file
@@ -0,0 +1,19 @@
|
|||||||
|
clear
|
||||||
|
mu=0.1;
|
||||||
|
x_0=linspace(-0.8, -0.15, 2)
|
||||||
|
y_0=zeros(1,2)
|
||||||
|
vx_0=linspace(-2, 2, 2)
|
||||||
|
vy_0=zeros(1,2)
|
||||||
|
ci=[1-mu-0.05 0 0.005 0.5290]
|
||||||
|
t0=[0;0]
|
||||||
|
T=[2;2]
|
||||||
|
tspan=2
|
||||||
|
arg1={@f;@f}
|
||||||
|
%tspan={[0 2],[0 2]};
|
||||||
|
arg=[mu;mu]
|
||||||
|
[X]=arrayfun(RK4_par,t0,T,x_0',y_0',vx_0',vy_0',arg)
|
||||||
|
% [X]=arrayfun(@f,[0;1],[ci;ci],[mu;mu]);
|
||||||
|
%Y=RK4(@f,tspan,ci,mu);
|
||||||
|
% figure
|
||||||
|
% plot(Y(:,1),Y(:,2))
|
||||||
|
% Y(end,1)
|
||||||
Reference in New Issue
Block a user