clear all
clc
global nummodes num_par
options.Resize='on';
%%% read measured data from experiments
testcasenumber=menu('Test case:','Case 1','Case 2');
switch testcasenumber
    case 1
        load measured_case1
        defaultanswer1={'10','6'};
        defaultanswer2={'4'};
        defaultanswer3={'0.004','6.4000e-007'};
        prompt1 = {'Enter number of samples:','Enter number of modes for updating:'};
        answer1 = inputdlg(prompt1,'',1,defaultanswer1,options);
        answer11=str2double(answer1);
        numsamples=answer11(1);nummodes=answer11(2);
        Varoutexp=outputdat(1:numsamples,1:nummodes);
        %Hz ---> (rad/s)^2
        Varoutexp=(2*pi*Varoutexp).^2;
        % mean vector and covariance matrix of measured data
        expFreq=mean (Varoutexp)';
        Vartest=cov (Varoutexp);
        % initial values of mean and covariance of updating parameters
        prompt2 = {'Enter number of updating parameters:'};
        answer2 = inputdlg(prompt2,'',1,defaultanswer2,options);
        num_par=str2double(answer2);
        Par=zeros(num_par,1);
        Varp=zeros(num_par,num_par);
        for i=1:num_par
            answer3 = inputdlg({['P' num2str(i)],['var P' num2str(i)]},'',1,defaultanswer3,options);
            answer33 = str2double(answer3);
            Par(i,1)= answer33(1);
            Varp(i,i,1)=answer33(2);
            clear answer3 answer33
        end
        std(:,1)=sqrt(diag(Varp(:,:,1)));
        %%%%%%%%%%%%%% iteration for stochastic model updating  %%%%%%%%%%%%%%%%
        % maximum number of iteration
        iteration=50;

        for itr=1:iteration
            [egvec egval M K]=fem (Par(:,itr));
            egvala=diag(egval);
            egvala=egvala(4:end);
            egveca=egvec(:,4:end);
            [Jac]=sens (egveca,egvala,Par(:,itr));
            CZZ=Jac*Varp(:,:,itr)*Jac';
            CTZ=Varp(:,:,itr)*Jac';
            anaFreq=egvala(1:6);
            %weighting value obtained from L_curve
            lambda=1e10;
            W1=diag([1/anaFreq(1),1/anaFreq(2),1/anaFreq(3),1/anaFreq(4),1/anaFreq(5),1/anaFreq(6)]);
            Tbar=inv(Jac'*W1*Jac+lambda*eye(4))*Jac'*W1;
            Par(:,itr+1)=Par(:,itr)+Tbar*(expFreq-anaFreq);
            Varp(:,:,itr+1)=Varp(:,:,itr)-CTZ*Tbar'-Tbar*CTZ'+Tbar*(CZZ+Vartest)*Tbar';
            std(:,itr+1)=sqrt(diag(Varp(:,:,itr+1)));
            if norm(Par(:,itr+1)-Par(:,itr))<1e-10 & norm(std(:,itr+1)-std(:,itr))<1e-10
                iteration=itr;
                break
            end

        end
        %ploting convergence figures for case 1
        figure (1)
        plot((0:iteration), Par(1,:),(0:iteration), Par(2,:),(0:iteration), Par(3,:),(0:iteration), Par(4,:));
        xlabel('Iterations');ylabel('Mean value of parameters');grid on;hold on;
        legend('t_1','t_2','t_3','t_4');

        figure (2)
        plot((0:iteration), std(1,:)./Par(1,:),(0:iteration), std(2,:)./Par(2,:),(0:iteration), std(3,:)./Par(3,:),(0:iteration), std(4,:)./Par(4,:));
        xlabel('Iterations');ylabel('COV of parameters');grid on;hold on;
        legend('COV (t_1)','COV (t_2)','COV (t_3)','COV (t_4)')
    case 2
        cd case2
        load measured_case2
        load options
        optimtool

end

