clc
close all
clear all

%___POMIARY___
pom1=load('POM0CM0.txt'); %___pomiar z ustawieniem cewek symetrycznym
pom2=load('POM1CM0.txt'); %___pomiar z ustawieniem cewek z przesunięciem od ekranu o 1cm
% pom4=load('POM2CM0.txt'); %___pomiar z ustawieniem cewek z przesunięciem od ekranu o 2cm
% pom5=load('POMDP0.txt'); %___pomiar z ustawieniem cewek symetrycznym i dłuższym przewodem pomiarowym
% pom6=load('SKOS0.txt'); %___pomiar z ustawieniem cewek z przesunięciem skośnym
freq=logspace(1,7,600);
Rb=12.625;
abs_cor_sym=Rb.*10.^(pom1(:,1)/20);
ph_cor_sym=pom1(:,2);

abs_cor_asym=Rb.*10.^(pom2(:,1)/20);
ph_cor_asym=pom2(:,2);

for k=1:length(freq)   
  f1(k)=real(abs_cor_sym(k,1)*exp(i*ph_cor_sym(k,1)*pi/180));
  f2(k)=imag(abs_cor_sym(k,1)*exp(ph_cor_sym(k,1)*i*pi/180));
  f1a(k)=real(abs_cor_asym(k,1)*exp(i*ph_cor_asym(k,1)*pi/180));
  f2a(k)=imag(abs_cor_asym(k,1)*exp(ph_cor_asym(k,1)*i*pi/180));
  f(1,k)=f1(k)+i*f2(k);
  fa(1,k)=f1a(k)+i*f2a(k);
end

figure(1)
loglog(freq,f,freq,fa,'r');grid;

% liczba biegunów
NPoles = [0 10];

%% rational fitting and frequency response
fit = rationalfit(freq,f,'NPoles',NPoles,'TendsToZero',true);
resp = freqresp(fit,freq);
A = fit.A; C = fit.C; D = fit.D;

figure(2)
loglog(freq,abs(f),freq,abs(resp),'r');grid;