% clear;
% clc;

fprintf('Bode diagram of the RC-circuit\n\n');
fprintf('Figure 1: Bode diagram theory+asymptote + average measurement\n')
fprintf('Figure 2: Bode diagram average measurement+SD + theory\n');
fprintf('Figure 3: Bode diagram average measurement+SD\n\n');
fprintf('N.B. SD=0 is not displayed in the figures!\n\n');

% fprintf('\nWatch out! This script searches the file in "D:\"\n\n');
% 
% f_name=input('Filename: ','s');
% f_name1=sprintf('%s%s',f_name,'_theorie');
% f_name2=sprintf('%s%s',f_name,'_meting');
% 
% wd=cd;
% cd 'D:\'
% 
% [th]=dlmread(f_name1,' ',0,0);
% [met]=dlmread(f_name2,' ',0,0);
% 
% cd(wd);



f_th=th(1,:);
amp_th=th(2,:);
fase_th=th(3,:);
asym=th(4,:);
if mean(met(1,:))==0;
    nr_met=0;
end;

if mean(met(1,:))~=0;
    size_met=size(met);
    nr_met=(size_met(1)-1)/2;
    nr_punt=size_met(2);
    f_m=met(1,:);
    
    for i=1:nr_met;
        for j=1:nr_punt;
            amp_met(i,j)=met(i+1,j);
            fase_met(i,j)=met(i+1+nr_met,j);
        end;
    end;
    if nr_met~=1;
        for i=1:nr_punt;
            av_amp(i)=mean(amp_met(:,i));
            av_log_amp(i)=20*log10(av_amp(i));
            av_fase(i)=mean(fase_met(:,i));
        end;
        sd_amp=std(amp_met,0,1);
        for i=1:nr_punt;
            if sd_amp(i)==0;
                sd_log_amp(i)=0;
            elseif av_amp(i)==1;
                av_amp(i)=1.0001;
            else
                sd_log_amp(i)=2/log(av_amp(i))*sd_amp(i);
            end;
        end;
        sd_fase=std(fase_met,0,1);
    else 
        av_log_amp=20*log10(amp_met);
        av_fase=fase_met;
        sd_log_amp=zeros(1,nr_punt);
        sd_fase=zeros(1,nr_punt);
    end;
end;

figure(1);
set(clf,'PaperType','A4');
set(clf,'Color',[1,1,1]);
%set(clf,'PaperPosition',[1.13386 3.97146 5.5 4.125]);
subplot(2,1,1);
if mean(met(1,:))~=0;
    graf=semilogx(f_th,amp_th,'r',f_th,asym,'g',f_m,av_log_amp,'.-b','LineWidth',1);
    legend([graf],'theory','asymp.','measurement',3);
else 
    graf=semilogx(f_th,amp_th,'r',f_th,asym,'g','LineWidth',1);
    legend([graf],'theory','asymp.',3);
end;
ylabel('Amplitude [dB]','Fontsize',10)

% grid bepaling
startf=floor(log10(f_th(1)));
sizef=size(f_th);
endf=ceil(log10(f_th(sizef(2))));
xf(1)=10^startf;
i=startf+1;
count=2;
while i<=endf;
    xf(count)=10^(i);
    count=count+1;
    i=i+1;
end;
set(gca, 'XTick',xf,'Fontsize',9);

if mean(met(1,:))~=0;
    min_amp_th=floor(min(amp_th));
    min_amp_met=floor(min(av_log_amp-abs(sd_log_amp)));
    min_amp=min(min_amp_th,min_amp_met);
    max_amp_th=ceil(max(amp_th));
    max_amp_met=ceil(max(av_log_amp+abs(sd_log_amp)));
    max_amp=max(max_amp_th,max_amp_met)+2;
else
    min_amp=floor(min(amp_th));
    max_amp=ceil(max(amp_th))+2;
end;
    
axis([10^startf,10^endf,min_amp,max_amp]);
grid on;
title('Bode diagram RC-circuit','Fontsize',10,'Fontweight','bold');

subplot(2,1,2);
if mean(met(1,:))~=0;
    graf=semilogx(f_th,fase_th,'r',f_m,av_fase,'.-b','LineWidth',1);
    legend([graf],'theory','measurement',3);
