omega_u=linspace(0,10,1000);
%b=linspace(10,100,1000);


F=10;
m=10;
k=100;
omega_0=sqrt(k/m);
b_c=2*sqrt(k*m)
n1=10;
b=(b_c-0.1)/n1;
b_sim=linspace(1,b_c-1,1000);
iplot=3;
N=n1;
%A=A_functie(F,m,omega_u,omega_0,b);


%plotten van de 10 grafieken met verschillende b
%-------------------------------------------------

Amaximas=zeros(n1,3); % 1ekolom = omega_u, 2ekolom = A_max, 2ekolom = b's

for i1=1:n1,
    A=A_functie(F,m,omega_u,omega_0,b*i1);
    Omega_uA=[omega_u',A'];
    [Y1,I1]=max(Omega_uA);	%beste rij voor elke kolom
    Amaximas(i1,:)=[Omega_uA(I1(2),:) b*i1];
    
    
    figure(1);
    plot(omega_u,A);
    zoom on;
    ylabel('Amplitude');
    xlabel('\omega_u (pulsatie van aandrijvende kracht)');
    title('Amplituderesonantiecurve');
    zoom on;
    hold on;
end;

saveas(gcf,'BijkomendeOpg1','tif');		%opslaan van figuur voor import in andere programma's
%saveas(gcf,'BijkomendeOpg1','fig');		%opslaan van figuur voor toekomstige her-import in Matlab

%maximas tekenen
%-------------------------------------------------
figure(1)
plot(Amaximas(:,1),Amaximas(:,2),'or');
ylabel('Amplitude en A_{max} (voor verschillende b)');
zoom on;

saveas(gcf,'BijkomendeOpg2a','tif');		%opslaan van figuur voor import in andere programma's
%saveas(gcf,'BijkomendeOpg2a','fig');		%opslaan van figuur voor toekomstige her-import in Matlab

%curve A_max in functie van b
%-------------------------------------------------

figure(2)
semilogy(Amaximas(:,3),Amaximas(:,2),'or');
xlabel('b (dempingssterkte)');
ylabel('A_{max}(b) punten');
zoom on;
hold on;

YY=[Amaximas(:,3),Amaximas(:,2)];
save 'Maximas10CurvesifvB.txt' YY -ascii -tabs

saveas(gcf,'BijkomendeOpg2aA_max(b)','tif');		%opslaan van figuur voor import in andere programma's
%saveas(gcf,'BijkomendeOpg2aA_ar','fig');		%opslaan van figuur voor toekomstige her-import in Matlab

%2.b simulatie volgens vgl. 6
figure(2)
semilogy(b_sim,A_ar_functie(F,b_sim,k,m));
xlabel('b (dempingssterkte)');
ylabel('A_{max}(b) punten eb A_{ar}(b) kromme');
zoom on;

saveas(gcf,'BijkomendeOpg2bA_max(b)enA_ar(b)','tif');		%opslaan van figuur voor import in andere programma's
%saveas(gcf,'BijkomendeOpg2aA_ar','fig');		%opslaan van figuur voor toekomstige her-import in Matlab

%2.c niet lineaire fit van A_max(b)-A_ar(b,k,m)
    iplot=2;
    n2=5000;
    n3=5000;
    ktestvec=linspace(k/10,k*5,n2)';
    mtestvec=linspace(m/10,m*5,n3)';
    chi2mat=zeros(n2,n3);
    
    for i1=1:n2,
        for i2=1:n3,  
          ktest=ktestvec(i1);
          mtest=mtestvec(i2);
          A_arTheoretisch=A_ar_functie(F,Amaximas(:,3),ktest,mtest); %Theoretische functie
          chi2=sum((abs(A_arTheoretisch-Amaximas(:,2))).^2);
          chi2mat(i1,i2)=chi2;		
            if iplot==1,
               figure(11)
               subplot(2,1,1);
               loglog(Iexp,Rloadexp,'o',Iexp,RloadTheoretisch,'*');
               zoom on;
               ylabel('R_{load}');
               xlabel('Stroom (Ampere)');
               set(gca,'FontSize',13)
               title(['par 1 attempt ',int2str(i1),' of ',int2str(n2),'  par 2 attempt ',int2str(i2),' of ',int2str(n3)],'FontSize',13);
            end;%end if iplot==1
        end;%end for i2
    end;%end for i1
  
