// Copyright J. Fourcade (volucres.fr)
// Ce logiciel est régi par la licence CeCILL soumise au droit français et
// respectant les principes de diffusion des logiciels libres. Vous pouvez
// utiliser, modifier et/ou redistribuer ce programme sous les conditions
// de la licence CeCILL| telle que diffusée par le CEA, le CNRS et l'INRIA 
// sur le site "http://www.cecill.info".

// Calcul d'une line array : indice de directivité

// Données d'entrées ----------------------------------------------------------
clear xs;
clear ys;
clear theta;
clear g;

// Nombre de source 
n = 11;

// Espacement entre source
de = 0.9;

// Largeur de la ligne
d = (n-1)*de ;

// Vitesse du son
C = 344 ;

// Fréquence
f = 100;

// Nombre de division angulaire
na = 180;

// Distance de calcul
rc = 9 ;

//Color : bleu=11, red=21, green=13, purple=24, yellow=32
cl = 13 

// Retard sur les elements de la ligne (ms)
// La taille de lr doit être la même que celle du nombre de source
tau = [2 ; 1 ; 0 ; 0 ; 1 ; 2 ];
//tau = [0 ; 1.59 ; 3.77 ; 3.77 ; 1.59 ; 0 ];
//tau = [0 ; 1.59/1.2 ; 3.77/1.2 ; 3.77/1.2 ; 1.59/1.2 ; 0 ];
//tau = [0 ; 1.59/1.2 ; 2.8 ; 2.8 ; 1.59/1.2 ; 0 ];
tau = [0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0 ;0 ;0];
tau = [2.5 ; 2 ; 1.5 ; 1 ; 0.5 ; 0 ; 0.5 ; 1 ; 1.5 ;2 ; 2.5];
//tau = [0 ; 0.5 ; 1 ; 1.5 ; 2 ; 0 ; 2 ; 1.5 ; 1 ;0.5 ; 0];


// Calculs --------------------------------------------------------------------

// Positions des sources
ym0 = d/2 ;
i=0;
for j=1:n
    i=i+1
    ys(i)=(j-1)*de-ym0;
    xs(i)=0;
end

// Pas d'angle
dtheta = %pi/2/na;

// Nombre d'onde
kn = 2*%pi*f/C;

// Longueur d'onde
lambda = C/f;

// Calcul pression dans l'axe
// Boucle sur les sources
paxe = 0;
for j=1:n

  // Calcul distance
  r = sqrt((xs(j)-rc)^2+ys(j)^2);
                
  // Contribution à la pression de la source i
  paxe = paxe + exp(-%i*kn*r)/r*exp(-%i*kn*C*tau(j)/1000);
end;

// Boucle sur l'angle
Q = 0;
for k = 1:na 
    
    // Calcul angle
    theta(k) = (k-0.5)*dtheta;
    
    // Calcul position microphone
    xm = rc*cos(theta(k));
    ym = rc*sin(theta(k));
         
    // Boucle sur les sources
    pts = 0;
    for j=1:n

        // Calcul distance
        r = sqrt((xs(j)-xm)^2+(ys(j)-ym)^2);
                
        // Contribution à la pression de la source i
        pts = pts + exp(-%i*kn*r)/r*exp(-%i*kn*C*tau(j)/1000);
    end;
    
    // Facteur de directivité
    g(k) = (abs(pts)/abs(paxe))^2 ;
    //g(k) = abs(pts/paxe)^2 ;
    
    // Facteur de directivité analytique
    gt(k) = cos(%pi*d*sin(theta(k))/lambda)^2;
    
    // Maj somme
    Q = Q + %pi/(2*na)*g(k)*cos(theta(k)); 
end

// Facteur de directivité (facteur 2 parce que rayonnement 2pi)
Q = 2/Q;
disp("Facteur de directivité numerique "+string(Q))
disp("Directivity index numerique "+string(10*log10(Q))+" dB")

// Polar plot
g = 10*log10(g)+50;
polarplot(theta,g,style=cl)
//legends(['40';'60';'100'],[11,21,13],opt="ur")

  
  


