// Mesure des parametres d'un HP
// J. Fourcade (volucres.fr)

// Données d'entrées ----------------------------------------------------------
// Constantes physiques
rho = 1.18 ;    // Densité atmosphérique (kg/m^3)
c = 345 ;       // Vitesse du son (m/s)

// Directory où sont les fichiers de mesures
d = C:\Users\Pascal\Desktop;

// Mode d'ajustement
// ReAjust = 1 : Re libéré
// Rejust = 0 ; Re figé
ReAjust = 1;

// Paramatre HP
Re = 5.5 ;    // Resistance de la bobine (ohm)
Sd = 1225 ;   // Surface active de la membrane (cm^2)
dm = 0;        // Masse ajouté à la bobine mobile (gr). 
                // Si dm egale zero pas de calcul de Vas
                XXX ; 
// Nom du fichier de mesure du HP 
// Format : frequence, amplitude, ...
// Pas de commentaires dans le fichier
// Le fichier fic2 n'est pas utilisé si dm=0
fic1 = "XXX.txt";
fic2 = "Altec 416 8A 24851 dm.txt"

// Valeurs initiales des paramètres à identifier
// a modifier si probleme de convergence du moindre carré
Fs = 37 ;    // Frequence de résonnance (hz)
Qes = 0.32 ;    // Facteur de qualité électrique
Qms = 7 ;    // Facteur de qualité mécanique

// Calculs --------------------------------------------------------------------
// Conversion des paramètres en SI
dm=dm/1e3;
Sd=Sd/1e4;

// Variable de Laplace
s =poly(0,'s');

// Fonction d'ajustement (calcul des résidus)
function err=Ze(param, m) 
    Fs = param(1);
    Qes = param(2);
    Qms = param(3);
    if length(param) == 4 then
        Re = param(4);
    end
    
    err = Zm-abs(Re*repfreq(ZHP(Fs,Qes,Qms),frq))';
endfunction

// Ajustement sur le HP -------------------------------------------------------
// Lecture du fichier de mesure
HP=read (d+fic1, -1, 3);

// Liste des frequences et valeurs de l'impédance mesurées
frq=HP(:,1);
Zm=HP(:,2);
m=size(frq,1);

// Paramètres initiaux
if ReAjust == 1 then
    pinit=[Fs Qes Qms Re];
else
    pinit=[Fs Qes Qms];
end;
p=size(pinit,1);

// Ajustement par moindres carrés
[sol,v,info] = lsqrsolve(pinit,Ze,m)
disp("Code retour moindres carrés HP : "+string(info))

// Calcul parametres
Fs=sol(1);
ws=2*%pi*Fs;
Qes=sol(2);
Qms=sol(3);
if ReAjust == 1 then
    Re=sol(4);
end;

// Impédance calculée
Zc=abs(Re*repfreq(ZHP(Fs,Qes,Qms),frq))';
    
clf
subplot(121);
plot2d(frq,Zm,5,logflag="ln")
plot2d(frq,Zc,11,logflag="ln")
xtitle("Impédance (ohm)")
legends(['Mesurée';'Calculée'],[5,11],opt="ur")
subplot(122);
plot2d(frq,Zm-Zc,-5,logflag="ln")
xtitle("Résidus (ohm)")

// Ajustement du HP avec la masse additionnelle -------------------------------
if dm<>0 then
    // Lecture du fichier de mesure
    HP=read (d+fic2, -1, 3);

    // Liste des frequences et valeurs de l'impédance mesurées
    frq=HP(:,1);
    Zm=HP(:,2);
    m=size(frq,1);
    
    // Paramètres initiaux
    pinit=[Fs Qes Qms]

    // Ajustement par moindres carrés
    [sol,v,info] = lsqrsolve(pinit,Ze,m)
    disp("Code retour moindres carrés HP avec masse : "+string(info))
    Fsd=sol(1)
    wsd=2*%pi*Fsd;
    Qesd=sol(2);
    Qmsd=sol(3);
end

// Impressions résultats
disp("Estimation des paramètres du HP : ")
disp("Fs  = "+string(Fs)+" hz")
disp("Qes = "+string(Qes))
disp("Qms = "+string(Qms))
disp("Re  = "+string(Re)+" ohm")

if dm<>0 then
    // Calcul Vas et Mms, Cms, Rms, Bl
    Mms=dm/(Fs*Qesd/(Fsd*Qes)-1);
    Cms=1/(Mms*ws^2);
    Vas=rho*c^2*Sd^2*Cms;
    Rms=1/(ws*Cms*Qms);
    Bl=sqrt(Re/(ws*Cms*Qes));
    disp("Mms  = "+string(Mms*1e3)+" g")
    disp("Cms  = "+string(Cms)+" m/N")
    disp("Rms  = "+string(Rms)+" kg/s")
    disp("Bl   = "+string(Bl)+" Tm")
    disp("Vas  = "+string(Vas*1e3)+" l")
    
     // Impression des parametres du HP avec masse additionelle
    disp("Estimation des parametres du HP avec masse additionnelle : ")
    disp("Fsd  = "+string(Fsd)+" hz")
    disp("Fs/Fsd   = "+string(Fs/Fsd))
    disp("Qesd = "+string(Qesd))
    disp("Qmsd = "+string(Qmsd))

    // Calcul Cmsd, Rmsd
    Mmsd=Mms+dm;
    Cmsd=1/(Mmsd*wsd^2);
    Rmsd=1/(wsd*Cmsd*Qmsd);
    disp("Cmsd  = "+string(Cmsd)+" m/N")
    disp("Rmsd  = "+string(Rmsd)+" kg/s")
end;