%% 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 gamma_i=calcolo_frequenze_massa_silent(condizioniCC,EJ,m,nmax,LL)

% valori del prodotto gamma 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)
    gamma=gammavet(k);
    detvet(k)=det(condizioniCC(gamma));
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

nmodi=length(ii);

tol=1.e-10;

%% ciclo sulle frequenze proprie e modi di vibrare (solo la prima!)
for imodo=1:1
    %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;
    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
    gamma_i(imodo)=gammas;
%     freprop(imodo)=1/2/pi*(gammas)^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
%     [~,icolo]=max(max(rapp));
%     [val,iriga]=max(rapp(:,icolo));
%     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)];
%     disegna_massa(ndom,imodo,gammas,coef,freprop(imodo),LL)
end

end