%% Ricerca delle frequenze proprie e modi di vibrare di una trave
% come input la matrice H delle condizioni al contorno
% frequenze proprie energy harvester
% vincolo come cantilever, frequenze al variare della sporgenza
function gammaL_i=calcolo_frequenze(condizioniCC,EJ,m,L,nmax)

% valori del prodotto gamma*L per la ricerca preliminare degli zeri del
% det(H) - plottaggio della funzione det(H)
gammavet=pi*[0:0.01:nmax];

for k=1:length(gammavet)
    gammaL=gammavet(k);
    detvet(k)=det(condizioniCC(gammaL));
end

%% Visualizzazione ricerca zeri di det(H)
% ricerca punti di minimo locale della funzione log(det(H)+epsilon)
% epsilon=1e-20 sommato all'argomento di log per garantire esistenza
% logaritmo

derlog=log(abs(detvet(2:end)+1e-20))-log(abs(detvet(1:end-1)+1e-20));
segno=derlog(2:end).*derlog(1:end-1);
jj=find(segno < 0);
ii2=find(derlog(jj) < 0);
ii=jj(ii2)+1;
figure
plot(gammavet,log(abs(detvet+1e-20)),gammavet(ii),log(abs(detvet(ii)+1e-20)),'o');grid

%pause
nmodi=length(ii);
pause
tol=1.e-10;

%% ciclo sulle frequenze proprie e modi di vibrare
for imodo=1:nmodi
    %raffina la soluzione dell'autovalore
    gamma1=gammavet(ii(imodo)-5);
    det1=det(condizioniCC(gamma1));
    gamma2=gammavet(ii(imodo)+5);
    det2=det(condizioniCC(gamma2));
    dgamma=gamma2-gamma1;
    %     disp('inizio ricerca zero');
    %     pause
    while dgamma > tol
        gammas=0.5*(gamma1+gamma2);
        detdet=det(condizioniCC(gammas));
        if det1*detdet >=0
            gamma1=gammas;
            det1=detdet;
        else
            gamma2=gammas;
            det2=detdet;
        end
        dgamma=gamma2-gamma1;
    end
    
    % Calcola la frequenza propria
    gammaL_i(imodo)=gammas;
    freprop(imodo)=1/2/pi*(gammas./L)^2*sqrt(EJ/m);
    
    % Calcola la matrice H
    H=condizioniCC(gammas);
    n=size(H);
    ndom=round(n/4);
    % definisce le sottomatrici (n-1)x(n-1) della matrice H e calcola il det
    % normalizzato (in rapp): tutte le sottomatrici con det NON zero sono
    % ammissibili
    for iriga=n:-1:1
        for icolo=n:-1:1
            Hrid=[H(1:iriga-1,1:icolo-1) H(1:iriga-1,icolo+1:n);H(iriga+1:n,1:icolo-1) H(iriga+1:n,icolo+1:n)];
            rapp(iriga,icolo)=det(Hrid)/trace(H);
        end
    end
    [val,icolo]=max(max(rapp));
    [val,iriga]=max(rapp(:,icolo));
    %     iriga=2;
    Hrid=[H(1:iriga-1,1:icolo-1) H(1:iriga-1,icolo+1:n);H(iriga+1:n,1:icolo-1) H(iriga+1:n,icolo+1:n)];
    N=-[H(1:iriga-1,icolo);H(iriga+1:n,icolo)];
    x=Hrid\N;
    coef=[x(1:icolo-1,1);1;x(icolo:end,1)];
    %     coef=[0;1;0;0];
    disegna(ndom,imodo,gammas,coef,freprop(imodo))
    pause
end
end