else
    graf=semilogx(f_th,fase_th,'LineWidth',1);
    legend([graf],'theory',3);
end;
xlabel('Frequency [Hz]','Fontsize',10);
ylabel('Phase [rad]','Fontsize',10);
set(gca, 'XTick',xf,'Fontsize',9);

if mean(met(1,:))~=0;
    min_fase_th=floor(min(fase_th));
    min_fase_met=floor(min(av_fase-abs(sd_fase)));
    min_fase=min(min_fase_th,min_fase_met);
    max_fase_th=ceil(max(fase_th));
    max_fase_met=ceil(max(av_fase+abs(sd_fase)));
    max_fase=max(max_fase_th,max_fase_met);
else
    min_fase=floor(min(fase_th));
    max_fase=ceil(max(fase_th));
end;  
    
axis([10^startf,10^endf,min_fase,max_fase]);
grid on;

if nr_met~=0;
    figure(2);
    set(clf,'PaperType','A4');
    set(clf,'Color',[1,1,1]);
%    set(clf,'PaperPosition',[1.13386 3.97146 5.5 4.125]);
    subplot(2,1,1);
    graf=semilogx(f_m,av_log_amp,'.-b',f_th,amp_th,'-r','LineWidth',1);
    legend([graf],'measurement','theory',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_log_amp(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            xl1=[10^(log(f_m(i))/log(10)-width_sd),10^(log(f_m(i))/log(10)+width_sd)];
            yl1=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)+sd_log_amp(i)];
            yl2=[av_log_amp(i)-sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    ylabel('Amplitude [dB]','Fontsize',10);

    % grid bepaling
    startfm=floor(log10(f_m(1)));
    sizefm=size(f_m);
    endfm=ceil(log10(f_m(sizefm(2))));
%    xfm(1)=10^startfm;
    xfm(1)=f_m(1);
    i=startfm+1;
    countm=2;
    flag=0;
    while flag~=1;
        if xfm(countm-1)<f_m(nr_punt);
            xfm(countm)=10^(i);
            countm=countm+1;
            i=i+1;
        else flag=1;
        end;
    end;
    xfm(countm-1)=f_m(nr_punt);
    %    set(gca,'XTick',xfm,'Fontsize',9);
    
    fmin_th=f_th(1);
    count_min=1;
    while fmin_th<f_m(1);
        count_min=count_min+1;
        fmin_th=f_th(count_min);
    end;

    size_fth=size(f_th);
    size_fm=size(f_m);
    fmax_th=f_th(size_fth(2));
    count_max=size_fth(2);
    while fmax_th>f_m(size_fm(2));
        count_max=count_max-1;
        fmax_th=f_th(count_max);
    end;

    amp_th2=amp_th(count_min:count_max);
    min_amp_met=floor(min(av_log_amp-abs(sd_log_amp)));
    min_amp_th2=floor(min(amp_th2));
    max_amp_met=ceil(max(av_log_amp+abs(sd_log_amp)));
    max_amp_th2=ceil(max(amp_th2));
    min_amp_met2=floor(min(min_amp_met,min_amp_th2))-2;
    max_amp_met2=ceil(max(max_amp_met,max_amp_th2))+2;

    set(gca,'Fontsize',9);
    axis([xfm(1),xfm(countm-1),min_amp_met2,max_amp_met2]);
    grid on;
    title('Bode diagram RC-circuit','Fontsize',10,'Fontweight','bold');

    subplot(2,1,2);
    
    graf=semilogx(f_m,av_fase,'.-b',f_th,fase_th,'-r','LineWidth',1);
    legend([graf],'measurement','theory',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_fase(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_fase(i)+sd_fase(i),av_fase(i)-sd_fase(i)];
            xl1=[10^(log10(f_m(i))-width_sd),10^(log10(f_m(i))+width_sd)];
            yl1=[av_fase(i)+sd_fase(i),av_fase(i)+sd_fase(i)];
            yl2=[av_fase(i)-sd_fase(i),av_fase(i)-sd_fase(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    
    xlabel('Frequency [Hz]','Fontsize',10);
    ylabel('Phase [rad]','Fontsize',10);
%    set(gca, 'XTick',xfm,'Fontsize',9);
 
    min_fase_met=floor(min(av_fase-abs(sd_fase)));
    max_fase_met=ceil(max(av_fase+abs(sd_fase)));
    fase_th2=fase_th(count_min:count_max);
    min_fase_th2=floor(min(fase_th2));
    max_fase_th2=ceil(max(fase_th2));
    min_fase_met2=floor(min(min_fase_met,min_fase_th2))-1;
    max_fase_met2=ceil(max(max_fase_met,max_fase_th2));
    
    set(gca,'Fontsize',9);
    axis([f_m(1),f_m(nr_punt),min_fase_met2,max_fase_met2]);
    grid on;
end;

if nr_met~=0;
    figure(3);
    set(clf,'PaperType','A4');
    set(clf,'Color',[1,1,1]);
%    set(clf,'PaperPosition',[1.13386 3.97146 5.5 4.125]);
    subplot(2,1,1);
    graf=semilogx(f_m,av_log_amp,'.-b','LineWidth',1);
%    legend([graf],'meting',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_log_amp(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            xl1=[10^(log(f_m(i))/log(10)-width_sd),10^(log(f_m(i))/log(10)+width_sd)];
            yl1=[av_log_amp(i)+sd_log_amp(i),av_log_amp(i)+sd_log_amp(i)];
            yl2=[av_log_amp(i)-sd_log_amp(i),av_log_amp(i)-sd_log_amp(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    ylabel('Amplitude [dB]','Fontsize',10);

    % grid bepaling
    startfm=floor(log10(f_m(1)));
    sizefm=size(f_m);
    endfm=ceil(log10(f_m(sizefm(2))));
%    xfm(1)=10^startfm;
    xfm(1)=f_m(1);
    i=startfm+1;
    countm=2;
    flag=0;
    while flag~=1;
        if xfm(countm-1)<f_m(nr_punt);
            xfm(countm)=10^(i);
            countm=countm+1;
            i=i+1;
        else flag=1;
        end;
    end;
    xfm(countm-1)=f_m(nr_punt);
%    set(gca,'XTick',xfm,'Fontsize',9);
    
    min_amp_met=floor(min(av_log_amp-abs(sd_log_amp)))-2;
    max_amp_met=ceil(max(av_log_amp+abs(sd_log_amp)))+2;

    set(gca,'Fontsize',9);
    axis([xfm(1),xfm(countm-1),min_amp_met,max_amp_met]);
    grid on;
    title('Bode diagram RC-circuit (measurement)','Fontsize',10,'Fontweight','bold');

    subplot(2,1,2);
    
    graf=semilogx(f_m,av_fase,'.-b','LineWidth',1);
%    legend([graf],'meting',3);
    hold on;
    width_sd=(log10(f_m(2))-log10(f_m(1)))/3;
    for i=1:nr_punt;
        if sd_fase(i)~=0;
            x=[f_m(i),f_m(i)];
            y=[av_fase(i)+sd_fase(i),av_fase(i)-sd_fase(i)];
            xl1=[10^(log10(f_m(i))-width_sd),10^(log10(f_m(i))+width_sd)];
            yl1=[av_fase(i)+sd_fase(i),av_fase(i)+sd_fase(i)];
            yl2=[av_fase(i)-sd_fase(i),av_fase(i)-sd_fase(i)];
            plot(x,y,'-b');
            plot(xl1,yl1,'-b');
            plot(xl1,yl2,'-b');        
        end;
    end;
    hold off;
    
    xlabel('Frequency [Hz]','Fontsize',10);
    ylabel('Phase [rad]','Fontsize',10);
%    set(gca, 'XTick',xfm,'Fontsize',9);
 
    min_fase_met=floor(min(av_fase-abs(sd_fase)));
    max_fase_met=ceil(max(av_fase+abs(sd_fase)))+2;
 
    set(gca,'Fontsize',9);
    axis([f_m(1),f_m(nr_punt),min_fase,max_fase]);
    grid on;
end;