%Laden van gegevens

    XX=load('BijkomendOpg1Data.txt');
    %Rloadexp=XX(:,1);		%1e kolom
    Iexp=XX(:,2)/1000;			%2de kolom
    Vexp=XX(:,3);			%3e kolom
    AdVexp=XX(:,4);			%4de kolom
    RdVexp=XX(:,5);			%5e kolom
    N=length(Iexp);

%Experimentele Kirchoff-functie
    Rloadexp=Vexp./(Iexp);

    iplot=2;
    n1=5000;
    n2=5000;
    Vemftestvec=linspace(5,20,n1)';
    Rinwtestvec=linspace(500,2500,n2)';
    chi2mat=zeros(n1,n2);

%bepaling van chi2 bij elke combinatie van fitparameters

    for i1=1:n1,
        for i2=1:n2,  
          Vemf=Vemftestvec(i1);
          Rinw=Rinwtestvec(i2);
          RloadTheoretisch=NortonRLoadFunction(Iexp,Vemf,Rinw); %Theoretische functie in fucntie van Iexp!
          chi2=sum((abs(RloadTheoretisch-Rloadexp)).^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(n1),'  par 2 attempt ',int2str(i2),' of ',int2str(n2)],'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

    Vemfbest=Vemftestvec(I1(I2));
    Rinwbest=Rinwtestvec(I2);

    figure(12)
    imagesc(Rinwtestvec,Vemftestvec,chi2mat);
    colorbar('FontSize',12);
    xlabel('R_{inwendig} (\Omega)','FontSize',13);
    ylabel('V_{emf} (V)','FontSize',13);
    title(['plot of chi2 values - best V_{emf}: index ',int2str(I1(I2)),' value:',num2str(Vemfbest), '  best R_{inw}: index ',int2str(I2),' value: ',num2str(Rinwbest)],'FontSize',13);
    hold on;
    plot(Rinwtestvec(I2),Vemftestvec(I1(I2)),'*k')
    hold off;
    set(gca,'Ydir','Reverse');

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

    n12=50;
    Rinwtestvec1=linspace(0.9*Rinwbest,1.1*Rinwbest,n12)';
    chi2vec=zeros(n12,1);

    for i2=1:n12,
      Vemf=Vemfbest;
      Rinw=Rinwtestvec1(i2);
      RloadTheoretisch=NortonRLoadFunction(Iexp,Vemf,Rinw); %Theoretische functie in fucntie van Iexp!
      chi2=sum((abs(RloadTheoretisch-Rloadexp)).^2);
      chi2vec(i2)=chi2;		
        if iplot==2,
           figure(11)
           loglog(Iexp,Rloadexp,'o',Iexp,RloadTheoretisch,'*');
           zoom on;
           ylabel('R_{load} (\Omega)','FontSize',13);
           xlabel('I (A)','FontSize',13);
           set(gca,'FontSize',13)
           hold on;%%%%%%%%%%%%%%%%%%%

           
           title(['par 2 attempt ',int2str(i2),' of ',int2str(n12)],'FontSize',13);
           set(gca,'FontSize',13)
           %hold on;%%%%%%%%%%%%%%%%%%%
        end;%end if iplot==1
    end;%end for i2

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


    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',Rinwbest,sigpar,N);
    astring=sprintf('parameterfitwaarde: %4.3e +/- %4.3e/sqrt(%i)\n',Rinwbest,sigpar,N);

    figure(14);
    plot(Rinwtestvec1,chi2vec,'o',Rinwtestvec1,chi2parabool,'-k')
    zoom on;
    xlabel('R_{inwendig}','FontSize',13);
    ylabel('chi2','FontSize',13);
    title(['chi2 afhankelijkheid van R_{inwendig} bij beste waarde van V_{uit} ',astring],'FontSize',13);
    set(gca,'FontSize',13)