% GENERATEUR DE TURBULENCE SYNTHETIQUE
% On considère des vortex isolés modèles, formés d'un coeur de rankine de 
% rayon rvortex et entourés d'un anneau négatif de vorticité sur une 
% largeur rvortex. La vitesse maximale de ces vortex est notée uvortex. 
% Ils sont placés de façon totalement aléatoire dans un 
% rectangle de longueur Lx et de largeur Ly. Un capteur ponctuel (ADV ou
% LDV) mesurant tous les dtacq (fréquence d'acquisition en 1/dtacq) est 
% ensuite déplacé à vitesse constante U le long de la direction x,
% et les mesures de vitesse (deux composantes) sont extraites de cette
% translation pour simuler le signal temporel obtenu par ce moyen.
% Cela correspond également à une mesure en un point fixe, au travers
% duquel les tourbillons passent en étant advectés à la vitesse moyenne U
% (hypothèse de turbulence gelée de Taylor).

% PARAMETRES
U=1.0;              % vitesse de l'écoulement (ou de déplacement du capteur)
dtacq=0.001;        % intervalle de temps entre deux mesures du capteur
Lx=200.0;           % longueur de la zone remplie de tourbillons
Ly=4.0;             % largeur de la zone remplie de tourbillons
Nvortex=500;        % nombre de tourbillons pour le remplissage
uvortex=1.0;        % vitesse maximale des tourbillons
rvortex=1.0;        % taille des vortex (coeur et anneau négatif de vorticité)

% CALCULS
Ntps=Lx/U/dtacq;    % nombre total de mesures
xyvortex=zeros(Nvortex,2);
for i=1:Nvortex
    xloc=rand(1)*Lx;
    yloc=rand(1)*Ly;
    xyvortex(i,1)=xloc;
    xyvortex(i,2)=yloc;
end
utab=zeros(Ntps,1);
vtab=zeros(Ntps,1);
for itps=1:Ntps
    xloc=itps*dtacq;
    yloc=Ly/2.0;
    uloc=0.0;
    vloc=0.0;
    for i=1:Nvortex
        dloc=sqrt((xloc-xyvortex(i,1))^2+(yloc-xyvortex(i,2))^2);
        xp=xloc-xyvortex(i,1);
        yp=yloc-xyvortex(i,2);
        if ((dloc<rvortex))
            uloc=uloc+(dloc/rvortex*uvortex)*yp/dloc;
            vloc=vloc+(dloc/rvortex*uvortex)*xp/dloc;
        elseif (dloc<(2.0*rvortex))
            uloc=uloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*yp/dloc;            
            vloc=vloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*xp/dloc;
        end
    end
    for i=1:Nvortex
        dloc=sqrt((xloc-xyvortex(i,1)-Lx)^2+(yloc-xyvortex(i,2))^2);
        xp=xloc-xyvortex(i,1)-Lx;
        yp=yloc-xyvortex(i,2);
        if ((dloc<rvortex))
            uloc=uloc+(dloc/rvortex*uvortex)*yp/dloc;
            vloc=vloc+(dloc/rvortex*uvortex)*xp/dloc;
        elseif (dloc<(2.0*rvortex))
            uloc=uloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*yp/dloc;            
            vloc=vloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*xp/dloc;
        end
    end
    for i=1:Nvortex
        dloc=sqrt((xloc-xyvortex(i,1)+Lx)^2+(yloc-xyvortex(i,2))^2);
        xp=xloc-xyvortex(i,1)+Lx;
        yp=yloc-xyvortex(i,2);
        if ((dloc<rvortex))
            uloc=uloc+(dloc/rvortex*uvortex)*yp/dloc;
            vloc=vloc+(dloc/rvortex*uvortex)*xp/dloc;
        elseif (dloc<(2.0*rvortex))
            uloc=uloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*yp/dloc;            
            vloc=vloc+((-1.0/3.0*dloc/rvortex+4.0/3.0*rvortex/dloc)*uvortex)*xp/dloc;
        end
    end
    utab(itps)=uloc;
    vtab(itps)=vloc;
end
tabtps=(1:Ntps)*dtacq;
figure;
plot(tabtps,utab);
hold on
plot(tabtps,vtab,'g');