%bepaling van beste fitparameters voor minimale chi2 

    [Y1,I1]=min(chi2mat);	%beste rij voor elke kolom
    [Y2,I2]=min(Y1);	%beste kolom van alle kolommen

    kbest=ktestvec(I1(I2));
    mbest=mtestvec(I2);

    figure(3)
    imagesc(mtestvec,ktestvec,log(chi2mat));
    colorbar('FontSize',12);
    xlabel('m','FontSize',13);
    ylabel('k','FontSize',13);
    title(['plot of log(chi2) values - best k: index ',int2str(I1(I2)),' value:',num2str(kbest), '  best m: index ',int2str(I2),' value: ',num2str(mbest)],'FontSize',13);
    hold on;
    plot(mtestvec(I2),ktestvec(I1(I2)),'*k')
    hold off;
    set(gca,'Ydir','Reverse');

saveas(gcf,'BijkomendeOpg2cNietLinFitInKenM','tif');		%opslaan van figuur voor import in andere programma's
%saveas(gcf,'BijkomendeOpg2aA_ar','fig');		%opslaan van figuur voor toekomstige her-import in Matlab

%bepaling van statistische onzekerheid op fitparameters (zonder rekening te houden met covariantie)
%m

    n12=50;
    mtestvec1=linspace(0.9*mbest,1.1*mbest,n12)';
    chi2vec=zeros(n12,1);

    for i2=1:n12,
      mtest=mtestvec1(i2);
      A_arTheoretisch=A_ar_functie(F,Amaximas(:,3),kbest,mtest); %Theoretische functie
      chi2=sum((abs(A_arTheoretisch-Amaximas(:,2))).^2);
      chi2vec(i2)=chi2;		
    end;%end for i2

    P= polyfit(mtestvec1,chi2vec,2);
    chi2parabool=polyval(P,mtestvec1);


    sigpar=sqrt(abs(-P(2)^2+4*P(1)*P(3)))/(2*P(1));
    parfitmin=-P(2)/(2*P(1));
    fprintf('parameterfitwaarde: %4.3e +/- %4.3e/sqrt(%i)\n',mbest,sigpar,N);
    astring=sprintf('parameterfitwaarde: %4.3e +/- %4.3e/sqrt(%i)\n',mbest,sigpar,N);

    figure(5);
    plot(mtestvec1,chi2vec,'o',mtestvec1,chi2parabool,'-k')
    zoom on;
    xlabel('m','FontSize',13);
    ylabel('chi2','FontSize',13);
    title(['chi2 afhankelijkheid van m bij beste waarde van k ',astring],'FontSize',13);
    set(gca,'FontSize',13)
    
%bepaling van statistische onzekerheid op fitparameters (zonder rekening te houden met covariantie)
%k

    n12=50;
    ktestvec1=linspace(0.9*kbest,1.1*kbest,n12)';
    chi2vec=zeros(n12,1);

    for i2=1:n12,
      ktest=ktestvec1(i2);
      A_arTheoretisch=A_ar_functie(F,Amaximas(:,3),ktest,mbest); %Theoretische functie
      chi2=sum((abs(A_arTheoretisch-Amaximas(:,2))).^2);
      chi2vec(i2)=chi2;		
    end;%end for i2

    P= polyfit(ktestvec1,chi2vec,2);
    chi2parabool=polyval(P,ktestvec1);


    sigpar=sqrt(abs(-P(2)^2+4*P(1)*P(3)))/(2*P(1));
    parfitmin=-P(2)/(2*P(1));
    fprintf('parameterfitwaarde: %4.3e +/- %4.3e/sqrt(%i)\n',kbest,sigpar,N);
    astring=sprintf('parameterfitwaarde: %4.3e +/- %4.3e/sqrt(%i)\n',kbest,sigpar,N);

    figure(7);
    plot(ktestvec1,chi2vec,'o',ktestvec1,chi2parabool,'-k')
    zoom on;
    xlabel('k','FontSize',13);
    ylabel('chi2','FontSize',13);
    title(['chi2 afhankelijkheid van k bij beste waarde van m ',astring],'FontSize',13);
    set(gca,'FontSize',13)
    
