clc
clear all
close all
 
L=pi;                       % dimensione della singola cavità
np=1000;                    % numero di particelle
tmax=80;                    % tempo di simulazione
dt=0.1;                     % intervallo di tempo
npt=tmax/dt + 1;            % numero di istanti temporali
DATI = zeros(npt,5,np);     % archiviazione dati
xyi=2*L*rand(np,2);                             % coordinate iniziali casuali
uvi(:,:)= [sin(xyi(:,1)).*cos(xyi(:,2)) ,...
     -cos(xyi(:,1)).*sin(xyi(:,2))] ;           % velocità iniziale

% parametri di simulazione
R = 0;
W = 0.5;
A = 2;
B=10; 
% Q=0.5;
% A = R*B/2;
% W = (1/R-3/2)*Q;

% risoluzione traiettoie per le np particelle 
for i=1:np
    
    [tempo,x_y_u_v] = ode45(@MaxeyRiley, 0:dt:tmax , [xyi(i,1) xyi(i,2) uvi(i,1) uvi(i,2)]); 
    DATI(:,:,i) = [tempo, x_y_u_v]; 
    clc
    i

end

% calcolo velocità media <V2> all'istante finale
v2_m = 0;
for i=1:np
    v2_m = v2_m + DATI(npt,5,i);
end
v2_m = v2_m/np;

% reimmisione nel dominio delle particelle che ne escono
for i=1:np
    for j=1:length(tempo)
        
        % x>2L
        if DATI(j,2,i) > 2*L
            DATI(j,2,i) = DATI(j,2,i) - 2*L * fix(DATI(j,2,i)/(2*L));
        end   
        
        % x<0
        if DATI(j,2,i) < 0
            DATI(j,2,i) = DATI(j,2,i) - 2*L * floor(DATI(j,2,i)/(2*L));
        end  
        
        % y>2L
        if DATI(j,3,i) > 2*L
            DATI(j,3,i) = DATI(j,3,i) - 2*L * fix(DATI(j,3,i)/(2*L));
        end  
            
        % y<0
        if DATI(j,3,i) < 0
            DATI(j,3,i) = DATI(j,3,i) - 2*L * floor(DATI(j,3,i)/(2*L));
        end  
            
    end

end

% salvataggio dati
lift = 0;
R = 0;
W = 0.5;
A = 2;

S = zeros(npt,np);
for i=1:5
   S(:,:)=DATI(:,i,:);
   nome = strcat( strcat( strcat( 'R', strcat( num2str(R) , strcat( '_B' , strcat( num2str(B) , strcat('_Q', strcat(num2str(Q),'_') ) ) ) ) ) , strcat( strcat('lift', num2str(lift) ) , strcat('__',num2str(i)) ) ) ,'.txt' );
   save(nome,'S','-ascii');  
end

disp(sprintf('*** FINE ***'))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% EQUAZIONE DI MAXEY & RILEY %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function df=MaxeyRiley(b,f,R,A,W)

df=zeros(4,1);
lift = 0;

df(1)=sin(f(1))*cos(f(2));
df(2)=-sin(f(2))*cos(f(1));

df(1)=f(3);
df(2)=f(4);
df(3)=-A*f(3)+...
    A*sin(f(1))*cos(f(2))+...
    1/2*R*(f(3)*cos(f(1))*cos(f(2))-...
    f(4)*sin(f(1))*sin(f(2)))+...
    R*sin(f(1))*cos(f(1))+...
    lift*0.727*sqrt(R*A)*((-f(4)-cos(f(1))*sin(f(2)))*...
    (2*sin(f(1))*sin(f(2))))/sqrt(2*abs(sin(f(1))*sin(f(2))));
df(4)=-A*f(4)-...
    A*sin(f(2))*cos(f(1))+...
    1/2*R*(f(3)*sin(f(1))*sin(f(2))-...
    f(4)*cos(f(1))*cos(f(2)))+...
    R*sin(f(2))*cos(f(2))+...
    A*W+...
    lift*0.727*sqrt(R*A)*((f(3)-cos(f(2))*sin(f(1)))*...
    (2*sin(f(1))*sin(f(2))))/sqrt(2*abs(sin(f(1))*sin(f(2))));

end
